Simulation study on parametric dependence of whistler-mode hiss generation in the plasmasphere

We conduct electromagnetic particle simulations to examine the applicability of nonlinear wave growth theory to the generation process of plasmaspheric hiss. We firstly vary the gradient of the background magnetic field from a realistic model to a rather steep gradient model. Under such variation, the threshold amplitude in the nonlinear theory increases quickly and the overlap between threshold and optimum amplitude disappears correspondingly, the nonlinear process is suppressed. In the simulations, as we enlarge the gradient coefficient of the background magnetic field, waves generated near the equator do not grow through propagation. By examining the range of suitable values of inhomogeneity factor S (i.e., |S|<2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|S|<2$$\end{document}), we find the generation of wave packets is limited to the equatorial region when the background field is steep, showing a good agreement with what is indicated by critical distance in the theory. We then check the dependence of generation of hiss emissions on different hot electron densities. Since the overlap between threshold and optimum amplitude vanishes, the nonlinear process is weakened when hot electron density becomes smaller. In the simulation results, we find similar wave structures in all density cases, yet with different magnitudes. The existence of suitable S values implies that the nonlinear process occurs even at a low level of hot electron density. However, by examining JE\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_E$$\end{document} that closely relates to the wave growth, we find energy conveyed from particles to waves is much limited in small density cases. Therefore, the nonlinear process is suppressed when hot electron density is small, which agrees with the theoretical analysis.


Introduction
Plasmaspheric hiss is a kind of whistler-mode electromagnetic wave with a typical frequency range from ∼ 100 Hz to several kHz and an amplitude range from a few pT to hundreds of pT. It broadly occurs in the Earth's inner magnetosphere and significantly influences the dynamics of energetic electrons. In the past few decades, many observations were conducted (e.g., Smith et al. 1974;Ondoh et al. 1983;Singh et al. 1999;Meredith et al. 2004;Li et al. 2015;Liu et al. 2020) and several generation mechanisms of plasmaspheric hiss were proposed. Although most of those generation mechanisms are based on the linear growth theory (Kennel and Petschek 1966), there still exist debates over whether seeds of the instability derive inside or outside of the plasmasphere. For instance, Chen et al. (2014) put forward that seeds of whistlermode hiss emissions may come from internal thermal fluctuations. Draganov et al. (1992) suggested that an external source from lighting may cause the generation of whistler-mode hiss. Santolík et al. (2006) proposed that plasmaspheric hiss may originate from whistler-mode chorus emissions generated near the geomagnetic equator outside the plasmasphere. Based on both spacecraft observations along with ray tracing study, Bortnik et al. (2008) suggested that chorus waves may propagate into the plasmasphere and evolve into hiss. However, through observations from the Van Allen Probes EMFISIS (Electric and Magnetic Field Instrument Suite and Integrated Science) instrument combined with simulations, Hartley et al. (2019) found that only a very small fraction of chorus satisfies the required conditions to enter the plasmasphere.
Based on high-resolution waveform data from EMFISIS, Summers et al. (2014) first identified and examined embedded fine structures of plasmaspheric hiss, which were previously considered as structureless and incoherent. The result showed that hiss emissions possess short rising and falling tone structures similar to the sub-packets of chorus. Omura et al. (2015) showed that the generation of plasmaspheric hiss shares the similar mechanism that was originally used in explaining the generation of whistler-mode chorus (Omura et al. 2008). Nakamura et al. (2016) found a good agreement between observed frequency sweep rates for rising tone elements and the nonlinear theory while applying the theory to hiss emissions. Hikishima et al. (2020) conducted particle simulations and succeeded in reproducing hiss emissions following the nonlinear theory.
To further check the applicability of nonlinear wave growth theory to the generation of hiss emissions, we conduct self-consistent electromagnetic particle simulations with various parameters and compare the results with theoretical analyses. In section "Simulation method", we describe the simulation method implemented in detail. In section "Dependence on the gradient of background magnetic field", we discuss the dependence of generation of hiss emissions on the gradient of the background magnetic field. In section "Dependence on hot electron density", we analyze how the variation of hot electron density influences the generation of plasmaspheric hiss. In section "Summary and discussion", we give the summary and discussion.

Simulation method
A one-dimensional simulation model based on KEMPO (Kyoto university ElectroMagnetic Particle cOde) code is implemented (Matsumoto and Omura 1985;Hikishima et al. 2009). In the model, the electrostatic field parallel to the background magnetic field is omitted, while the magnetic field is given by B h = B 0 (1 + ah 2 ) as an approximation of the dipole field line near the equator. Here, B 0 is the field strength at the equator, a and h are the gradient coefficient of the parabolic magnetic field and the distance away from the equator, respectively. In the simulations, we set B 0 as unity and B h calculated at the boundary of the field line as a parameter to derive a. Normally, the gradient coefficient a = 1.86 × 10 −7 (� e0 /c) 2 , where e0 is the equatorial gyrofrequency of electron and c is speed of light, approximating the real dipole field of the earth at L = 4 . In this paper, we specifically assume a = 1.86 × 10 −7 (� e0 /c) 2 as Case 1, a = 3.74 × 10 −6 (� e0 /c) 2 as Case 2, and a = 3.74 × 10 −5 (� e0 /c) 2 as Case 3.
An open boundary condition along with damping regions at both boundaries is implemented to handle outgoing waves (Umeda et al. 2001), where the range of damping regions at each side is 6.25% system length. Inside the simulation region, two species of cold and hot electrons are initialized. The cold component is used to create a background plasma environment supporting wave propagation and the hot one interacts with waves through cyclotron resonance. For Case 1, the plasma frequencies of cold and hot species are specified as ω pe = 15� e0 and ω ph = 0.3� e0 , respectively. In the simulations, we keep the ratio between charge and mass of both species constant, using the relation of the plasma frequency ω p and density n to obtain values of particle charge, where ǫ 0 is vacuum permittivity, m is the mass of a super-particle, and n = N p /h s . Here, N p is the number of super-particles representing cold or hot electrons ( N p cold /N p hot in Table 1) and h s is the system length.
We allot the average numbers of cold and hot particles at each grid to an approximate value of 1024 and 2048, respectively. In the initialization, both species obey a subtracted Maxwellian distribution in the entire system region (Omura 2021). The equation of the subtracted Maxwellian distribution is written as (2): Here, u and u ⊥ are the parallel and perpendicular momenta of particles, U t and U t⊥ are the parallel and perpendicular thermal momenta of particles, respectively. Parameters ρ (0 ≤ ρ ≤ 1) and β (0 < β < 1) specify the relative height and width of the momentum distribution subtracted from a Maxwellian distribution, respectively. In the simulations, to keep the consistency with Thermal velocity of hot electrons at the equator: v th� , v th⊥ 0.24c, 0.28c the previous study (Hikishima et al. 2020), we assume ρ = 1 , and β = 0.3 . For the hot species, the thermal velocities are set as v th� = 0.24c and v th⊥ = 0.28c , creating an initial temperature anisotropy where v and v ⊥ are the parallel and perpendicular velocities of particles.
We adopt the initialization method that generates velocities of hot particles at the equator first and then allocate their spatial positions to the whole simulation region following the first adiabatic invariant (Nogi et al. 2020). Specifically, we calculate mirror points of particles according to their parallel and perpendicular velocities at the equator, then distribute them over the regions between mirror points. For particles with mirror points out of the simulation range, the method sets a boundary condition to ensure that all of them are located inside the region. In this approach, hot particles form the subtracted Maxwellian distribution and undergo bounce motions at the same time. In the simulations, particles that reach the boundary will be reflected with a randomized gyrophase of perpendicular velocity. Figure 1 shows the initial velocity distribution of hot electrons and its variation at t = 20480[ −1 e0 ] in Case 1. In the simulations, we count particles that fall into the loss cone as missing. From Fig. 1b, we see the value of the velocity distribution at the low pitch angle part becomes smaller, representing that some particles are scattered into the loss cone. The decrement at the high pitch angle part indicates that the free (c) Case 3 (b) Case 2 (a) Case 1

Dependence on the gradient of background magnetic field
We change the gradient of the background magnetic field in both theoretical analyses and simulations to check the correlation between them. From the nonlinear wave growth theory (Omura 2021), the ideal nonlinear process happens at the inhomogeneity factor S ≈ −0.4 , as shown in (3), where s 0 , s 1 , and s 2 are three parameters defined in (38) to (40) in Omura (2021), and ω , w , and e are the frequency, amplitude of waves, and the gyrofrequency of electrons, respectively. From (3), the whole nonlinear wave growth process can be divided into two stages. The first stage is due to the frequency variation ∂ω/∂t , and the second stage is due to the gradient ∂� e /∂h . Specifically, the first stage is dominant near the equator where the gradient is 0, and wave packets are formed and growing at this stage. After the formation of wave packets, their frequencies will not change through the propagation. When the wave packets propagate from the equator to higher latitudes, the second stage plays a more important role in keeping S ≈ −0.4 . The second stage is also called the convective growth stage. Case 2 a=3.74 10 -6 ( e0 /c) 2 Case 1 a=1.87 10 -7 ( e0 /c) 2 half system length Assuming the gradient term to be 0 and letting S = −0.4 , we derive the optimum amplitude that represents an amplitude that gives the maximum nonlinear wave growth condition. In the first stage of the nonlinear process, waves cannot grow much greater than the optimum amplitude. Similarly, if the frequency variation term is omitted and letting S = −0.4 , we obtain the threshold amplitude, meaning the minimum amplitude that a wave packet can grow locally as an absolute instability and enter an ideal convective growth stage. We refer to the equations of the normalized optimum and threshold amplitudes in Omura (2021) and write them in (4) and (5): and ã = ac 2 /� 2 e0 are normalized optimum amplitude, threshold amplitude, phase velocity, group velocity, average perpendicular (4) ] , gradient coefficient a in three cases are a 1.87 × 10 −7 (� e0 /c) 2 ; b 3.74 × 10 −6 (� e0 /c) 2 , and c 3.74 × 10 −5 (� e0 /c) 2 , respectively. Typical wave structures with amplitudes over 10 −4.5 [B 0 ] are outlined momentum, parallel thermal momentum, wave frequency, hot plasma frequency, resonant velocity and gradient coefficient of background field, respectively. Parameters ξ and χ are defined by ξ 2 = ω(� e0 − ω)/ω 2 pe and χ 2 = 1/(1 + ξ 2 ) , respectively, specifying the dispersion relation. The Lorentz factor is given by where v and v ⊥ are the parallel and perpendicular velocity of particles. In addition, Q is the depth of an electron hole due to the depletion of trapped resonant electrons in the velocity phase space, and τ represents the ratio between the nonlinear transition time and the nonlinear trapping time. In this Here, column a corresponds to Case 1 ( a = 1.87 × 10 −7 (� e0 /c) 2 ), column b corresponds to Case 2 ( a = 3.74 × 10 −6 (� e0 /c) 2 ), column c corresponds to Case 3 ( a = 3.74 × 10 −5 (� e0 /c) 2 ). The red dashed lines represent critical distances in Case 2 and Case 3 study, we assume Q = 0.1 and τ = 0.5 to keep the consistency with the previous research.
From (4) and (5), the gradient coefficient a only affects the threshold amplitude, which is also illustrated in Fig. 2 that the overlap between threshold and optimum amplitudes vanishes in large gradient cases. The theoretical result demonstrates that in Case 2 and Case 3, wave amplitudes after the first stage of the nonlinear wave growth are not large enough to enter an ideal convective growth stage, therefore, S = −0.4 may not be satisfied at high latitudes and the nonlinear convective growth is suppressed there. Figure 3 shows the spatial and temporal evolution of forward-propagating waves under different gradients of the magnetic field in the simulations. As the gradient increases, we find that the wave amplitudes at high latitudes are gradually decreasing especially in the first half period of the simulation run ( 0 ∼ 2 × 10 4 [ −1 e0 ] ), declaring a good agreement with the theory. It is worth noting that the group velocity of wave packets slightly increases for the large gradient case, leading to a curved form of wave path, as shown in Fig. 3c. This is because, in Case 3, we specified the background gradient that is 200 times greater than that in Fig. 3a (i.e., Case 1).
In the theory, we define the transition point of the two stages as a critical distance h c , at which the frequency variation term and the gradient term in (  are equal. Inside the critical distance, the first stage is dominant, and we can assume wave packets are generated there. In Fig. 4, focusing on the typical frequency of hiss at around 0.05 e0 , we find h c in Case 1 is greater than half of the system length, meaning wave packets can be formed inside the whole simulation region. For Case 2 and Case 3, however, their critical distances are smaller, indicating that the generation of wave packets is limited to the equatorial region.
In Fig. 5, we obtain dynamic spectra of forward-propagating waves at various locations of the simulation region for all three gradient cases, and outline structures with wave amplitudes over 10 −4.5 [B 0 ] . In Case 1, some typical wave packets are generated from h = −100[c −1 e0 ] .
In Case 2, there are some typical structures at h = −50[c −1 e0 ] . However, we can only find obvious wave packets from h = 0[c −1 e0 ] in Case 3. The result is consistent with the theoretical critical distance.
From h = −100[c −1 e0 ] to −50[c −1 e0 ] in Case 2 and Case 3 of Fig. 5, wave amplitudes are growing. This is caused by the linear growth process. In Fig. 6, we examine the effectiveness of the linear and nonlinear process on the growth of wave amplitudes. In Fig. 6a, after rapid growth at the initial adjustment stage due to thermal fluctuations of hot electrons, the averaged wave amplitudes in the whole simulation region indicate that Case 2 and Case 3 have similar growth rates to the linear growth rate at the initial stage, while Case 1 undergoes a stronger nonlinear process, leading to a higher overall amplitude.  We refer to the equation (92) in Omura (2021) and calculate the nonlinear growth rates at h = −100, −50 , and 0[c −1 e0 ] in Case 1, as shown in Fig. 6b. As the wave amplitude becomes larger, the nonlinear growth rate will decrease correspondingly. This is illustrated in Fig. 6b that as wave packets approach the equator, the nonlinear growth rate continues decreasing and reaches the level of the linear growth rate.
In Fig. 7, we implement a band-pass filter of 0.04 ∼ 0.06 e0 to forward-propagating waves within the time range of 0 ∼ 2 × 10 4 [ −1 e0 ] and calculate their instantaneous frequencies and inhomogeneity factor S. From wave amplitudes and frequencies in the first and second row, the area where wave packets are generated becomes much closer to the equatorial region from Case 1 to Case 3, showing a good agreement with the theory. From Tobita and Omura (2018), the nonlinear process can exist even when the inhomogeneity factor S is down to −2 . To show the variation of S structures, we adjust the colorbar range to 0 to 5 in the bottom panels of Fig. 7. We notice that S values which correspond to the nonlinear process (i.e., |S| < 2 ) widely locate in the whole simulation region in Case 1, yet are much limited to the equator in Case 2 and Case 3. This is consistent with the critical distance in each case, as indicated by the red dashed lines in the figure. The result implies that gradients in Case 2 and Case 3 are too large to satisfy the condition of nonlinear process at high latitudes, therefore not only their first stage of the nonlinear process has a very  Fig. 11 limited range, also their convective growth stage is strictly suppressed. Wave packets at high latitudes in Case 2 and Case 3 only undergo a weak linear growth process.

Dependence on hot electron density
We vary the hot electron density in both theoretical analyses and simulations. The hot plasma frequency ω ph is used to represent the hot electron density. Particularly, we assume ω ph = 0.15� e0 as Case 4 and ω ph = 0.075� e0 as Case 5 to make comparisons with Case 1. In the Fig. 11 Spatial and temporal evolution of wave amplitudes (top panels), instantaneous frequencies (middle panels) of forward-propagating waves ( t = 0 ∼ 2.04 × 10 4 [ −1 e0 ] ) that are after the band-pass filter of 0.04 ∼ 0.06 e0 , and their corresponding S values (bottom panels). Column a corresponds to Case 1 ( ω ph = 0.3� e0 ); column b corresponds to Case 4 ( ω ph = 0.15� e0 ); column c corresponds to Case 5 ( ω ph = 0.075� e0 ). The red dashed lines represent the critical distance in Case 5 simulations, we set q/m to be constant. To keep the unchanged level of thermal fluctuations, i.e., the charge value q needs also to be constant. We assume the electromagnetic thermal fluctuation is determined by the hot electrons rather than the cold electrons by setting up the very small thermal velocity and the large number of super-particles for cold electrons. Therefore, in Eq. (1), since the hot plasma frequency ω ph in Case 4 and Case 5 are 1/2 and 1/4 of the value in Case 1, we need to change the number density of hot particles n to 1/4 and 1/16 of the value in Case 1 for Case 4 and Case 5, respectively.
From (4) and (5), ω ph is related to both optimum and threshold amplitude. In Fig. 8, we give theoretical results of the optimum and threshold amplitudes at different hot electron density levels. As the level goes down, the optimum amplitude gets lower, while the threshold value becomes higher, leading to the disappearance of overlap between the optimum and threshold amplitudes. The theoretical result denotes that the first stage of the nonlinear process becomes weaker, and wave packets can hardly enter the convective growth stage in a lower density case. In Fig. 9, we find the critical distances in Case 1 and Case 4 are greater than the system length, meaning that the first stage of the nonlinear process can exist in the whole simulation region. As for Case 5, its critical distance is about one-quarter of the system length, thus the first stage of the nonlinear process can happen within half of the simulation region around the equator. Figure 10 shows the dynamic spectra of forward-propagating waves at different locations for all three hot electron density levels. Although their backgrounds are at the same thermal fluctuation level, the wave amplitudes are very different. We cannot find hiss emissions in either Case 4 or Case 5, which shows a good agreement with the theoretical result. We find similar wave structures at the initial stages of all cases and locations, as indicated in Fig. 10. Based on the analyses in the previous section, although at small hot electron density levels, these structures are generated through the frequency variation, i.e., the first stage of the nonlinear wave growth process. To check this point, we implement a band-pass filter of e0 ] and calculate their instantaneous frequencies and inhomogeneity factor S, as shown in Fig. 11. From the nonlinear wave growth theory, a suitable value of inhomogeneity factor S (i.e., |S| < 2 ) represents the existence of the nonlinear process. In the bottom panels of Fig. 11, we find some values of |S| are smaller than value 2 along with the forward-propagating wave packets in Case 4 and Case 5, although they are much more sparse compared to those in Case 1. With the reference of red dashed lines which represent the critical distance in Case 5, S values that correspond to the nonlinear process occur more inside the distance, declaring a good agreement with the theoretical result.
In Fig. 12, we calculate J E that is the component of resonant current parallel to the wave electric field E w for all density cases. We see that the value of J E in the low hot electron density case is much smaller. Since J E is closely related to the nonlinear wave growth (Omura 2021), in Case 4 and Case 5, there is not much energy transferred from particles to waves through nonlinear resonant motions. Therefore, although the nonlinear process does exist in the low hot electron density cases, it is strictly suppressed.
To further check the correlation between simulation and theoretical results, we compare wave amplitudes, frequencies, and calculated S values at the equator for all three hot electron density cases, as shown in Fig. 13. Each peak region in the wave amplitude profile represents a forward-propagating packet that is passing through the location. We find that almost all wave packets in Case 1 have corresponding S values that are smaller than or close to unity. As for Case 4 and Case 5, we can hardly find small values of S, meaning most wave packets are not undergoing the nonlinear wave growth process.

Summary and discussion
We conducted a series of electromagnetic particle simulations to check the dependence of generation and nonlinear wave growth process of hiss emissions on simulation parameters. We find the following: Fig. 13 a Wave amplitudes of forward-propagating waves ( t = 0 ∼ 2.04 × 10 4 [ −1 e0 ] ); b instantaneous frequencies, and c corresponding absolute S values in different cases of the hot electron density. Data with wave amplitudes less than 2 × 10 −5 are eliminated. The green dashed line in c represents unity 1. As the gradient of the background magnetic field increases, the wave amplitude generated becomes smaller because of the suppression of the nonlinear process. For a realistic gradient environment (i.e., Case 1), simulation results demonstrate that the generation and growth of plasmaspheric hiss due to the nonlinear process is possible. 2. For the large gradient cases, the convective growth stage plays the key role in the wave growth. However, due to large values of the gradient at high latitudes, nonlinear resonant motions cannot be retained, thus waves stop growing. 3. As the hot electron density decreases, wave growth is suppressed drastically. The intensity of the energyconversion process between particles and waves (denoted by J E ) is varied, fewer wave packets can enter the nonlinear process in the small hot electron density cases.
The present simulation results confirm the correlation between the generation of plasmaspheric hiss and the nonlinear wave growth theory, indicating that hiss emissions are originated from local thermal fluctuations near the equatorial region of plasmasphere and undergo a nonlinear wave growth process through the formation of electron holes or hills. However, because of the limited size and duration of hiss elements, it is difficult to analyze the holes or hills. Although we find good agreements between the theoretical and simulation results, the threshold amplitude we analyzed is derived for rising tone structures. There is no theoretical expression of threshold amplitude for falling tone structures yet. Wu et al. (2020) recently reported that the background magnetic field inhomogeneity is not required for chirping of chorus, but a negative inhomogeneity may cause falling tone structures. The case for a homogeneous background magnetic field should be investigated more. The precipitation or nonlinear scattering of hot electrons by coherent plasmaspheric hiss is also an interesting research topic. We leave them as future studies.