Excitation of Internally Driven ULF Waves by the Drift‐Bounce Resonance With Ring Current Ions Based on the Drift‐Kinetic Simulation

We investigated the excitation mechanism and global distribution of internally driven Pc3‐Pc5 ultralow‐frequency waves in the inner magnetosphere using a kinetic and self‐consistent simulation (Geospace Environment Modeling System for Integrated Studies‐Ring Current model). These waves can be excited through the drift‐bounce resonance with ring current ions injected into the inner magnetosphere during substorms. Initial phase space density of ring current ions is the butterfly‐like distribution with the asymmetry in pitch angle direction. We find two kinds of internally driven ultralow‐frequency waves. First, the dynamic power spectra of the field fluctuations show the excitation of poloidal, toroidal, and compressional Pc5 waves. We find that the excited waves are fundamental modes propagating westward with the azimuthal wave number of m ~ − 50. These waves are excited by the drift resonance with ions having the energy of 50–150 keV. Second, we can reproduce the excitation of poloidal Pc3 waves in the dusk region with m ~ − 30, which are driven by the bounce resonance with ions having energies of 20–50 keV for the first time. We investigate the energy transfer between ions and waves and find that poloidal Pc5 and Pc3 waves are excited mainly by positive energy gradient of the phase space density, while we find the inward gradient of the phase space density at L > 6. The velocity distribution of the growth rate suggests that Pc5 waves are driven by ions with pitch angle of around 90°, while Pc3 waves are excited by ions at oblique pitch angle.


Introduction
Pc5-Pc3 ultralow-frequency (ULF) waves (Pc5, 1.7-6.7 mHz; Pc4, 6.7-22 mHz; and Pc3, 22-100 mHz) are electromagnetic pulsations widely observed in the inner magnetosphere. Pc4-Pc5 ULF waves are generated by external sources such as solar wind dynamic pressure variations and the Kelvin-Helmholtz instability with a characteristic of low azimuthal wave number (m number) of m~1 (e.g., Claudepierre et al., 2008). There is an external source for Pc3 waves in upstream waves in the solar wind (Russell & Hoppe, 1983). Theoretical studies have shown that Pc5-Pc3 can be also excited by internal sources through a wave-particle interaction with ring current ions injected from the magnetotail during substorms (Southwood, 1976). Since Pc5 ULF waves can drive radial transport of radiation belt electrons (e.g., Elkington, 2006;Shprits et al., 2008), to investigate the excitation and global structure of ULF waves is important for understanding the variation of charged particles in the inner magnetosphere.
There are three oscillation modes for ULF waves: poloidal mode (the radial magnetic field fluctuation of dB r and the azimuthal electric field fluctuation of dE ϕ ), toroidal mode (dB ϕ and dE r oscillations), and compressional mode (the parallel magnetic field fluctuation of dB ‖ ). Poloidal mode ULF waves are able to make wave-particle interaction with ring current ions through the azimuthal electric field. This is called the drift-bounce resonance (Southwood, 1976). The resonance condition is where ω is the wave angular frequency, m is the m number of poloidal mode waves, n is an integer (corresponding to the number of wavelengths across which particles drift during one bounce motion in the wave rest frame), ω d is the drift angular frequency, and ω b is the bounce angular frequency. Note that ion drift and bounce frequencies are both energy and pitch angle dependent (Hamlin et al., 1961).
6W qB eq L 2 R 2 E 0:35 þ 0:15 sin α eq À Á ; (2) where W is the ion energy, L is the L value, R E is the Earth radius, B eq is the equatorial magnetic field magnitude, m p is the ion mass, q is the electric charge, and α eq is the equatorial pitch angle of ions. There are two solutions for the drift-bounce resonance: the drift-bounce resonance (mω d~ωb ) and the bounce resonance (ω~ω b ). Standing wave structure associated with the drift-bounce resonance is an odd mode (symmetric field line oscillations between the magnetic equator) when the integer n is even, while it is an even mode (asymmetric field line oscillations between the magnetic equator) when n is odd (Southwood & Kivelson, 1982).
There have been a number of observational studies discussing the drift-bounce resonance excitation of internally driven ULF waves. While there have been several studies which observed fundamental poloidal waves associated with the drift resonance (e.g., Dai et al., 2013;Yamamoto et al., 2018), previous observations overwhelmed the drift-bounce resonance excitation of second harmonic poloidal waves (e.g., Le et al., 2017;Min et al., 2017;Oimatsu et al., 2018;Takahashi et al., 1985Takahashi et al., , 2018Yamamoto et al., 2019). These observations suggest that internally driven ULF waves excited by the drift-bounce resonance tend to have a poloidal polarization and propagate westward with high m numbers (|m|~70-200) (e.g., Dai et al., 2013;Takahashi et al., 2018). On the other hand, toroidal mode can be excited by the mode conversion from poloidal mode waves in inhomogeneous plasmas (Mann & Wright, 1995).
Previous simulation studies of internally driven ULF waves have been based on magnetohydrodynamic (MHD) equations of cold plasma and investigated the field line resonance excitation of ULF waves (e.g., Degeling et al., 2018). However, these MHD models do not include the energy-dependent drifts of ions and self-consistent development of the fields, which are essential for internally driven ULF waves. Recently, Yamakawa et al. (2019) detected internally driven Pc5 waves by using a global drift-kinetic model, which assumed an isotropic phase space density (PSD) ion injection at the night side with north-south symmetry. These Pc5 waves are fundamental mode waves propagating westward excited by the drift resonance with ions having the energy of 80-120 keV. However, the drift-bounce resonance excitation of ULF waves was not detected in the symmetric PSD case.
In this study, we investigate the excitation of internally driven ULF waves associated with the drift-bounce resonance based on a global drift-kinetic simulation for the first time. We use Geospace Environment Modeling System for Integrated Studies-Ring Current (GEMSIS-RC) model (Amano et al., 2011) as well as Yamakawa et al. (2019). Section 2 describes the simulation model and setups. In section 3, we present simulation results of ULF waves and investigate the growth rate of ULF waves. Section 4 discusses simulation results by comparing with Yamakawa et al. (2019), and section 5 concludes this study.

Model Description
GEMSIS-RC model is a kinetic simulation code to solve five-dimensional drift-kinetic equation for the ion PSD f(x, v || , μ) and Maxwell's equations in a self-consistent manner under the assumption that the magnetic moment of μ is conserved (Amano et al., 2011). Note that x and v || represent position in the 3-D space and the parallel velocity, respectively. Since this model is capable of assuming the force-imbalanced state and including the ion drift and bounce motions unlike previous MHD models, it enables us to simulate internally driven ULF waves driven by ring current.
The coordinate system in GEMSIS-RC model follows the modified dipole coordinate (Kageyama et al., 2006). Note that the coordinate x 1 is the direction along the dipole field line, x 2 is the direction in which the L value changes, and x 3 is the azimuthal direction. The simulation domain and the grid number are the same as those used in Yamakawa et al. (2019) except for the number of x 3 grids, that is, we use the parameters of the minimum L value of L min = 3.6 R E , the maximum L value of L max = 6.6 R E , the minimum radial distance of R min = 3.0 R E , the minimum angle from the pole of θ min = π/3, and the stretching parameter a ¼ 100 R 2 E .
The numbers of grids are N 1 = 64, N 2 = 32, and N 3 = 1,024 (see Yamakawa et al., 2019, for detail). We enhance the number of x 3 grids, since we want to reproduce internally driven ULF waves with higher m number than that of Yamakawa et al. (2019) (|m|~20). This is because ULF waves in spacecraft observations usually have the m number of |m| > 50 (e.g., Dai et al., 2013;Takahashi et al., 2018). For the velocity space, we use the same grid structure as Yamakawa et al. (2019): −4,500 ≤ v || km/s ≤ 4,500 and 10 −4 ≤ μ keV/nT ≤ 10 0 with the numbers of grids, N v‖ = 32 and N μ = 32, respectively.
The initial condition of electromagnetic fields is only the dipole magnetic field with a dipole moment of M E = − 8.05 × 10 15 Tm 3 , and the electric field is set to be zero everywhere at the initial state. In this simulation, we include the kinetics of ions using five-dimensional PSD of f(x, v || , μ), while we treat electrons as fluid with no mass under the assumption that electrons move very quickly and satisfy charge neutrality condition. The initial pressure distribution is set to be localized at the nightside to simulate ion injection from the magnetotail as well as Yamakawa et al. (2019); the equatorial density is where L 0 = 5.0 R E , ΔL = 0.5 R E , ϕ 0 = π, and Δϕ = π/16. We assume the same density distribution at the initial state as Yamakawa et al. (2019). Although the density and temperature seem higher compared to statistical studies of the injection (e.g., Denton et al., 2005;Lavraud et al., 2005), we use the above initial condition to save the calculation time. Note that we include only protons for the ion PSD. In contrast to symmetric isotropic (T ‖ = T ⊥ = 16 keV ≡ T 0 ) PSD with loss cones case in Yamakawa et al. (2019), the initial PSD is the butterfly distribution with asymmetry in pitch angle direction. Butterfly distribution has been reported by previous in situ spacecraft observations (e.g., Shi et al., 2016). We add the asymmetry of the PSD since the drift-bounce resonance excitation of ULF waves is not found if we put symmetric butterfly distribution. We also assume that the asymmetry of the PSD is possible when injected ions are affected by the tilt of the dipole field. The initial equatorial PSD is given by

Journal of Geophysical Research: Space Physics
Note that k is the Boltzmann's constant. Figures 1a and 1b show the velocity distribution and pitch angle dependence at each energy of the initial equatorial PSD at L = 5, MLT = 0. If we integrate the PSD over the velocity space, we can get the number density of N 0 . The initial PSD inside loss cones is zero (see Figure 1a), and we adiabatically map the equatorial PSD along the field line to the off-equatorial region (Amano et al., 2011;Yamakawa et al., 2019). Figure 1c shows the equatorial P ⊥ distribution at t = 0 s.
We use the same boundary condition as Yamakawa et al. (2019), no field fluctuation at all boundaries, which means that we have finite layers in x 1 and x 2 boundaries in order to prevent the wave reflection at the boundaries, where field fluctuations are damped by an artificial smoothing. We assume that particles undergo mirror reflection at ionospheric (x 1 ) boundaries and are not injected at radial (x 2 ) boundaries. In addition to this hot population of ions, there is a background cold plasma, which does not move and has no MLT dependence with the density given by where N bg = 10 5 cm −3 , Z bg = 6.37 × 10 6 m, H bg = 3.16 × 10 3 m, and β bg = 0.1. The density is based on the empirical model (Morioka et al., 2008;Sato, 1998) By using this background density, the equatorial Alfven velocity becomes about 500-700 km/s. We use the constant time step of t = 0.05 s for all variables.

Excitation and Global Distribution of ULF Waves
In this section, we present simulation results of the excited ULF waves. First, we show the transport of energetic ions. Next, excitation and global distribution of ULF waves are presented. Figures 1c and 1d show the equatorial P ⊥ distributions at t = 0 s ( Figure 1c) and at t = 1,200 s ( Figure 1d). Since the initial condition of ions is in the force-imbalanced state, ions in the initial high-pressure region move toward duskside and dawnside. Then, ions undergo westward magnetic drifts, which lead to the pressure distribution with dawn-dusk asymmetry as shown in Figure 1d as well as the symmetric PSD case (Yamakawa et al., 2019).  Figure 2g illustrates the time evolution of dB r and dE ϕ components using high-pass filter (>10 mHz) in order to detect high-frequency oscillations. We detect internally driven ULF waves with two different frequencies.
First, the dynamic power spectra in Figures 2h-2l show the clear excitation of poloidal, toroidal, and compressional Pc5 waves with the frequency of 2-6 mHz. The wave frequency gradually decreases with time as well as the symmetric PSD case (Yamakawa et al., 2019). We estimate the m numbers to be m = −50 for dB r , m = −48 for dB ϕ , and m = −52 for dB ‖ with westward propagation. Note that we calculate the m numbers at one local point by analyzing the phase differences between the two nearest grid points in the azimuthal direction. We have confirmed that the m numbers do not change so much (Δm~1) when we calculate the m number by analyzing the phase differences at other grid points (e.g., two second nearest grid points). We find that the field line mode structure is a symmetric odd mode (see Figure S1 in the supporting information). The frequency of poloidal Pc5 waves indicates that they are fundamental mode waves driven by the drift resonance. Compressional Pc5 waves are also excited as shown in Figure 2e. Figures 2e and 2f indicate that the plasma pressure is anticorrelated with the magnetic pressure. We will discuss the excited Pc5 waves in more detail in Figure 3 and next section.

Journal of Geophysical Research: Space Physics
Second, we detect the excitation of Pc3 variations (22-40 mHz) unlike the symmetric PSD case (Yamakawa et al., 2019) as shown in Figures 2h-2l. We refer to these oscillations as poloidal Pc3 waves, since dB r has larger amplitudes than dB ϕ . We estimate the m number of poloidal Pc3 waves using the same method, giving the value of m = −25. From the resonance condition of Equation 1, we estimate that Pc3 waves are excited by the bounce resonance (n = 1, ω~ω b ) because of the wave frequency and the m number. We have confirmed that the field line mode structure is an asymmetric even mode (see Figure S2). In addition, dE ϕ leads dB r by 90°around t = 1,200 s as illustrated in Figure 2g, and this is the characteristic of the second harmonic poloidal mode oscillations in northern hemisphere.
We find that Pc5-Pc3 waves are excited near the equator. Figure 3 shows the spatial distribution of Pc5 waves at t = 1,200 s at MLAT = 1°. We detect the excitation of Pc5 waves around 4.8 < L < 6 and 12 < MLT < 18 at t = 1,200 s for all mode waves as shown in Figure 3. If we integrate these power spectra over the Pc5 range, we have confirmed that the wave power of compressional Pc5 waves is the strongest, and poloidal and (a)  Figure 2 shows that radial electric field perturbation of dE r is in phase with dB ‖ around t = 900 s, which indicates that toroidal perturbation of dE r is caused by the compressional perturbation of dB ‖ . These results suggest that the toroidal mode excitation is secondary and triggered by the excitation of compressional Pc5 waves. Compressional Pc5 waves have the stronger wave power than poloidal and toroidal Pc5 waves. We discuss compressional Pc5 waves in section 4. Figure 3a shows that poloidal Pc5 waves are excited at 4.5 < L < 6 and L > 6. There is a region where the poloidal wave power is weak around MLT = 16. We consider that this is because ions give the energy to toroidal and compressional Pc5 waves, since the power spectra of these two mode waves are the strongest around MLT = 16 as shown in Figures 3b and 3c. The wave frequency of the excited Pc5 waves in Figures 3d-3f is in the range of 2-6 mHz. The frequency is higher in the duskside. This is due to the higher m number in the duskside as shown in Figures 4h and 4i, since the wave frequency of the drift resonance is proportional to the m number (Southwood et al., 1969). Note that the excited Pc5 waves propagate westward and can be seen in the dawnside and nightside after t = 1,200 s.  condition (Amano et al., 2011), it is excluded from the discussion of the drift-bounce resonance in the next section. We discuss poloidal Pc3 waves in the duskside (17 < MLT < 20) and Pc5 waves in the dayside dusk region in the next section.

Growth Rate of Poloidal Waves
We have identified the excitation of Pc5 waves through the drift resonance (n = 0) and Pc3 waves through the bounce resonance (n = 1) in the previous section. In this section, we investigate the growth rate and energy transfer of poloidal Pc5 and Pc3 waves. The growth rate of poloidal waves was derived by Southwood et al. (1969).
where f is the ion PSD, W is the ion energy, and v res is the velocity of resonant ions. If the growth rate is positive (negative), ions lose (gain) energy, and poloidal waves grow (decay). Here, the growth rate can be separated into these two terms.

Journal of Geophysical Research: Space Physics
Since dW dL is negative for the westward propagating ULF waves, positive energy gradient of the PSD (positive γ 1 ) or the inward gradient of the PSD (positive γ 2 ) is suitable for the waves to grow from Equations 9-11 (Yamakawa et al., 2019). Then, we discuss each term and total contribution of the growth rate for poloidal Pc5 and Pc3 waves. Note that the growth rate of the Pc3 waves tend to be higher than that of the Pc5 waves, since dW dL is in proportion to the wave frequency (Southwood et al., 1969).
First of all, we investigate ion energy and pitch angle having positive growth rate of Equation 8. Figure 5 shows the velocity distribution of each term and total growth rate as functions of pitch angle and the magnetic moment at t = 1,200 s at L = 5.5, MLT = 18, and MLAT = 1°(the same point as Figure 2). From Equations 1-3 and the frequency of poloidal waves, we can draw the theoretical resonance curves in  . Note that we choose n = 0 for poloidal Pc5 waves and n = 1 for poloidal Pc3 waves. Also, we use the wave frequency of ω in Figure 3d and the m number in Figure 3g for the Pc5 waves and the wave frequency of ω in Figure 4d and the m number in Figure 4g for the Pc3 waves. Although whether ions can make the drift and drift-bounce resonant interaction with ULF waves depend on the sign of dE ϕ along the trajectory of ions, we would rather focus on ion energy and pitch angle giving the energy to waves at MLAT = 1°. There are several points in the velocity space where the growth rate is positive along the resonance curve for both the Pc5 and Pc3 waves as shown in Figures 5c and 5f. This is because of positive energy gradient of the PSD (positive γ 1 ), while the outward gradient of the PSD (negative γ 2 ) suppresses the growth of the waves as illustrated in Figures 5a and 5b and 5d and 5e at this local point. We find the energy and pitch angle contributing to the wave growth at the point, the energy of W = 40-50 keV and pitch angle of α = 80-100°for the Pc5 waves (drift resonance), while W = 30-40 keV and α = 35-65°, 115-150°for the Pc3 waves (bounce resonance). Figure 6 shows the histogram of ion pitch angle and energy, which have the most contribution to the wave growth (the point in the velocity space where the growth rate is maximum along the theoretical resonance curve). Note that we select the spatial point where the poloidal wave power is greater than 5 × 10 −3 nT 2 /Hz for the Pc5 waves and 5 × 10 −6 nT 2 /Hz for the Pc3 waves. Red (blue) color means the number of spatial points where γ 1 (γ 2 ) term is dominant. As indicated in Figures 6a and 6b, Pc5 waves are driven by ions around pitch angle of α~90°(60-120°), while Pc3 waves are driven by ions at oblique pitch angle (α~50°, 130°). This result represents the difference between the drift and bounce resonances, which is consistent with the theory (Southwood & Kivelson, 1982). Figure 6c indicates that the free energy for the Pc5 waves is almost in the range of 50-150 keV, which is comparable to the result of the symmetric PSD case (Yamakawa et al., 2019). As for the Pc3 waves, ions of 20-50 keV can contribute to the excitation of the waves as illustrated in Figure 6d. We find that both the Pc5 and Pc3 waves are driven mainly by positive energy gradient of the PSD rather than the inward gradient of the PSD. Number of spatial points (d) Figure 6. The histogram of pitch angle distribution contributing to the growth of (a) Pc5 and (b) Pc3 poloidal mode at MLAT = 1°at t = 1,200 s. Also, the energy of ions contributing to the wave growth of (c) Pc5 and (d) Pc3 are shown in bottom two panels. The color of each bin represents the number of spatial points where γ 1 term (red) and γ 2 term (blue) are dominant. Note that we select the spatial point where the poloidal wave power is higher than 5 × 10 −3 nT 2 /Hz for the Pc5 waves and 5 × 10 −6 nT 2 /Hz Pc3 waves.
We calculate each term of the growth rate in Equations 9 and (10) at each spatial point in order to discuss the wave growth as well as the symmetric PSD case (Yamakawa et al., 2019). Figure 7 shows the global distribution of each term and total growth rate for the drift resonance (Figures 7a, 7c, and 7e) and the drift-bounce resonance (Figures 7b, 7d, and 7f) at t = 1,200 s at MLAT = 1°. We calculate γ 1 and γ 2 using the energy and the magnetic moment of the point in the velocity space where the growth rate is maximum along the theoretical curves satisfying the resonance condition, which are shown in Figure 5.
We find the inward gradient of the PSD (positive γ 2 ) in Figure 7c at L > 6 or in the front part of the drifting protons (12 < MLT < 17). This inward gradient can cause the excitation of poloidal Pc5 waves in Figure 3a in this region. We also find the inward gradient of PSD (positive γ 2 ) in the region of 17 < MLT < 24 and L > 6 as shown in Figure 7c. However, we do not see the excitation of poloidal Pc5 waves in Figure 3a but poloidal pulsations with the frequency about 1 mHz in this region. We consider that the injected ions do not have sufficient energy to excite Pc5 waves in this region. On the other hand, Figures 7a and 7c indicate that positive energy gradient of the PSD (positive γ 1 ) mainly contributes to the growth of poloidal Pc5 waves at L < 6. These characteristics of the growth rate for poloidal Pc5 waves are also seen in the symmetric PSD case (Yamakawa et al., 2019). Figure 7f indicates that the region where poloidal Pc3 waves are excited have the positive growth rate in the dusk side. While at L > 6 where Pc3 waves are excited, we detect the inward

Discussion
The simulation results in section 3 show the excitation of two kinds of internally driven ULF waves: the Pc5 waves (the drift resonance, n = 0) and the Pc3 waves (the bounce resonance, n = 1). In this section, we compare the simulation results with symmetric PSD injection case (Yamakawa et al., 2019).
As shown in Figure 3, we detect the excitation of poloidal, toroidal, and compressional Pc5 waves. The power spectra of Pc5 waves are about 1 order of magnitude higher than the symmetric PSD case (Yamakawa et al., 2019). If we simulate under the initial condition of the symmetric PSD as well as Yamakawa et al. (2019) with the azimuthal grid numbers of N 3 = 1,024, the power spectra of Pc5 waves are comparable to this study, while we cannot detect the excitation of poloidal Pc3 waves. These results indicate that the increase of the power spectra of the Pc5 waves is because of the enhancement of the spatial resolution along the azimuthal direction (from N 3 = 128 to N 3 = 1,024). This enhancement also causes the increase of the m number from |m|~20 to |m|~50. The increase of the m number is less than that of the grid numbers along the azimuthal direction. This result indicates that the m number in our simulation is well converged. In contrast to Pc5 waves, the excitation of Pc3 waves was not detected in symmetric PSD case (Yamakawa et al., 2019). Therefore, the inclusion of the asymmetry of the PSD causes the drift-bounce resonance excitation of poloidal Pc3 waves in the simulation.
Figures 3a-3c show that compressional Pc5 waves have about 3 orders of magnitude higher power spectra than poloidal and toroidal Pc5 waves. We also confirm that the Pc5 waves have the fundamental mode structure as shown in Figure S1. This type of pulsations is also seen in the symmetric PSD case (not shown in Yamakawa et al., 2019). Here we discuss the polarization of the excited Pc5 waves. Southwood et al. (1969) has proposed that poloidal mode waves are associated with the drift-bounce resonance in a homogeneous field. Recently, Yamamoto et al. (2018) observed the excitation of compressional pulsations propagating westward with the m number of m = − 50 around the magnetic equator. The waves were odd mode waves and in the drift resonance with energetic ions. It is interesting to see that these wave properties are actually seen in the simulation as shown in Figure 3. Yamamoto et al. (2018) consider that the Pc5 waves are mainly compressional because dB r of fundamental mode is small, and compressional mode becomes prominent around the magnetic equator (Takahashi et al., 1992). Since we show the global distribution of Pc5 waves at MLAT = 1°in Figure 3, this interpretation can qualitatively explain the excited Pc5 polarization. It also should be noted that we do not include the field line resonance, since field fluctuations are damped by an artificial smoothing, which prevents the wave reflection at ionospheric boundaries in this simulation (Amano et al., 2011). The lack of the field line resonance can suppress the power of shear Alfven waves. As shown in Figure 4a, we detect the second harmonic poloidal mode in the dusk and night regions. However, the amplitude of poloidal Pc3 waves in Figure 2 is quite small compared to the observations of the second harmonic mode. We have confirmed that if we enhance the initial ion density, the power spectra of the Pc5 and Pc3 waves increase, while the basic properties such as the wave frequency and the m number do not depend on the initial density distribution so much. In order to discuss the excitation of ULF waves more quantitatively by comparing with observations, we will need to include the feedback from the ionosphere and consider ion density and PSD, which are similar to those observed by spacecraft observations by coupling with empirical models. Figures 6a and 6c show that ions of 50-150 keV and α~90°contribute to the wave growth of Pc5 waves. This is consistent with the theory that ions with pitch angle around 90°can most effectively make drift resonant wave-particle interaction with a symmetric wave mode (Southwood & Kivelson, 1982). We confirm that the energy is higher in the dayside, and this is due to energy-dependent drifts of energetic protons. Figures 6b  and 6d indicate that Pc3 waves are in the bounce resonance with protons of 20-50 keV at oblique pitch angle. We consider that Pc5 waves are distributed in MLT 12-18 in Figure 3a whereas Pc3 waves are excited in MLT 17-20 in Figure 4a since ion energy responsible for the excitation of Pc5 waves is higher than that for Pc3 waves. Figure 2 shows that the fluctuations of electromagnetic fields are dissipating from t = 1,200 s. We suggest that one reason is the numerical dissipation of GEMSIS-RC model and the other is due to the growth rate of ULF waves. We have confirmed that the growth rate of Pc5 waves in Figure 7e turns from positive 10.1029/2020JA028231

Journal of Geophysical Research: Space Physics
to negative at t = 1,200 s at the spatial point selected in Figure 2 (L = 5.5 and MLT = 18). This is due to the negative γ 2 term (outward gradient of the PSD) as shown in Figure 7c. We consider that the outward gradient of the PSD around MLT 18-20 is caused by the transport of low-energy ions (20-50 keV), which suppresses the wave excitation in Figure 2. Figure 7 shows that the inward gradient of the PSD at L > 6 and positive energy gradient of the PSD at L > 4.5 provide the free energy to excite poloidal waves. We suggest that the energy-and L-shell-dependent drifts of injected protons affect the value of the spatial and energy gradient of the PSD. If we change the type of the initial injection, the spatial distribution of ULF waves can be altered.
The simulation results indicate that GEMSIS-RC model is capable of reproducing both drift and drift-bounce resonance excitations of ULF waves by describing energy-dependent drift and bounce motions of energetic ions. However, this model lacks of some important processes such as the ionospheric current, the convection electric field, and the dynamics of the plasmasphere (Amano et al., 2011;Yamakawa et al., 2019). Therefore, these physical processes will be included for the quantitative discussion of internally driven ULF waves.

Conclusion
We investigate the excitation of internally driven ULF waves by using a drift-kinetic (GEMSIS-RC) simulation. We assume that ions injected from the magnetotail have the butterfly-like velocity distribution with asymmetry in pitch angle direction. We find the internal excitation of not only Pc5 waves but also Pc3 waves. Pc5 waves are driven by the drift resonance, and their frequencies gradually decrease from 6 to 2 mHz as well as symmetric PSD case (Yamakawa et al., 2019). We can reproduce the excitation of poloidal Pc3 waves associated with the bounce resonance in the dusk region (22-40 mHz) for the first time. By analyzing the velocity distribution of the growth rate, poloidal Pc5 waves are excited by ions around pitch angle of 90°, while poloidal Pc3 waves are driven by ions at oblique pitch angle. We investigate the energy transfer between ions and poloidal waves. While we find the inward PSD gradient at L > 6, poloidal Pc5 and Pc3 waves are excited mainly due to positive energy gradient of the PSD. We suggest that this study will be the basis of the future self-consistent modeling of internally driven ULF waves.

Data Availability Statement
The data will be available at the UTokyo Repository (http://doi.org/10.15083/00079433).