Reply to comment by R. M. Thorne and R. B. Horne on Khazanov et al. [2002] and Khazanov et al. [2006]

Abstract

[1] It is well-known that the effects of electromagnetic ion cyclotron (EMIC) waves on ring current (RC) ion and radiation belt (RB) electron dynamics strongly depend on such particle/wave characteristics as the phase-space distribution function, frequency, wave-normal angle, wave energy, and the form of wave spectral energy density. The consequence is that accurate modeling of EMIC wave spatial-temporal-spectral distributions and RC particles requires robust inclusion of the interdependent dynamics of wave growth/damping and wave propagation/refraction/reflection, along with wave tunneling and mode conversion and particles. Such a self-consistent RC-EMIC wave model is being progressively developed by Khazanov et al. [2002, 2003a, 2006, 2007]. This model is based on a system of coupled kinetic equations for the RC and EMIC wave power spectral density with explicit inclusion of the ray tracing equations. [2] The theoretical formalism for RC ions and RB electrons is well established and is based on gyroaveraged and bounce-averaged kinetic equations that have been developed over the years and systemized in the book by Khazanov [1979]. The application of this formalism to RC ions was continued by Khazanov and Kozyra at the University of Michigan during 1991–1994 and formed a large part of the Ph.D. dissertation work of Fok [1993] and of Jordanova [1995]. Khazanov et al. also applied this formalism to the study of ionosphere-plasmasphere transport of suprathermal electrons [Khazanov et al., 1992, 1994], and global photo [Khazanov et al., 1996; Khazanov and Liemohn, 2002] and plasma sheet [Khazanov et al., 1998] electron transport. This formalism was also generalized to study relativistic electron transport [Khazanov et al., 1999, 2000], as well as different aspects of RC and RB electron formation using various magnetospheric electric and magnetic field topologies [Khazanov et al., 2003b, 2004b, 2004c]. All these above-mentioned studies are the heritage of our RC-EMIC wave model that was presented by Khazanov et al. [2006]. [3] Thorne and Horne [2007, hereinafter referred to as TH2007] call the Khazanov et al. [2002, 2006] results into question in their Comment. The points in contention can be summarized as follows. TH2007 claim that (1) important damping of waves by thermal heavy ions is treated incorrectly in our model, and Landau damping during resonant interaction with thermal electrons is not included; (2) EMIC wave damping due to RC O+ is not included in our simulation of the 2–7 May 1998 storm; (3) nonlinear processes limiting EMIC wave amplitude are not included in our model; (4) growth of the background electromagnetic fluctuations to a physically significant amplitude must “occur during a single transit of the unstable region” with subsequent damping in the vicinity of the bi-ion latitude, and consequently the bounce-averaged wave kinetic equation employed in the code is not valid. Our reply will address each of these points as well as other criticisms mentioned in the Comment. [4] TH2007 is focused on two of our papers that are separated by 4 years. Significant progress in the self-consistent treatment of the RC-EMIC wave system has been achieved during those years. The paper by Khazanov et al. [2006] presents the latest version of our model, and in this Reply we refer mostly to this paper. [5] EMIC wave damping due to thermal heavy ions and electrons has always been included in our studies (see Khazanov et al. [2002], section 2; Khazanov et al. [2003a], section A3.2; Khazanov et al. [2006], sections 3.1, 4.1, and 5.2; Khazanov et al. [2007]). Particularly, an essential part of our newest study by Khazanov et al. [2007] uses the same formalism and is devoted to energy deposition to thermal plasmaspheric electrons due to Landau damping of EMIC waves. Khazanov et al. [2006, equation (22)] explicitly include resonant absorption by thermal ions and electrons in the RC-EMIC wave model. In other words, the damping rate in the right-hand side of the equation is a result of integrating the local resonant damping rate along the entire ray phase trajectory (r, θ), including that in the vicinity of the bi-ion frequencies, over the wave bounce period. So, the TH2007 statement that “the important damping of waves by thermal heavy ions near the bi-ion location is not explicitly treated by the bounce-averaged kinetic equation (1)” is incorrect. The thermal plasma in our model is currently treated independently from the self-consistent dynamics of the RC-EMIC wave system. It means, particularly, that we do not take into account a change of the thermal electron and ion temperatures due to resonant Landau and cyclotron wave damping. The core plasma is assumed to be Maxwellian with a temperature of 1 eV for both electrons and ions. For that temperature, only electron Landau damping is important, and wave damping by thermal ions is negligible. [6] The TH2007 statement that “the reflection coefficient was essentially assumed to be unity in the [Khazanov et al., 2006] formulation” also appears to result, in part, from a misunderstanding because the tunneling of EMIC waves through the region of the bi-ion reflection is another loss mechanism also included in our RC-EMIC wave model. For the plasma parameters adopted in the one the particular study by Khazanov et al. [2006], only a negligible portion of the EMIC wave energy was found able to tunnel across the reflection region thereby allowing us to exclude the effect of tunneling in the wave kinetic equation for this case and safely assume that reflection from the bi-ion region is essentially perfect. The later has nothing to do with EMIC wave resonant damping by thermal heavy ions/electrons, which in our modeling is able to operate as a wave energy absorber. Therefore the criticism listed in statement (1) in section 1 of this Reply is simply incorrect. [7] During the main phase of major storms RC O+ may dominate [e.g., Hamilton et al., 1988; Daglis, 1997] and, as a result, contribute to strong damping of the He+-mode EMIC waves [Thorne and Horne, 1997]. Although there is no doubt that this damping process is critically important, we have serious concern over the ability of the RC model used in the paper of Thorne and Horne [1997] to adequately represent the situation during the main and early recovery phase of a storm. Let us provide observational results that strongly support our statement. Braysy et al. [1998] reported observations of EMIC waves obtained by the Freja satellite and provided remarkable results and conclusions. Particularly, they observed oxygen band waves for about 7 h during the later part of the main phase of the 2–8 April 1993 storm. Since the estimated drift time for RC O+ is only 2–4 h, one would expect to find oxygen band waves at different MLTs. However, all oxygen waves were found in the evening-midnight MLT sector and, in particular, none were observed in the prenoon sector. This implies a very asymmetric O+ RC during the main phase and suggests that the RC oxygen ion loss rate is considerably faster than the drift speed. As emphasized by Braysy et al. [1998], these results are difficult to explain in terms of charge exchange and Coulomb scattering and suggest that the production of EMIC waves contributes significantly to RC O+ decay during the main and early recovery phases. In other words, owing to generation of the oxygen band EMIC waves, most RC O+ precipitates before reaching the dusk MLT sector. [8] These observations clearly demonstrate that to adequately take into account He+-mode energy absorption by RC O+, the O+-mode EMIC waves should concurrently be included in global simulations. While O+-mode EMIC waves are not yet included in our model, this will be completed in the near future. In any case, Table 1 in the work of Thorne and Horne [1997] was generated from a simulation without oxygen band waves, and it is unlikely that the listed RC O+ parameters adequately represent the situation during the studied storm, especially in the dusk sector, for which all the calculations were presented. In addition, Thorne and Horne [1997] used single bi-Maxwellian fits to the simulated RC O+ and H+ distribution functions prior to calculating growth/damping rates. As shown below, this method incorrectly predicts wave growth/damping and the resulting impact on the RC. [9] Next, let us evaluate the validity of excluding He+-mode damping by RC O+ in the 2–7 May 1998 storm simulation that was presented in our paper [Khazanov et al., 2006]. Using the RC kinetic model of Jordanova et al. [1998], Farrugia et al. [2003] found that RC O+ content did not exceed 30% during the main phase of this storm. This estimate was obtained from a global simulation similar to that used by Thorne and Horne [1997], which did not include oxygen band waves. Therefore as follows from the conclusion by Braysy et al. [1998], Farrugia et al. [2003] likely overestimated the RC O+ content during the event. In addition, the calculations of Thorne and Horne [1997] clearly demonstrated that the above RC O+ percentage cannot significantly suppress the He+-mode amplification and only slightly influences the resulting growth; inclusion of 26% O+ in the RC population causes the net wave gain decrease by 20% only. It is for this reason that we chose to initially exclude RC O+ in our particular simulation of 2–7 May 1998. Therefore despite the fact that RC O+ can dominate during the main phase of major storms and causes an additional wave damping, the criticism summarized in statement (2) above should not be a significant factor in simulations of the 2–7 May 1998 storm studied by Khazanov et al. [2002, 2006]. [10] Let us now address the criticism that nonlinear processes, which limit EMIC wave amplitude, are excluded from our model. The nonlinear interaction of large amplitude EMIC waves, for example, the modulational instability that results in generation of solitons and is described by the derivative nonlinear Schrodinger equation [Gamayunov and Khazanov, 1995, and references therein], leads to phase correlation, and in such a system the wave-ion interaction is quite different in comparison with a quasi-linear approach. Another possible mechanism of nonlinear EMIC wave saturation is due to lower hybrid wave generation with subsequent nonlinear scattering and Landau damping on thermal plasma [e.g., Gamayunov et al., 1992; Khazanov et al., 1997, 2004a]. In order to describe the latter nonlinear process, a full kinetic particle-in-cell (PIC) code should be used (for initial results, see Singh and Khazanov [2004] and Singh et al. [2007]), and the hybrid model suggested for use by TH2007 is not an appropriate tool for modeling that process because the hybrid model treats electrons as a massless fluid. [11] At present, our model is based on quasi-linear equations, and the validity of the quasi-linear approach has been carefully monitored from the first version of the model to the most recent (see equation (8) and following text in the work of Khazanov et al. [2003a]). The quasi-linear validity criterion employed in our model is based on the test particle simulation of Kuramitsu and Hada [2000], who showed when quasi-linear diffusion is consistent with nonlinear diffusion. The quasi-linear EMIC wave saturation takes place during most of the storm time, and the introduced criterion restricts wave energy during the main and early recovery phases only when nonlinear EMIC wave saturation likely takes place. During further development of our model, the possibilities for strong nonlinear wave-particle and/or wave-wave interactions will be taken into account, as needed, by using, for example, PIC simulation. Note that the criterion currently used in our model restricts wave amplitudes to values of about 10 nT, which is consistent with the maximum observed EMIC wave amplitudes [e.g., Erlandson and Ukhorskiy, 2001]. Consequently, although nonlinear processes are not yet explicitly included in our model, the nonlinear saturation level imposed in the model is reasonable. [12] The TH2007 statement that growth of background fluctuations to physically significant amplitude must occur during a single transit of the unstable region is based on the calculations of Thorne and Horne [1997]. First, let us note that we are able to reproduce the Thorne and Horne results for path-integrated gain using our code and their modeling parameters (not shown here). In the paper of Thorne and Horne [1997], the RC H+ and O+ distribution functions were obtained from simulation of the November 1993 magnetic storm using the Michigan RC-Atmosphere interaction Model (RAM) [Kozyra et al., 1997] and then fitted by bi-Maxwellian distribution functions (see Table 1 in their paper). It is not clear whether the Michigan RAM particle distributions in their paper were obtained with or without feedback from He+-mode EMIC waves. In other words, it is not clear whether the RAM simulation included wave-ion scattering or did not. Let us examine these two possibilities using results from the Michigan RAM and our model. [13] Assuming that the calculations presented below are rather general and should not depend on a particular storm, we refer to the 2–7 May 1998 magnetic storm (see Khazanov et al. [2006] for more details). Turning off EMIC wave damping by thermal plasma (which is normally included in our RC-EMIC wave model), we calculate the maximum equatorial growth rate for He+-mode EMIC waves using a bi-Maxwellian fit to a simulated RC H+ distribution function. Figure 1 shows the ratio of the local equatorial growth rate calculated using a bi-Maxwellian fit, γbi-max, to a corresponding growth rate obtained from the simulated distribution function without any approximation, γexact. This figure was generated at 80 h after 0000 UT on 1 May 1998 during the early recovery phase. Figure 1 (left) shows results from the Michigan RAM without including wave feedback, Figure 1 (middle) corresponds to the Michigan RAM results with an empirical wave model included as described by TH2007 in the section 2 of the Comment (see also [Kozyra et al., 1997]), and Figure 1 (right) represents the results from our RC-EMIC wave model. The corresponding EMIC wave distributions from our model and the Michigan RAM can be found in the work of Khazanov et al. [2006, Figures 6 and 8]. Without wave feedback, the bi-Maxwellian fit most often overestimates the local equatorial growth rate by at least a factor of two. This overestimation increases dramatically (up to a factor of 30) if RC-EMIC wave scattering is included in the global simulation, indicating different response of γexact and γbi-max on wave feedback. Note that comparing the wave growth and damping rates for a bi-Maxwellian and bi-kappa distributions under magnetospheric conditions, Xue et al. [1996] came to a similar conclusion. The local RC H+ distribution function is not only affected by the local wave distribution but also depends on the prehistory. As follows from Figure 1, the wave feedback depends strongly on the EMIC wave model used, where there is an overestimation by a factor 9 for the Michigan RAM, and by a factor 30 for our model. Moreover, EMIC wave feedback can even cause γexact to be negative while γbi-max > 0 (see gray color in Figure 1). [14] Figure 1 clearly demonstrates that use of a bi-Maxwellian fitted distribution for calculating growth rate routinely overestimates value predictions, especially during and/or after geomagnetically active periods when EMIC wave feedback is crucial for the fine structure of the RC distribution. This conclusion is not only true for hour 80 shown in Figure 1 but for the entire 2–7 May 1998 storm simulation (not shown). The presented theoretical result is strongly supported by observations [Anderson et al., 1996a]. These authors analyzed the proton cyclotron instability in the Earth's outer magnetosphere, L > 7, using Active Magnetosphere Particle Tracer Explorers/Charge Composition Explorer (AMPTE/CCE) magnetic field, ion, and plasma wave data. They found that magnetospheric hot proton distributions, from 1 keV to 50 keV, are not well characterized by a single bi-Maxwellian distribution. By fitting a sum of several bi-Maxwellians to the data, Anderson et al. [1996a] improved the analytical fit to the observations, reducing the residual between the fit and the data by factor of 4 to 30. The conclusion of Anderson et al. is that determination of T⊥ and T∥ by a single bi-Maxwellian moment calculation is inadequate for EMIC wave instability analysis. As a result, the full RC distribution function on global magnetospheric spatial and temporal scales is required for a realistic assessment of wave excitation, and consequently, the overall wave effect on the RC populations. [15] It is shown above that use of a bi-Maxwellian fit to the RC H+ (and also to the RC O+) distribution can overestimate the local equatorial growth rate by an order of magnitude. The factor 10 was obtained using the Michigan RAM simulation, and re-evaluation using the full phase space particle distribution is sufficient to decrease the wave gain that was otherwise obtained by Thorne and Horne [1997] (the equatorial growth rate estimated from their results is about 6 × 10−1S−1) to values well below what is needed to account for wave growth during a single transit of the unstable region. [16] The above theoretical and observational evidences suggest that waves generally cannot grow substantially during a single transit of the unstable region [see also Demekhov, 2007], and some alternative generation models are needed [e.g., Guglielmi et al., 2001; Demekhov, 2007]. On the other hand, He+-mode EMIC waves are well guided along a magnetic field line and experience “fast” quasi-periodic bouncing between surfaces of the O+-He+ bi-ion hybrid frequency in opposite hemispheres [Horne and Thorne, 1993]. The plasmapause and/or dayside plume are the most favorable regions for EMIC wave generation in the inner magnetosphere [e.g., Horne and Thorne, 1993; Fraser et al., 2005] because net refraction is suppressed there. The wave normal angle in those regions oscillates about θ = π/2 slowly approaching that value. The ray path for the He+-mode in the vicinity of the plasmapause was illustrated in Figure A5 of Khazanov et al. [2006], and the “fast” and “slow” timescales were found to be τfast ∼ 102 s and τslow ∼ 103 s, respectively. Another timescale characterizing the wave evolution is a typical growth time and for geomagnetically active periods this was estimated to be τgrowth = 1/γ ∼ 103 s or slightly less. (Note that because waves cannot grow substantially during a single transit of the unstable region, the resulting γ includes both the energy source due to interaction with the hot RC and the energy sink due to absorption by thermal and hot plasmas and must be evaluated on a time scale of the wave bounce period [Khazanov et al., 2006].) The presented timescale hierarchy (along with the above theoretical and observational evidences) suggests that the bounce-averaged approximation employed in our RC-EMIC wave model is valid. So a physical model of EMIC wave bouncing between the off-equatorial magnetic latitudes corresponding to the ion-ion hybrid frequency in conjugate hemispheres, with tunneling across the reflection zones and subsequent strong absorption in the ionosphere, may be a potential candidate for the wave amplification model in the inner magnetosphere. [17] TH2007 state that wave observations contradict our modeling conditions. It is stated that EMIC waves are confined to a narrow wave normal angle in the equatorial active zone. TH2007 further state that directional wave observations contradict the existence of bouncing EMIC waves near the magnetic equator, while supporting their model for single pass wave growth with following severe damping near the bi-ion location. Our interpretation of in situ wave observations contradicts these assertions. We find that observational studies support the existence of bouncing wave packets at latitudes below the bi-ion frequency and unidirectional flow into the ionosphere above these latitudes. We find that EMIC waves routinely exhibit a mix of oblique and field-aligned wave normal angles near the magnetic equator, consistent with wave growth over multiple bounces between the bi-ion frequency latitudes. [18] Thorne and Horne [1997] suggested a model in which the wave normal angle is confined to less than 10 degrees over the entire near equatorial unstable region with subsequent severe wave damping near the bi-ion location. The restriction suggests that EMIC wave ellipticity is close to −1 in the entire near equatorial active zone. In other words, waves are near circularly left-handed polarized. This expectation contradicts observations [e.g., Fraser and Nguyen, 2001; Meredith et al., 2003; Loto'aniu et al., 2005; Anderson et al., 1992]. Observations find that wave events near the magnetic equator are evenly distributed from left-hand polarized to near linearly polarized with some right-hand polarized admixture, and there is a clear tendency for the polarization to become more linear with increasing magnetic latitude. Because it is not applicable in this case, the observation of a significant number of linearly polarized events near the equator cannot be explained by polarization reversal from left-handed through linear to right-handed at the crossover frequency, as discussed for other events by Young et al. [1981] (quasi-field aligned waves can have a linear polarization if the Young et al. mechanism takes place). Therefore the observed linear polarization suggests that waves will often be highly oblique inside the unstable region near the equator. Using the more reliable wave step polarization technique, Anderson et al. [1996b] and Denton et al. [1996] analyzed data from the AMPTE/CCE spacecraft and presented the first analysis of near linearly polarized waves for which the polarization properties have been determined. They indeed found a significant number of wave intervals with the wave normal angle θkB > 70°. The above observations cannot be reconciled with the wave amplification/damping model of Thorne and Horne [1997] (θkB < 10° in their scenario) but have a natural explanation in the framework of our RC-EMIC wave model. The new results presented by Khazanov et al. [2007] demonstrate that occurrences of the oblique and field-aligned wave normal angle distributions appear to be nearly equal near the magnetic equator with slight dominance of oblique events, consistent with observations. [19] Now we consider the TH2007 statement that the observational study by Loto'aniu et al. [2005] is consistent with the theoretical prediction of He+-mode growth and damping by Thorne and Horne [1997] and invalidates the concept of wave packet bouncing between off-equatorial magnetic latitudes corresponding to the ion-ion hybrid frequency. Loto'aniu et al. [2005] used magnetic and electric field data from CRRES to obtain the Poynting vector for Pc 1 EMIC waves. They reported bidirectional wave energy propagation, both away and toward the equator, for 26% of the events observed below 11° ∣MLat∣, and unidirectional energy propagation away from the equator for all events outside ±11° of the equator. Engebretson et al. [2005] found a similar EMIC wave energy propagation dependence but with mixed direction within approximately ±20° MLat, and consistently toward the ionosphere for higher magnetic latitudes. These observations lead Engebretson et al. [2007] to the conclusion that “the mixed directions observed in the above studies near the equator is evidence of wave reflection at the off-equatorial magnetic latitudes corresponding to the ion-ion hybrid frequency. Waves that reflect would then set up a standing (bidirectional) pattern in the equatorial magnetosphere. Waves that tunnel through would tend to be absorbed in the ionosphere and not be able to return to equatorial latitudes.” This conclusion by Engebretson et al. [2007] is in agreement with the physical picture underlying our RC-EMIC wave model [see also Guglielmi et al., 2001] and does not contradict, as will be shown next, the Loto'aniu et al. [2005] observations. [20] Loto'aniu et al. [2005] state that their statistical results show the unidirectional events outside ±11° of the equator. Upon close inspection of the published work, it is apparent that a majority of the unidirectional events are actually observed outside of ±18° of the equator because of data gaps between −18° and −14° and between −24° and −21° (the data gap in the northern hemisphere is an orbital effect). Loto'aniu et al. [2005] estimate the bi-ion frequency location at ∣MLat∣ ≈ 15°–20°, which is consistent with the 10°–20° ∣MLat∣ from [Rauch and Roux, 1982; Perraut et al., 1984]. Consequently, if there are heavy ions and waves are generated below the corresponding “bi-ion” latitude, they must be able to tunnel through the reflection zone to the latitudes observed (or pass through this zone freely if waves are guided). Although there are no concurrent observations, let us consider spectrograms 3a and 3b from Loto'aniu et al. [2005] as typical. We can see that high-latitude events have much less power than low-latitude ones. This is consistent with tunneling from a low-latitude source region to high latitudes through the region of bi-ion reflection. (Note that low-frequency events shown in Figure 3a of Loto'aniu et al. [2005] are likely generated at high latitudes.) The implication from those observations is that waves are not strongly damped before/after reflection contrary to a remark by TH2007. Inconsistency remains with identification of the transition latitude between bidirectional and unidirectional wave propagation in the two observational studies by Loto'aniu et al. [2005] and Engebretson et al. [2005]. At least partly, this inconsistency may be due to unavailability of wave observations at specific latitudes in the work of Loto'aniu et al. [2005] and/or differences in heavy ion content between the two studies. [21] Observations presented by Loto'aniu et al. [2005] below 11° ∣MLat∣ show that 26% of the events support the concept of wave packets bouncing between the off-equatorial magnetic latitudes corresponding to the ion-ion hybrid frequency. The events in Figure 3b of Loto'aniu et al. [2005] were observed at MLat ≈ −10.5°, that is, near the edge of the equatorial unstable region, and bidirectional wave energy propagation for packets b-h was observed. All these packets were mostly linearly polarized and, as a result, waves were highly oblique. As noted by Loto'aniu et al. [2005], on average, simultaneous compressional Pc 5 wave amplitudes were less than 0.3 nT over the EMIC wave events and it is unlikely that bidirectional pattern is due to a plasma property modulation by Pc 5. It is very difficult to generate highly oblique waves locally, and there is no active region below the satellite location. So the equatorward wave packets are likely reflected below the satellite at a latitude corresponding to the O+-He+ bi-ion frequency. If this reflection point is located well below MLat ≈ −10.5°, there is a conflict with the CRRES statistics because it did not observe the equatorially directed wave energy fluxes above 11° ∣MLat∣. However, as we pointed out above, this inconsistency may be due to unavailability of wave observations at specific latitudes in the work of Loto'aniu et al. [2005]. [22] An alternative explanation for not observing upward wave Poynting flux is because it is below the observational threshold. That explanation was suggested by Demekhov [2007]. Indeed, Loto'aniu et al. [2005] reported downward Poynting flux in the range of 1.3–10 μW/m2/Hz with an uncertainty of 0.1 μW/m2/Hz. The implication is that reflection of 1% of the wave energy would not have been observable. Given 1% reflection, a one-pass integrated gain Γ > 2.3 would be needed for the CRRES instrument to observe the reflected wave. This gain is quite realistic, and more sensitive measurements appear to be needed to confirm the presence or absence of reflected Pc 1 waves. Yet another explanation for not observing upward directed waves is because the upward wave is masked by the stronger down going wave [Demekhov, 2007]. Pearl wave packet durations are close to the repetition period, meaning that the packet length can be even longer than the field-line length between reflection points. In this case, the reflected portion of the wave packet may often be present at the same time as the large-amplitude downward waves. Off the equator and toward the reflection point the larger wave will dominate, obscuring observation of the weaker upward reflected wave. Near the equator wave amplification will take place along with a smaller wave group velocity and it will be more likely that waves in both directions will be observed with similar amplitude. Although a quantitative testing of these scenarios is required, the arguments presented by Demekhov [2007] seem reasonable. [23] First, let us comment on the statement by TH2007 that “An implicit assumption for the applicability of equation (1) is that after reflection, some wave energy is returned to the unstable region near the equator with propagation vectors aligned close to the ambient magnetic field direction to allow further amplification”. There is no implicit assumption for returning waves to the unstable region in our model. Khazanov et al. [2006, equation (22)] explicitly include the ray tracing equations in the RC-EMIC wave model. In other words, the effects of EMIC wave propagation and refraction are explicitly included in the equation which drives the wave power spectral density. The growth/damping rate in the right-hand side of the equation (22) is a result of averaging of the local growth/damping rates along the entire ray phase trajectory (r, θ) over the wave bounce period. The second term on the left-hand side of the equation is due to the fact that reflected waves return to the equator with more oblique wave normal angles, and it takes into account the wave energy outflow from the region of small wave normal angles to θ = 90°. The wave B-field distributions presented by Khazanov et al. [2006] are very well organized by the plasmapause and/or plume locations because the density gradient at those locations counteracts refraction caused by the magnetic field gradient and curvature. The net refraction is suppressed there, and after reflection, waves are able to return to the near equatorial unstable region with more or less field-aligned wave normal angles. This allows wave packets to spend more time in the phase region of amplification (not highly oblique wave normal angles in the near equatorial zone) resulting in wave growth dominantly at the plasmapause and/or in enhanced plasma density created by the dayside plume. This particular result is in complete agreement with the results of Horne and Thorne [1993] and Fraser et al. [2005]. [24] Second, it is stated in TH2007 that “Previous calculations of path-integrated wave gain during storm conditions [e.g., Thorne and Horne, 1997; Jordanova, 2005] are sufficient to drive waves to the observed amplitudes during propagation through the unstable equatorial region.” This suggests that RC H+ pitch angle diffusion will be changed from a regime of weak diffusion to a regime of strong diffusion during less than half of the EMIC wave bounce period that is about the same as the RC H+ bounce period. Such a change will cause strong modification of the RC distribution function. The approach promoted in TH2007 consequently makes Jordanova's bounce-average RC formalism inadequate to obtain the RC distribution function. While the bounce-averaged approach is inadequate to represent single-pass wave growth to high amplitude, it is entirely appropriate when growth to high amplitude, including nonlinear conditions, takes place on successive transits of EMIC waves through the equatorial unstable region. Under this circumstance, as presented in our publication, the RC phase space distribution function only slightly changes during the bounce period. [25] The main points of this Reply can be summarized as follows. [26] 1. The EMIC wave damping by thermal heavy ions and electrons, including the regions near the bi-ion location, have always been explicitly included in all our studies. [27] 2. The RC O+ can be neglected in the simulation of 2–7 May 1998 presented by Khazanov et al. [2006]. [28] 3. Our model is based on quasi-linear equations and the validity of this approach has been monitored in all versions of the model. Quasi-linear EMIC wave saturation takes place during most of the storm time, and the controlling criterion restricts wave amplitudes at values about 10 nT during the main and early recovery phases only when nonlinear saturation takes place. [29] 4. The insistence that wave growth takes place during a single transit of the unstable region is based on an approach [Thorne and Horne, 1997] that overestimates growth rates because of approximating particle distributions with a single bi-Maxwellian. It suggests confinement of the wave normal angle of propagating waves to less than 10 degrees, which contradicts observations near the magnetic equator. Rapid wave growth also violates the assumptions on which the bounce-averaged RC kinetic equation is based. [30] 5. The observation of EMIC waves by Loto'aniu et al. [2005] and Engebretson et al. [2005] at latitudes above the estimated reflection/tunneling points demonstrates that waves are not subject to severe damping before (and so after) they get reflected/tunneled. [31] 6. Contrary to the TH2007 statement that observational evidence contradicts our modeling results, the observations of Loto'aniu et al. [2005] and Engebretson et al. [2005] are consistent with our modeling. The explanation given by Engebretson et al. [2007] is for a physical model of EMIC wave bouncing between the locations of the ion-ion hybrid frequency at conjugate latitudes with tunneling across the reflection zones and subsequent strong absorption in the ionosphere. [32] To conclude, we welcome this discussion because it draws focus to the details of what is needed to accurately model the RC-EMIC wave processes in Geospace. The issues raised by TH2007 represent important differences in long standing published research that need to be resolved before the community can coherently advance in this field. We maintain the validity of the basic concept of the RC-EMIC wave model presented by Khazanov et al. [2006] through the discussion and evidence provided in this Reply. [33] This research was performed while K. Gamayunov held a NASA Postdoctoral Program appointment at NASA/MSFC. Funding in support of this study was provided by NASA grant UPN 370-16-10, NASA HQ POLAR project, and NASA LWS program. [34] Zuyin Pu thanks the reviewers for their assistance in evaluating this paper.