Chorus wave-normal statistics in the Earth’s radiation belts from ray tracing technique

. Discrete ELF/VLF (Extremely Low Fre-quency/Very Low Frequency) chorus emissions are one of the most intense electromagnetic plasma waves observed in radiation belts and in the outer terrestrial magnetosphere. These waves play a crucial role in the dynamics of radiation belts, and are responsible for the loss and the acceleration of energetic electrons. The objective of our study is to reconstruct the realistic distribution of chorus wave-normals in radiation belts for all magnetic latitudes. To achieve this aim, the data from the electric and magnetic ﬁeld measurements onboard Cluster satellite are used to determine the wave-vector distribution of the chorus signal around the equator region. Then the propagation of such a wave packet is modeled using three-dimensional ray tracing technique, which employs K. R¨onnmark’s WHAMP to solve hot plasma dispersion relation along the wave packet trajectory. The observed chorus wave distributions close to waves source are ﬁrst ﬁtted to form the initial conditions which then propagate numerically through the inner magnetosphere in the frame of the WKB approximation. Ray tracing technique allows one to reconstruct wave packet properties (electric and magnetic ﬁelds, width of the wave packet in k -space, etc.) along the propagation path. The calculations show the spatial spreading of the signal energy due to propagation in the inhomogeneous and anisotropic magnetized plasma. Comparison of wave-normal distribution obtained from ray tracing technique with Cluster observations up to 40 ◦ latitude demonstrates the reliability of our approach and applied numerical schemes.


Introduction
The dynamics of the Earth's outer radiation belt have been intensively studied since pioneering work in the late 1950s and early 1960s (see e.g. van Allen and Frank, 1959;Dungey, 1963), and yet their understanding is still limited. One of the most important processes that govern its dynamics is assumed to be the resonant wave-particle interaction. It is thought to contribute to electron acceleration (up to ∼ MeV energies) and losses into the ionosphere (Lyons and Thorne, 1973), especially during geomagnetic disturbances (see Horne, 2002;Chen et al., 2007;Shprits et al., 2008b). Since increases in the energetic electron flux are hazardous to satellites (Baker, 2002), manned space missions, or can drastically affect atmospheric chemistry (Thorne, 1977), it is crucial to improve the existing numerical models that specify and predict the dynamics of radiation belts (Bourdarie et al., 1996;Li et al., 2001;Glauert and Horne, 2005;Summers et al., 2007;Fok et al., 2008;Mourenas et al., 2012). These models are based on the quasi-linear approach where resonant wave-particle interactions are described in terms of particle pitch angle and energy diffusion (Kennel and Petschek, 1966;Trakhtengertz, 1966;Lyons, 1974;Lyons and Williams, 1984) to determine the time scale for electron losses and acceleration. Using this approach, it has been shown that whistler mode waves are one of the most efficient drivers of electron losses in radiation belts (e.g. Lyons et al., 1972;Albert, 2003;Horne et al., 2005b). In addition, electron resonant interaction with ULF waves also results in effective radial diffusion due to violation of third adiabatic invariant (see recent review by Shprits et al., 2008a), and radial transport can be enhanced by scattering of the second adiabatic invariant (see e.g. Ukhorskiy et al., 2011, and references therein).
However, since original paper by Lyons et al. (1971) and until now (see e.g. Glauert and Horne, 2005;Summers et al., 2007;Albert, 2007), the majority of calculations of the diffusion coefficients have been based on the assumption of parallel whistler wave propagation, relative to the background magnetic field. According to this approximation, the distribution functions of wave frequency and angle θ between the magnetic field and wave-vector are taken as Gaussian functions with two fixed parameters: mean value and a constant variance. These parameters can be different for various MLT, K p index and for several ranges of magnetic latitude (see e.g. Summers et al., 2007;Shprits, 2009;Ni et al., 2011, and references therein). Moreover, mean wave amplitude also varies with MLT (Horne et al., 2005a). Another assumption is that the wave power along the L-shell is constant (Glauert and Horne, 2005;Summers et al., 2007). These simplifications are used due to the lack of observational data of these distributions.
An analysis of chorus whistler wave data of Cluster satellites during the 2001-2009 years has been performed by Agapitov et al. (2011aAgapitov et al. ( , 2012, where the measurements were concentrated in the outer radiation belt, i.e. between 4R E and 7R E (here R E is the Earth's radius). The major result of this study is the evidence of the strong dependence of characteristics of θ-angle distribution upon magnetic latitude. Indeed, as the latitude increases, the mean value of the distribution, which is around 20 • at the magnetic equator, shifts towards larger angles and the variance (width) of the distribution also increases, confirming the tendencies observed by Hayakawa et al. (1986) and reported in the comprehensive review by Sazhin and Hayakawa (1992). Moreover, these results were confirmed near the equator by Li et al. (2011) by use of THEMIS spacecraft statistics. Thus, providing the proper model that allows us to supplement existing databases of wave distributions on frequency, wave vector and L-shell, the accuracy of the codes that compute the radiation belts diffusion coefficients can be significantly improved (Horne et al., 2005a;Shprits and Ni, 2009;Artemyev et al., 2012a,b). To this end, it is necessary to develop a numerical ray tracing code that can first reproduce wave distributions based on direct measurements and then to fill up the observed distributions by "complementary data" obtained from simulations.
The purpose of this paper is to model the distribution of chorus wave parameters, taking as a basis Cluster observations and completing them by results of numerical modeling, presented hereafter.
In the following sections, we shall describe the numerical code used in this study: the key parameters that are needed, and the models of the magnetic field and plasma density distributions inside the inner magnetosphere that are included. Then we present the set of parameters necessary for computation of ray traces. The results of calculations form the database that is then used to calculate statistical distributions of whistler wave-normal angles at different magnetic latitudes for different L-shells. Finally, the comparison of the obtained numerical k-vector distribution with direct observations is presented, followed by conclusions. The description of ray-tracing methodology and numerical code can be found in Appendix.

Numerical model description
In this section, we present the numerical model that computes ray trajectories of chorus emissions (ELF/VLF whistler modes) in the inner magnetosphere. The calculation is carried out making use of realistic plasma density model and pre-selected initial wave distributions. The source localization of the chorus waves is assumed to lie at the magnetic equator between 4.5R E and 7R E , according to LeDocq et al. (1998);Santolík et al. (2005a,b); Agapitov et al. (2011a); Li et al. (2011). By use of the code, we integrate wave trajectories through the Earth's inner magnetosphere, in the frame of the geometrical optics approach (i.e. ray tracing; see Appendix for more details). This code is valid for any linear electromagnetic wave mode in the anisotropic multicomponent (up to 6 species) hot plasma with slowly varying (in comparison with the wavelength) spatial parameters and moderate absorption. The presented code employs WHAMP (Waves in Homogeneous Anisotropic Multicomponent Plasma) code developed by Rönnmark (1982) to solve hot plasma dispersion relation, whose solution is then used to trace whistler ray trajectories.
This dispersion solver allows to improve the evaluation of wave characteristics such as wave damping. In particular, one can better describe the effect of absorption by the medium, which can result in the improvement of the calculations at lower altitudes.
Additional effects that can also be taken into account are specific features of particle distributions that may be present in radiation belts such as loss cone particle distribution, temperature anisotropy or relative drifts of different particle populations. One of the important available options is the possibility to include suprathermal energetic particles that play an important role in wave generation.
Magnetospheric chorus emissions are commonly observed between 2.5R E and 10R E . Near 2.5R E the Earth's core magnetic field is assumed to have a tilted dipolar structure with the dipole magnetic momentum corresponding to year 2005 epoch, and an analytical model (Olson and Pfitzer, 1977) includes the contribution of sources external to the Earth Ann. Geophys., 30, 1223-1233, 2012 www.ann-geophys.net/30/1223/2012/ (namely magnetopause, tail and ring currents) at larger altitudes. It is valid for all tilts of the Earth's dipole axis during rather quiet magnetosphere (K p ≤ 4), and has been optimized for the near-Earth region (from 2R E to 15R E ). We use the Global Core Plasma Model (Gallagher et al., 2000), denoted GCPM hereafter, to obtain smoothly varying densities of the magnetospheric plasma species throughout the entire volume of the inner magnetosphere. GCPM is an empirically derived core plasma density model of four plasma components -electrons, protons, oxygen ions and helium ions. It consists of separate models for the plasmasphere, plasmapause and polar cap, and merges with IRI (International Reference Ionosphere) 2007 model at low altitudes (Bilitza and Reinisch, 2008). Moreover, this model supplies complete three-dimensional spatial density distributions. It is important to note that, whereas the plasma density model we use is three-dimensional, we reduce it to a 2-D grid in radial distance and latitude for one Magnetic Local Time (MLT = 09:00); i.e. we neglect azimuthal variations. This procedure allows us to significantly reduce the computation time.
Treating numerically magnetospheric chorus wave propagation, it is important to describe properly the region of plasmapause, a very dynamic transition region where plasma density drops abruptly. In the plane of magnetic equator, the plasmapause is situated at the distance about 4R E for rather quiet magnetic conditions (K p = 4 here). Indeed, the steep gradients of density greatly affect the trajectories of the propagating waves, and thus need to be well defined by the code. Figure 1 shows the electron density as a function of radial distance for the set of magnetic latitudes, at MLT = 09:00, which is rather typical for chorus wave registration for a level of magnetospheric activity of K p = 4 (see Agapitov et al., 2011a). The plasmapause is well defined here, as it is smooth and clearly present at low latitudes.
To complete the plasma model, one should define the characteristics of the particle distribution functions. In our calculations, we assume that all plasma species obey Maxwell's distributions with temperatures of 0.5 eV, which approximatively correspond to the Akebono data-based temperature model (Kutiev et al., 2002).

Reconstruction of chorus wave-vector distribution upon magnetic latitude
The main goal of our paper is the reconstruction of chorus wave-normal direction distribution in the inner magnetosphere in such a way that it can be compared with Cluster spacecraft observations. To do so, chorus wave power distribution is to be defined realistically in the source region supposed to be at equator. Then their trajectories are calculated by means of the ray tracing technique (see Appendix), using the realistic model of the inner magnetosphere described above. In this section, we show how the statistics of cho- rus wave-normal angle distribution as a function of latitude λ observed by Cluster can be reproduced using numerical calculations of ray propagation from the chorus source region.
Chorus waves typically appear as short coherent bursts of wave packets propagating in two distinct bands, scaling on the equatorial electron gyrofrequency e,equ . The lower band frequency range is known to be ω = 0.1−0.45 e,equ , and the upper band frequency range is ω = 0.5 − 0.7 e,equ (see e.g. Helliwell, 1969, 1976). However, here we present only results for the lower-band chorus because the Cluster STAFF-SA frequency range cannot completely cover the full frequency range of the upper band (Agapitov et al., 2011a). Additionally, upper-band choruses are substantially less intense than lower-band chorus (Meredith et al., 2001;Haque et al., 2010). The chorus source region is located in the vicinity of the magnetic equator (see also LeDocq et al., 1998;Santolík et al., 2005a,b) (sometimes the source can be also situated in non-equatorial minimum B-pockets during high magnetospheric activity (Tsurutani and Smith, 1977;Agapitov et al., 2010), but more detailed discussion of this topic is beyond the scope of the present study). Chorus waves are usually observed in the night, day and dawn sectors (Agapitov et al., 2011a;Meredith et al., 2003;Li et al., 2009) with a maximum in occurrence probability and wave intensity between 06:00 and 12:00 MLT (see e.g. Li et al., 2009;Agapitov et al., 2011a). Moreover, chorus waves are known to be generated outside the plasmasphere with the amplitude www.ann-geophys.net/30/1223/2012/ Ann. Geophys., 30, 1223-1233, 2012 maximum at L ≈ 7 (Meredith et al., 2003;Santolík et al., 2005b;Li et al., 2009). Finally, as mentioned above, recent studies of chorus wave-normal distributions using extensive sets of observations from different spacecraft (see Agapitov et al., 2011a;Li et al., 2011) have shown that the distributions of intense chorus waves (i.e. lower-band, rising elements) are generally non-Gaussian and have a non-zero mean value. In order to reconstruct chorus wave-normal distribution statistics, we have calculated the numerous ray trajectories in the inner magnetosphere (in terms of L-shell, latitude λ and longitude φ). Then these distributions of rays have been properly weighted in terms of L, θ and frequency at the magnetic equator, to accurately reproduce their statistical characteristics in the source region described above. The wave power distribution is described by the function G(f, θ, ϕ), where θ and ϕ are the latitudinal and azimuthal angles that describe the orientation of wave-normal in the field-aligned coordinate system, respectively, at a given frequency f (see wave distribution function (WDF) method introduced by Storey and Lefeuvre, 1974).
Making use of the code described in Appendix, we have computed a database consisting of ∼ 35 000 numerical trajectories of electromagnetic whistler waves generated in the vicinity of the magnetic equator and magnetospherically reflected at high latitudes.
The wave parameters included in this database are as follows: plasma density parameters have been chosen to reproduce the typical region of chorus occurrence (MLT = 09:00) for a level of magnetospheric disturbance, during which maximum chorus activity is observed (K p = 4.0). The date is 7 September 2002, at 00:30 UT. Rays are generated at the equator within the range 4.5 ≤ L 0 ≤ 7, with a step of L 0 = 0.1. At each starting point L 0 , we launch a set of waves with frequencies in the range 0.1 e,equ ≤ ω ≤ 0.5 e,equ , spaced at intervals of ω = 0.1 e,equ , to reproduce the lower-band chorus wave power distribution. The 3-D obliqueness of chorus wave propagation is represented by the two angles θ and ϕ, with θ being the latitudinal angle between k-vector and magnetic field, and ϕ the azimuthal angle. For each L 0 and ω, the range of initial angles ϕ 0 between k-vector and the background magnetic field B 0 in the equatorial XY plane is −180 • ≤ ϕ 0 ≤ +180 • with ϕ 0 = 10 • , ϕ 0 = 0 • pointing towards the Earth. Since the resonance cone angle θ res is frequency-dependent, the range 0 ≤ θ 0 ≤ θ res of initial angles θ 0 between k-vector and B 0 in the XZ plane is non-uniform, varying from 0 ≤ θ 0 ≤ 80 • for a frequency of 0.1 e,equ , and 0 ≤ θ 0 ≤ 55 • for 0.5 e,equ . The interval step between each θ 0 is θ 0 = 5 • .
This numerical ray distribution is uniform in terms of ω and also in terms of L 0 and θ 0 for chorus wave frequency range. It thus needs to be intensity-weighted at the injection points to represent the chorus source region. This weighting is divided into three individual components representing the dependence on L 0 , θ 0 and ω. Chorus wave power distribution in frequency is usually approximated by a Gaussian function with a peak value of ω ≈ 0.34 e,equ (see Burtis and Helliwell, 1976). Agapitov et al. (2011a) have shown that chorus waves have a distribution close to a Gaussian over L-shells for moderate activity in the inner magnetosphere, with a maximum value at about L 0 = 6.5. Finally, according to Agapitov et al. (2011a), the equatorial θ 0 distribution is not uniform but has a pronounced maximum at θ 0 ≈ 20 • .
Thus, at the magnetic equator, we use a weight function h 0 (θ 0 ) for θ 0 distribution which is simply a cross section of the experimental distribution in Fig. 2e from Agapitov et al. (2011a) As a result, each ray trajectory is weighted with an intensity given by the following weight function: where the subscript 0 denotes the initial (equatorial) value. It is worth noting that, for all ray trajectories, the selected region in space is situated in the range 4.5 < L < 7 at every latitude; thus, ray coordinates leaving it or entering from outside this L-shell range are not taken into account in the simulations.
In Fig. 2, one can see that the poleward ray wave-normals quickly tend to the resonance cone; the maximum value of the distribution starts near the equator at around 20 • and reaches nearly 80 • before 30 • latitude. Another important characteristic of the distribution is that the variance rests approximately constant at all latitudes up to 30 • . The behavior of the numerical distribution of poleward rays is thus very Ann. Geophys., 30, 1223-1233, 2012 www.ann-geophys.net/30/1223/2012/ Fig. 3. Spatial divergence of ray trajectories as a function of latitude, with respect to wave frequency. Each ray position intersecting one latitudinal plane (λ = const) is plotted by a point, whose color represents its frequency, from ω = 0.1 e,equ (dark blue) to ω = 0.5 e,equ (red) with ω = 0.1 e,equ . Rays are launched from a single starting point of coordinates [5R E , 0, 0], confined in a cone centered around the local magnetic field, with parameters 5 • ≤ θ 0 ≤ θ res and −180 • ≤ ϕ 0 ≤ +180 • . consistent with the observed distribution presented in Fig. 2e from Agapitov et al. (2011a), where the same wave-normal angle distribution is constructed using the Cluster STAFF-SA measurements from years 2001 to 2009. In the following, we use X = tan θ as parameter since it is widely used as an input for diffusion coefficient calculation (Lyons et al., 1971). The good agreement is confirmed by Fig. 4, where the probability distribution function (PDF) of X distribution of poleward rays as a function of latitude is presented. Indeed, the distributions obtained from experimental data (left panel) and as a result of numerical simulations (right panel) both exhibit the same tendency, i.e. a rapid increase of the mean value and variance with the growth of the latitude. The quantitative growth of these two parameters is discussed in the next section. It is worth noting here that the distributions are shown only up to 30 • latitude, since the experimental data for higher latitudes are not sufficient to perform reliable statistical analysis (this can be seen from the large error bars in Fig. 5 in the range 30 • ≤ λ ≤ 40 • ). However, from Fig. 2 one can notice some discrepancies with the k-vector distribution from Agapitov et al. (2011a) for equatorward rays, for which 90 • < θ ≤ 180 • . Indeed, although the variance and shape of the distribution are quite similar to those observed, the maximum of equatorial distribution for numerical rays is centered at θ ≈ 120 • , whereas according to Cluster data the peak is supposed to stand at θ ≈ 155 • . These rays can either be magnetospherically reflected chorus waves generated at the magnetic equator, or chorus waves generated from non-equatorial minimum-B pockets (which might not be taken into consideration in our simulations). These discrepancies can be explained by the fact that chorus waves are spatially deviated from their initial L-shell during their propagation and magnetospheric reflection in the magnetosphere, resulting in power attenuation of reflected signals (see e.g. Parrot et al., 2003Parrot et al., , 2004Agapitov et al., 2011b). To illustrate these effects, we show in Fig. 3 the spatial spreading of different ray trajectories for the plasma parameters mentioned above. We launch, from a single starting point of coordinates [5R E , 0, 0], a set of rays towards the Northern Hemisphere. These rays are confined in a cone centered around the local magnetic field with parameters 5 • ≤ θ 0 ≤ θ res and −180 • ≤ ϕ 0 ≤ 180 • . We vary wave frequency from 0.1 e,equ up to 0.5 e,equ with step 0.1 e,equ . The different wave frequencies are shown by color on the Fig. 3. The distribution of ray trajectories is symmetrical with respect to Y coordinates, since we neglect the azimuthal density variations. The ray trajectories are plotted in a different panel for each latitudinal plane up to λ = 50 • . The higher frequency waves are seen to reach lower X values faster than lower frequency waves, and are more concentrated around Y = 0, whereas waves with frequency of 0.1 e,equ are more regularly distributed along Y, forming a ring-like pattern as latitude increases. One can see here that the spatial spreading of poleward ray trajectories is very large. Thus, a coherent chorus wave packet (consisting of different frequencies) would greatly diverge as it propagates towards higher latitudes. Therefore, poleward or equatorward waves that originate from outside the range 4.5 ≤ L ≤ 7 can be included in statistical studies based on spacecraft measurements. However, the problem of geometrical spreading and magnetospheric reflection of chorus wave packets (with respect to plasma and magnetic field parameters for instance) requires a detailed study and will be addressed in future work.
This good statistical agreement between numerical and observed behavior of poleward chorus k-vector distributions corroborates the possibility of reconstructing the observed distribution upon all magnetic latitudes making use of our numerical database of ray trajectories.

Discussion
In our model, we have made use of realistic properties of the chorus wave distribution in the source region (Agapitov et al., 2011a;Li et al., 2009Li et al., , 2011 to obtain wave distributions in the inner magnetosphere by means of ray tracing. We have found that the distributions obtained have very similar statistical properties as observed onboard spacecraft in case studies (see e.g. Agapitov et al., 2011b), in statistical studies of STAFF-SA wave data measurements of Cluster armada , and similar to results of other simulations (see e.g. Bortnik et al., 2011). We have shown that our calculated distribution of chorus wave-normals is very similar to the characteristics obtained by Agapitov et al. (2011a). It is worth noting that the initial distribution of chorus wave normals in the vicinity of the equator is also close to the distributions often used in the numerical studies of particle diffusion. In previous calculations of energetic electron angular diffusion rates (Horne et al., 2005a;Glauert and Horne, 2005;Summers et al., 2007;Shprits and Ni, 2009;Artemyev et al., 2012a), these distributions are described in terms of Gaussian ansatz of the variable X = tan θ as where X m and X w are the mean value and the variance of the distribution, respectively. In most of the above-cited modelling works, X m = 0 and X w was taken as a constant. Our numerical model shows that the average value of the angle θ grows with the latitude that is in a good agreement with observations onboard Cluster satellite (Agapitov et al., 2011a), and shows that, at medium and high latitudes, the approximation based on constant value of the average angle becomes invalid. This difference in wave-normal distributions is known to be very important for the evaluation of the angular diffusion rates and can lead to an underestimation of the diffusion rates especially for electrons with small pitch angles, which spend a large portion of their time oscillating along the field Ann. Geophys., 30, 1223-1233 line at high latitudes (Artemyev et al., 2012a,b;Mourenas et al., 2012). The distribution of X obtained in our calculation is very close (see Fig. 4) to the results obtained from Cluster STAFF-SA measurements. The variations of parameters X m and X w from Cluster STAFF-SA and from our numerical calculation are shown in Fig. 5 (indicated by circles and by solid lines respectively) up to large (40 • ) magnetic latitudes. The weighting function defined by Eq. (1) is here applied to the numerical database of initial equatorial ray distribution directed poleward in order to reproduce observed distributions at equatorial latitudes. Figure 5 shows that the agreement for the mean value for all magnetic latitudes (top panel) is very good. The tendency for the variance (bottom panel) is also well reproduced as seen in this figure, where the global increase of the width of distribution with magnetic latitude is quite clear. However, some of the features of the variance seen in experimental data, such as the rapid growth up to ≈ 12 • latitude and the decrease at ≈ 25 • , are not reproduced in the simulations. There can be several reasons for this; for instance, it can be due to the fact that experimental data are averaged over dawn/day sector (02 < MLT< 14), and for a time period of ten years, during which the plasma density parameters change due to magnetospheric activity or due to the presence of waves that are generated on different L-shells, or due to some other reasons. We plan to carry out more studies to resolve this question in the future. Since chorus wave trajectories are very sensitive to parameters mentioned above, this probably may explain the difference of the dispersion of wave-normal angles with latitude as observed in experimental data, compared to numerical computations. For higher latitudes, this can also be explained by the lack of experimental data, as already mentioned above. From this comparison, the statistical distributions obtained from numerical computations using our database seem to be a reliable estimation for chorus wave parameter distribution dependence upon magnetic latitude. This complementary data could help to estimate parameters in the regions where observational statistics is not sufficient, for example at high latitudes, as seen on the large error bars from 30 • latitude in Fig. 5. www.ann-geophys.net/30/1223/2012/ Ann. Geophys., 30, 1223-1233, 2012 In this paper, we present a new 3-D ray tracing code for a hot magnetized plasma with moderate absorption, which has been developed including realistic models of plasma densities and magnetic field for the inner magnetospheric region. Ray tracing was carried out assuming that the wave source is situated at equator. The initial distribution of waves in the source region was chosen to correspond to statistical distribution dependencies on k-vectors, L-shells and frequencies obtained from observations. To this end, the weight functions corresponding to distributions inferred from observations were applied to initial set of rays. Then, making use of our numerical database, we reconstruct chorus wave-vector distributions as a function of magnetic latitude using weight functions as described above. The results of our calculations are in good agreement with statistical distributions found using ten years of observational data measured onboard Cluster spacecraft. The distribution notably exhibits large wavenormal angles at medium and high latitudes, which makes the quasi-longitudinal approximation no longer valid at these latitudes. This can lead to rather large errors in the evaluation of the diffusion rates that define particle losses very important for correct description of the radiation belt dynamics. Thus, our numerically obtained statistical database is quite important as a reliable starting point for chorus wave distributions reconstruction in the magnetospheric regions where the satellite data are poor. It can be used to fulfill the lack of observations by providing necessary parameters for the global simulations, in particular, at higher latitudes, where the effects of resonant interactions between waves and particles can play very important role determining the radiation belts dynamics (see Shklyar and Matsumoto, 2009;Artemyev et al., 2012a,b, for instance).

Ray tracing technique and numerical code description
Linear oscillations in a magneto-active plasma are determined by the linearized Vlasov equation together with the Maxwell equations for the variable fields E and B (for calculation details, see e.g. Sagdeev and Shafranov, 1961;Ginzburg, 1970;Akhiezer, 1975). Combining these two equations, it is possible to calculate the dielectric permittivity tensor and then solve the well-known dispersion relation (see Rönnmark, 1982, for instance).
To follow the trajectory of a given wave in the plasma, one should calculate wave parameters at each propagation point. This can be achieved making use of real Hamiltonian equations in a medium with moderate absorption (Suchy, 1981), in other words, by using eikonal equations which allow one to express the Maxwell equations for electromagnetic waves in terms of geometrical rays. This approximation is called the WKB or geometrical optics approximation and usually leads to the system of ordinary differential equations (ODEs) that can be then integrated numerically (e.g. see Haselgrove, 1954;Yabroff, 1961;Kimura, 1966;Inan and Bell, 1977;Shklyar and Jiříček, 2000;Bortnik et al., 2011).
The electromagnetic field should be presented in the following form: for all six components, where r is the position, t is the time, k is the wave vector and ω the wave frequency. Then, assuming that the field amplitude E{ B}(k, ω; r, t) varies slowly compared with the eikonal (written as a four-dimensional line integral), Then Maxwell's equations together with the constitutive equations are transformed into a system of eikonal equations.
First, the condition of solvability is nothing other than the dispersion relation Then, a solution of this equation by the method of characteristics (Suchy, 1981), assuming moderate absorption, leads to the following ODE system (real Hamilton equations for complex space-time coefficients): where A is the amplitude of the wave.
In the framework of the geometrical optics approach, density and magnetic field strength are assumed to be slowly varying parameters compared to wavelength scale. Thus, at any point of the considered plasma system, it is possible to construct the dispersion relation with complex coefficients slowly varying in space and time. Then the wave path and Ann. Geophys., 30, 1223-1233, 2012 www.ann-geophys.net/30/1223/2012/ the parameters along propagation trajectory can be reconstructed by integrating the ODE system (Eq. A4) by any standard method.
Having core plasma parameters and magnetic field distribution completely resolved in any point of the modeled volume (in our case 2-D (L, λ) grid), one is able to define the wave dispersion, and to obtain all the characteristics of the propagating whistler wave by solving the arising system of ordinary differential equations (Eq. A4). These are, namely, the points r of its trajectory, as well as the wavenumber k, the complex frequency ω, and the amplitude A in all points of the wave path.
Once the densities and the magnetic field have been tabulated on the grid, the wave packet has to be initialized; i.e. dispersion function has to be calculated in the equatorial "starting" point of the ray. In this study, the source is supposed to be on the magnetic equator and to have finite perpendicular scale in the equatorial plane.
To define the wave dispersion, we use an algorithm implemented in the well-known WHAMP code (Rönnmark, 1982), which calculates the dispersion function of the hot magnetized plasma.
Up to six different plasma components can be included into the code, each being specified by its density, temperature, particle mass, drift velocity along the magnetic field, the depth and the size of the loss cone and temperature anisotropy. Thus, the program is applicable to a really wide class of plasmas, and the method should in general be useful whenever a homogeneous magnetized plasma can be approximated by a linear combination of Maxwellian components.
The kinetic dielectric tensor is approximated by introducing a Padé approximant for the plasma dispersion function (Fried and Conte, 1961) when evaluating the susceptibility tensor. The resulting expression for the dielectric tensor is valid for all real k and very well suited for numerical evaluation. Then, a complex frequency ω, which satisfies the dispersion relation (Eq. A3), can be found using any iteration method. Here Newton's iteration method is employed (see e.g. Press et al., 1992). After each wave packet has been initialized (the frequencies and wavenumbers of chorus forming whistlers are calculated consistently for the plasma and magnetic field parameters in the supposed source region), its propagation trajectory in the magnetosphere is traced by integrating the ODE system (Eq. A4) by Runge-Kutta-Verner fifth-order initial value ODE solver.