Kinetic simulation of the electroncyclotron maser instability: relaxation of electron horseshoe distributions
^{1} Armagh Observatory, Armagh BT61 9DG, Northern Ireland
email: aku@arm.ac.uk
^{2} Institute of SolarTerrestrial Physics, 664033 Irkutsk, Russia
Received: 15 September 2010
Accepted: 17 November 2010
Context. The electroncyclotron maser instability is responsible for the generation of the auroral kilometric radiation of the Earth and similar phenomena at other magnetized planets of the Solar System. The recently discovered radio emission from ultracool dwarfs has many similarities with the planetary auroral radio emission. The in situ measurements in the terrestrial magnetosphere indicate that the radiation is produced by nonthermal electrons with a horseshoelike distribution. Kinetic simulations of the electroncyclotron maser instability for these distribitions have not yet been performed.
Aims. We investigate the amplification of plasma waves by the horseshoelike electron distribution, as well as the relaxation of this distribution due to the electroncyclotron maser instability. We determine the parameters of the generated plasma waves, the timescales of the relaxation process, and the conversion efficiency of the particle energy into waves.
Methods. We developed a kinetic relativistic quasilinear 2D code to simulate the coevolution of an electron distribution and highfrequency plasma waves. The code includes the processes of wave growth and particle diffusion, which are assumed to be much faster than other processes (particle injection, etc.). A number of simulations were performed for different parameter sets that seem to be typical of the magnetospheres of ultracool dwarfs (in particular, the plasma frequency is much less than the cyclotron one).
Results. Our calculations indicate that the fundamental extraordinary mode dominates strongly. The generated waves have a frequency that is slightly below the electron cyclotron frequency and propagate across the magnetic field. The final intensities of other modes are negligible. The conversion efficiency of the electron energy into the extraordinary waves is typically around 10%. Complete relaxation of the unstable electron distribution takes much less than a second.
Conclusions. Energy efficiency of the electroncyclotron maser instability is more than sufficient to provide the observed intensity of radio emission from ultracool dwarfs. On the other hand, the observed light curves of the emission are not related to the properties of this instability and reflect, most likely, the dynamics of the electron acceleration process and/or geometry of the radiation source.
Key words: radiation mechanisms: nonthermal / planets and satellites: aurorae / brown dwarfs / radio continuum: stars / radio continuum: planetary systems
© ESO, 2011
1. Introduction
The electroncyclotron maser instability (ECMI) can occur in a magnetized plasma with a nonequilibrium electron distribution function (a positive slope in perpendicular velocity is required) and results in an effective amplification of electromagnetic waves at the low harmonics of the electron cyclotron frequency (Wu & Lee 1979; Melrose & Dulk 1982; Winglee & Dulk 1986). In a relatively rarefied plasma (when the electron plasma frequency to cyclotron frequency ratio ω_{p}/ω_{B} ≪ 1), the dominant mode of astrophysical electroncyclotron masers is the extraordinary wave with a frequency close to the fundamental cyclotron frequency. The emission has a narrow spectral bandwidth, a narrow directivity pattern, and high (nearly 100%) degree of circular polarization. ECMI is responsible for the generation of the auroral kilometric radiation (AKR) in the magnetosphere of the Earth and, most likely, the similar phenomena at other magnetized planets of the Solar System (Zarka 1998; Treumann 2006). ECMI has also been applied to interpret certain types of solar and stellar sporadic radio bursts (Melrose & Dulk 1982; Fleishman & Melnikov 1998).
Recently, a number of late M stars and brown dwarfs (termed ultracool dwarfs, UCDs) were found to be the sources of unexpectedly intense radio emission at the frequencies of a few GHz (Berger 2005; Hallinan et al. 2006, 2007, 2008; Antonova et al. 2008). The emission includes a slowlyvarying weaklypolarized component, as well as short intense periodic bursts with almost 100% polarization, whose period seems to coincide with the rotational period of the star. While the slowlyvarying component may be explained by the incoherent gyrosynchrotron radiation, the periodic bursts require a coherent radiation mechanism such as ECMI. Measurements of magnetic fields for the stars of late spectral class (M4) using phaseresolved spectropolarimetry (Donati et al. 2006, Morin et al. 2008) indicate that the magnetic field has a dipolelike structure where the dipole axis is close to the rotational axis. Thus, the radio emission of UCDs seems to be similar to the auroral emission of the Solar System planets, but with a much stronger magnetic field (not less than 3000 G, to provide the highest observed emission frequency). The existence of such magnetic fields at cool stars (<M 9) is also confirmed by infrared measurements (e.g., Reiners & Basri 2007).
ECMI has been investigated analytically and using numerical simulations in a number of papers (see, e.g., the review of Treumann 2006, and references therein). The required unstable electron distributions (with a deficiency of particles with low transversal velocity) are formed naturally due to particle reflection from a magnetic field gradient. Using a linear approximation, it is easy to show that even relatively weak fluxes of accelerated electrons can result in large growth rates, so the waves can be amplified from the level of thermal fluctuations to the observed intensities at distances not exceeding a few kilometers (see, e.g., Melrose & Dulk 1982; Bingham & Cairns 2000; Bingham et al. 2001, as well as estimates in this article). However, the resulting wave energy is expected to be comparable to the particle energy; this leads to considerable particle diffusion on the waves, which, in turn, ensures that the linear approximation is no longer applicable. To interpret the cosmic radio emissions in terms of ECMI, we have to use a nonlinear model that takes into account the relaxation of the unstable electron distribution due to interaction with the excited waves. In particular, only nonlinear models can provide us with such an important parameter as the transformation coefficient of the particle energy into waves.
In early studies, ECMI was associated with the electron distribution of the losscone type. Kinetic simulations of the relaxation of the losscone were presented by Aschwanden (1990) and Fleishman & Arzner (2000). However, in situ measurements in the AKR sources have shown that the radiation should be produced (at least, in this particular case) by a ringlike or horseshoelike distribution (e.g., Delory et al. 1998; Ergun et al. 2000). These distributions are formed when electron beams (accelerated by an electric field) move into stronger magnetic field regions. The kinetic simulation of the relaxation of these distributions have not yet been done (although there are some particleincell simulations; see, e.g., the references in the review of Treumann 2006).
When simulating the relaxation of unstable electron distributions, the greatest challenge is to account for the spatial movement of the plasma waves and particles. Wave propagation out of the region occupied by the electrons with an unstable distribution naturally limits the wave amplification. On the other hand, the waves amplified in one part of the radiation source will cause relaxation of the electron distribution in other parts. Electron movement in an inhomogeneous magnetic field is necessary to form an unstable distribution (of the losscone or horseshoe type). To take into account all the aforementioned factors, one has to create a 3D model of the emission source and adjacent regions (e.g., a magnetic loop in the solar or stellar corona), which, in turn, requires enormous computational resources. The most popular solution of this problem is to neglect the spatial movement of waves and particles completely, so that the model is reduced to an initial value problem where an arbitrary unstable electron distribution is used as the initial condition. We can neglect spatial movements if the time of wave/particle escaping from the wave amplification region τ_{esc} far exceeds the typical diffusion time τ_{diff}. This approach (a diffusive limit) is used in the abovementioned papers on kinetic simulation, as well as in most particleincell simulations. This model is obviously far from reality since it cannot explain the longterm generation of the emission and, in addition, requires the almost instant formation of an unstable electron distribution. Nevertheless, a model that considers only the wave growth/damping and particle diffusion allows us to (i) investigate the qualitative behaviour of the relaxation of unstable electron distributions; (ii) determine the parameters of the produced waves; (iii) estimate the diffusion and relaxation timescales (which is necessary to make a conclusion about the validity of the diffusive limit); and (iv) determine the transformation coefficient of the accelerated particles energy into waves.
Measurements within the AKR sources have shown that the cold electron component is almost absent there. Under these conditions, the wave dispersion is determined mainly by the energetic electrons and the relativistic effects become important. Relativistic corrections to the dispersion relation result in a decrease in the cutoff frequency of the fast extraordinary mode; in addition, the fast and slow branches of the extraordinary mode can reconnect to form a single branch (Winglee 1983, 1985; Strangeway 1985, 1986; Robinson 1986, 1987; Le Quéau & Louarn 1989; Louarn & Le Quéau 1996). These effects allow the waves generated at the frequencies below the nonrelativistic cyclotron frequency to escape freely from the source into a vacuum. Solving the exact relativistic dispersion equation is a complicated task. However, previous studies (e.g., Robinson 1986, 1987) have demonstrated that in a sufficiently hot lowdensity plasma, the wave dispersion is similar to that in a vacuum (with the refraction index N → 1 and group speed 3_{gr} → c). This approximation seems to be valid for the sources of planetary auroral radio emission.
We developed a kinetic relativistic 2D code to simulate the coevolution of electron distributions and plasma waves in the diffusive limit. Different wave modes can be considered both separately and simultaneously. An electron distribution of the horseshoe type is considered to be the source of plasma oscillations. The calculations are made for conditions that seem to be typical of UCDs, although the results can be easily scaled to other emission sources (e.g., the Earth or Jupiter). Exact plasma parameters for the magnetospheres of UCDs are unknown; in particular, we do not know whether the cold plasma component exists or not. Therefore we consider two opposite cases: (i) when a lowtemperature plasma with a Maxwellian distribution dominates and (ii) when such a component is absent. In the former case, we use the cold plasma dispersion relation; in the latter case, the dispersion relation is assumed to be similar to that in a vacuum.
We note that the model used is restricted to the highfrequency elecromagnetic/magnetoionic waves and does not consider generation of the lowfrequency waves (e.g., acoustic ones). In the AKR sources, the amplitude of the lowfrequency waves can reach high levels resulting in the formation of solitary structures (such as electron and ion holes), which, in turn, can affect the generation of the radio emission (Pottelette et al. 2001; Mutel et al. 2006, 2007). However, these effects are beyond the scope of this paper.
The model used is described in Sect. 2. The initial conditions of the model (including the electron distribution function) are described in Sect. 3. The simulation results are presented in Sect. 4 and discussed in Sect. 5. The conclusions are drawn in Sect. 6. The calculation formulae (most of which can be found elsewhere) are given in Appendices.
2. Kinetic equations
If spatial movement of waves and particles is neglected, the coevolution of the electron distribution and waves in a plasma may be described in the general case by the system of equations (1)where is the energy density of oscillations of mode σ in the space of wave vectors, σ = 1,...,N, γ^{(σ)} is the growth rate for a given mode, k is the wave vector, f is the electron distribution function, is the diffusion tensor (describing electron scattering on the waves of mode σ), and p is the electron momentum. The equations in the system above are coupled implicitly, since the growth rate of plasma waves depends on the electron distribution function, while the diffusion tensor depends on the intensity and spectral distribution of the plasma waves. The expressions for the growth rate and elements of the diffusion tensor are given in Appendices B and C, respectively.
In this work, we explore the coevolution of the electron distribution and plasma waves using numerical simulations. The distributions f(p) and are defined on regular grids in (p,α) and (ω,θ)spaces, respectively, where α is the electron pitch angle, ω is the wave frequency, and θ is the wave propagation direction (with respect to the ambient magnetic field). For different modes, different grids are used. In most simulations, we considered two wave modes (e.g., ordinary + extraordinary). The grid size was chosen to be 60 × 60 data points for all considered distributions. In Aschwanden (1990), plasma waves were described using an irregular adaptive grid where the density of data points in phase space was approximately proportional to the initial growth rate. However, our simulations have demonstrated that for the ringlike or horseshoelike electron distributions, the region of positive growth in phase space can shift noticeably during the process of relaxation. Therefore, a regular grid is more suitable, and the considered area in (ω,θ)space has to be wider than the initial region of positive growth. The system of equations given in Eq. (1) is integrated with respect to time using the Gear formulae of fourth order (see Appendix D for details of the numerical code). All cyclotron harmonics affecting the waves growth/damping and particle diffusion are considered; however, as a rule, the first harmonic is the most dominant.
Since the simulated system is closed, both its energy and particle number must be conserved (the total system energy equals the sum of the energy of particles and the energies of all considered oscillation modes). Fulfillment of the conservation laws can be considered as a test of the selfconsistency in the model and the accuracy of the numerical code. In our simulations, at the late stage of the relaxation process, the particle number was conserved with the relative error ≲1.5 × 10^{3}, and the total energy of the system was conserved with the relative error ≲ 3 × 10^{3}. These estimates correspond to the models without the thermal plasma component; for models including the thermal plasma, the computation accuracy can be even higher.
3. Initial conditions
We assume that the initial electron distribution function has the form (2)where n is the total electron concentration, n_{b} is the concentration of nonthermal electrons, f_{0} is the Maxwellian distribution function of thermal electrons, and f_{b} is the distribution function of nonthermal electrons (both functions f_{0} and f_{b} are assumed to be normalized to unity). The cases of both n_{b} ≪ n (when the thermal component dominates) and n_{b} = n (when the thermal component is absent) are considered.
As stated above, we assume that the nonthermal electron distribution has a horseshoelike shape (the distributions observed in the terrestrial magnetosphere are shown, e.g., in Figs. 3 and 5 of Ergun et al. 2000). Instead of performing a detailed study of the formation of a horseshoelike distribution, we describe the initial electron distribution by the model function that is similar to the observed ones (3)where A is the normalization factor and μ = cosα. The shape of the distribution function is determined by such parameters as the typical electron momentum p_{b}, the electron dispersion in momentum Δp_{b}, the pitchangle boundary of the losscone α_{c} (or μ_{c} = cosα_{c}), and the losscone boundary width Δμ_{c}. An example of the distribution function given by Eq. (3) can be seen in Fig. 2a. For α_{c} = 0, we obtain an isotropic ringlike distribution.
The initial energy density of plasma waves is assumed to equal the level of thermal oscillations (4)where k_{B} is the Boltzmann constant and T_{0} is the effective plasma temperature. In addition, it is assumed that the energy density of plasma waves cannot fall below the thermal level during the process of wave/particle evolution, because of spontaneous radiation.
In all simulations, we assume that the nonthermal distribution function has Δp_{b}/p_{b} = 0.2 and Δμ_{c} = 0.2. The thermal component of the plasma (if present) is described by a Maxwellian distribution of temperature 10^{6} K. The initial temperature of plasma waves T_{0} equals 10^{6} K. The remaining parameters of the considered simulation models are given in Table 1; they were chosen to explore the influence of various factors on the process of wave/particle evolution. In all cases, the plasma density is relatively low, so that the ratio ω_{p}/ω_{B} varies in the range from 10^{3} to 10^{1}; the total electron concentration n is calculated using the plasma frequency ω_{p}. The relative concentration of the energetic electrons n_{b}/n varies from 10^{4} to 1.
Fig. 1 Initial growth rates of the extraordinary a) and ordinary b) waves. Simulation parameters correspond to the model 15 (Table 1), and the wave dispersion is assumed to be similar to that in a vacuum. 

Open with DEXTER 
Parameters of the different simulation models.
4. Results
4.1. Relaxation of the horseshoelike electron distribution
4.1.1. The case of when a thermal plasma component is absent
Firstly, we consider in detail an example of the coevolution of the electron distribution and plasma waves. We assume that a lowtemperature thermal component is absent (n_{b} = n). Investigations of the dispersion relations for weaklyrelativistic plasmas (Robinson 1986, 1987) have shown that the wave dispersion (for the waves propagating across the magnetic field,  cosθ  ≪ 1) becomes similar to that in a vacuum if the typical electron speed 3_{e} satisfies the condition (3_{e}/c) ≳ (ω_{p}/ω_{B})^{2}. This requirement is satisfied, e.g., for the particle energy E_{b} ≳ 3 keV and ω_{p}/ω_{B} ≲ 0.1. Thus we assume that the wave refraction index equals unity for both the ordinary and extraordinary modes. We also assume that the waves are elliptically polarized, the axial ratio of the polarization ellipse being T_{E} = cosθ for the extraordinary mode and T_{O} = −1/cosθ for the ordinary mode (T_{E}T_{O} = −1). The above relations follow from the magnetoionic theory when ω_{B}/ω → 1 and ω_{p}/ω → 0 (Melrose & Dulk 1991; Willes et al. 1994); however, we found that the simulation results are not very sensitive to the exact value of the axial ratio T_{σ} provided that for quasitransversal propagation  T_{E}  ≪ 1 (and, accordingly,  T_{O}  ≫ 1).
As an illustration, the following parameters were chosen: the magnetic field B = 1430 G that corresponds to the electron cyclotron frequency of f_{B} = 4 GHz (a typical value for the radio emission of ultracool dwarfs), the plasma to cyclotron frequency ratio ω_{p}/ω_{B} = 10^{3} that corresponds to the electron concentration n = n_{b} = 2 × 10^{5} cm^{3}, the typical energy of the energetic electrons E_{b} = 10 keV, and the losscone boundary α_{c} = 60°. In Table 1, this parameter set corresponds to the model 15.
Figure 1 shows the contours of growth rates of the ordinary and extraordinary modes at the initial moment (t = 0). Only the positive growth rates are shown; the contour levels are evenly distributed between zero and the maximal growth rate value for a given mode which is also shown at the figure. One can see that the most effective wave amplification takes place slightly below the fundamental cyclotron frequency. Both emission modes are generated mainly in the perpendicular direction with respect to the magnetic field, with a slight asymmetry due to the losscone feature. Growth rates decrease rapidly with increasing frequency, but wave amplification (in an oblique direction) can occur even above the cyclotron frequency; in this region, the emission directivity patterns become essentially asymmetric. One can also note that the maximal growth rate of the extraordinary mode exceeds that of the ordinary mode by more than two orders of magnitude. For both the ordinary and extraordinary modes, there are also amplification regions near higher harmonics of the cyclotron frequency, but with much lower growth rates (they are not shown in the figure). Thus, the fundamental extraordinary mode is the dominant one.
Fig. 2 Time evolution of the electron distribution for the model 15 (Table 1). 

Open with DEXTER 
Fig. 3 Time evolution of the extraordinary waves for the model 15 (Table 1). The contour levels are 0.0075, 0.015, 0.03, 0.06, 0.125, 0.25, 0.5, and 0.9 of the maximum wave energy density. 

Open with DEXTER 
Figure 2 shows the electron distribution function at different times (the chosen time moments correspond to the different stages of the relaxation process, see comment to Fig. 4). For some time after the beginning of the simulation, while the wave energy density is low, the distribution function remains very similar to the initial one (Fig. 2a). When the waves are amplified to a certain critical level, an effective diffusion of electrons on these waves begins. As a result, the electrons drift along the trajectories of p_{z} = const. towards lower values of p_{ ⊥ }, thus losing energy. At the early relaxation stage (Fig. 2b), this results in the formation of a plateau slightly below the maximum of the initial distribution function (with respect to the momentum) and around α ≃ 90° (with respect to the pitch angle). On the other hand, at the upper and lower boundaries of this plateau (with respect to the momentum), the gradient of the distribution function increases. At the middle relaxation stage (Fig. 2c), the “plateau” expands with respect to both the momentum and pitch angle, and we can see that this feature is not really flat: although the derivative ∂f/∂p_{ ⊥ } in this area approaches zero, the derivative ∂f/∂p_{z} remains nonzero, and therefore the distribution function gradually decreases with increasing p_{z}. Nevertheless, since the maser amplification or damping of waves depends mainly on the derivative with respect to the transversal component of the momentum, a contribution of this region of the distribution function to the growth rate is close to zero. At the upper boundary (with respect to the momentum), the “plateau” extends beyond the circle p = p_{b} (where the initial distribution function had a maximum) and the positive slope in p_{ ⊥ } disappears; at the lower boundary, this slope is still quite large. At the late relaxation stage (Fig. 2d), the “plateau” expands even further, so now we have ∂f/∂p_{ ⊥ } ≤ 0 almost everywhere, and amplification of waves nearly ceases. A positive slope remains only in the lowenergy region, but the resulting growth rate of the waves is very low. At even later times, no qualitatively new features arise; the simulations revealed that the “hole” around p = 0 slowly shrinks and the distribution function asymptotically approaches a stationary (saturated) state, while its rate of change decreases with time.
Figure 3 shows the distribution of the extraordinary mode energy density in (ω,θ)space at different times; in addition, the maximal values of the amplification coefficient Λ are shown (this coefficient equals the ratio of the energy density at a given time to its initial level). The letters identifying the panels correspond to those in Fig. 2, but panel (a) is omitted since the wave energy at the initial moment is negligible. For some time after the beginning of the simulation, the growth rate remains nearly constant, so the wave energy grows exponentially with time. As a result, at the onset of relaxation of the electron beam (Fig. 3b), the waves are concentrated in a relatively narrow range with respect to both the frequency and the propagation direction; this region is much narrower than the initial region of positive growth rate. When the logarithm of the amplification coefficient lnΛ reaches values of about 24–25, the waves begin to modify the electron distribution. This results in a sharp decrease (down to zero) in the growth rate in those regions where it was large at the beginning of the simulation, and also in an increase in the growth rate in the adjacent regions. These changes reflect the formation of a plateau on the electron distribution function (see Fig. 2b). As a result, the maximal intensity of waves nearly ceases to increase, but the wave distribution in phase space broadens (see Fig. 3c, which corresponds to the middle relaxation stage). At a later stage (Fig. 3d), the distribution of plasma waves becomes even broader. We note that the shape of the final wave distribution differs considerably from that of the initial region of positive growth; in particular, relaxation of the electron beam allows waves to be generated at noticeably lower frequencies than at the initial moment. On the other hand, the wave energy density in the frequency range above the cyclotron frequency remains negligible throughout the relaxation process, although the growth rate can be initially positive in this region.
Fig. 4 Time profiles of the total wave energy a) and average growth rate b). In panel a), the solid line is the result of numerical simulations and dashed line is a simplified functional fit. Simulation parameters correspond to the model 15 (Table 1). 

Open with DEXTER 
Figure 4 shows the time history of the integral characteristics of the extraordinary waves: the total energy density (integrated over frequency and propagation angle) and the average growth rate (only the positive values of the growth rate were considered). In general, the time history is very similar to that inferred by Aschwanden (1990). For some time after the beginning of the simulation (0 < t ≲ t_{ons}), the growth rate is constant and the wave energy grows exponentially but still remains very low. When the wave energy reaches a certain critical level (at t ≃ t_{ons}), relaxation of the unstable electron distribution begins and the growth rate starts to decrease. We note that this critical level is considerably lower than the final wave energy. After the onset of relaxation (at t_{ons} ≲ t ≤ t_{ss}), the total energy of waves continues to grow at an increasing rate, but the corresponding curve somewhat differs from an exponential one. At time t_{ss} (which is defined as the time of the steepest slope), the rate of change in the total wave energy ∂W/∂t reaches its maximum. Later (at t > t_{ss}), both the average growth rate and the rate of change in the total wave energy gradually decrease, approaching zero asymptotically; at the same time, the total wave energy asymptotically reaches the saturation level (W_{∞}). At this stage, the time profile of the total wave energy can be described with good accuracy by the wellknown saturation curve (5)where W_{ss} = W(t_{ss}). The distributions of the electrons and waves in Figs. 2, 3 correspond to the times t = 0, t = t_{ss} (early stage of relaxation), t = t_{ss} + τ_{sat} (middle stage), and t = t_{ss} + 3τ_{sat} (late stage). For the simulation parameters used in this section, relaxation begins at t_{ons} ≃ 6.5 μs or , where γ_{max} is the maximal growth rate of the extraordinary mode at the initial moment, the time of the steepest slope t_{ss} ≃ 9.45 μs or , and the relaxation timescale τ_{sat} ≃ 14.4 μs or . The final (at t → ∞) total energy density of the extraordinary waves equals W_{∞} ≃ 4.64 × 10^{4} erg cm^{3}, which amounts to 13.3% of the initial energy density of the accelerated particles. The derived timescales are similar (by an order of magnitude) to those given in Aschwanden (1990). However, the conversion efficiency of the particle energy into waves is now considerably higher because of the different type of unstable electron distribution (Aschwanden 1990, considered only the losscone).
As said before, the growth rate of the ordinary mode is much lower than that of the extraordinary mode. Simulations considering both emission modes simultaneously have shown that the ordinary mode is amplified by less than a factor of 1.25 relative to the level of thermal oscillations. Thus the resulting energy of the ordinary waves is negligible and they neither make a measurable contribution to the radio emission of planets and UCDs nor affect the electron distribution.
4.1.2. The case of when a thermal plasma component is present
We now investigate the case of when a thermal plasma component dominates. We use the same parameters of the magnetic field and energetic particles as in the previous section (f_{B} = 4 GHz, n_{b} = 2 × 10^{5} cm^{3}, E_{b} = 10 keV, α_{c} = 60°), but now a thermal plasma with the concentration of n_{0} = 2 × 10^{7} cm^{3} is present (n_{b}/n = 10^{2}), so that the plasma to cyclotron frequency ratio is ω_{p}/ω_{B} = 10^{2}. In Table 1, this parameter set corresponds to the model 9. For the wave parameters, we use the cold plasma dispersion relation (see Appendix A).
Fig. 5 Initial growth rates of the extraordinary a) and ordinary b) waves in the presence of a thermal electron component. Simulation parameters correspond to the model 9 (Table 1), and the cold plasma dispersion relation is used. 

Open with DEXTER 
Figure 5 shows the initial growth rates (at t = 0). In a cold magnetized plasma, the extraordinary mode is divided into two branches: the fast (X) mode exists above the cutoff frequency ω_{c}, while the slow (Z) mode exists below the resonance frequency ω_{r} (ω_{c} > ω_{r}). Under the conditions considered, both the cutoff and resonance frequencies almost coincide with the electron cyclotron frequency. However, the frequency gap between the amplification regions of the X and Zmodes is much larger than ω_{c} − ω_{r} since near the cyclotron frequency the waves are heavily damped by the thermal plasma component. The Zmode is generated across the magnetic field, while the Xmode is generated in an oblique direction (at θ ≃ 75°). The maximal growth rate of the Zmode exceeds that of the Xmode by more than an order of magnitude. The ordinary (O) mode has the lowest growth rate; its dispersion curve is continuous, but the amplification region is divided into two by the absorption band near the cyclotron frequency.
Simulations have shown that the time evolution of the electron distribution is very similar to that when the thermal plasma is absent (see Fig. 2). Similarly, the time evolution of the Zmode (which is the dominant mode) does not differ much from the evolution of the extraordinary mode shown in Fig. 3. Therefore, we do not present a complete time history of the waves and particles. Figure 6 shows the late relaxation stage. One can see that the “plateau” on the electron distribution function (formed by the nonthermal electrons) merges with the thermal electron population at this stage. The conversion efficiency of the particle energy into the Zmode waves (only the nonthermal component is considered) is about 11.8%, which is slightly less that for the model without the thermal plasma. The relaxation timescales are very similar to the values found in the previous section. We can conclude that in a lowdensity plasma (with ω_{p}/ω_{B} ≪ 1), the exact form of the dispersion relations is not very important, since both the vacuumlike and cold plasma dispersion relations provide almost the same results.
Fig. 6 Electron distribution function a) and the Zmode energy density b) at the late relaxation stage. The thermal component of the distribution function is shown only partially. Simulation parameters correspond to the model 9 (Table 1). 

Open with DEXTER 
Fig. 7 Time evolution of the Xmode waves for the model 9 (Table 1). The contour levels are 0.125, 0.25, 0.5, and 0.9 of the maximum wave energy density. 

Open with DEXTER 
We can see that the generation of the Zmode waves is very effective. However, these waves cannot escape from the generation region. Thus the question is of interest: can the intensities of the freely propagating modes (X and/or Omode) reach sufficiently high levels despite the lower growth rate? We simulated the simultaneous evolution of the different modes. Figure 7 shows the distribution of the Xmode energy density at different times. The time history is different from that of the Zmode since the Xmode generation is caused by the losscone feature. Since the losscone is destroyed during relaxation (because of diffusion of particles on the Zmode waves), the Xmode amplification in some regions of phase space is replaced by absorption, so that the region occupied by the waves shrinks with time. Simulations for times later than those displayed in the figure have shown that the maximal intensity of waves can also somewhat decrease. However, the main result for the Xmode is that its intensity is much lower than that of the Zmode: the parameter lnΛ does not exceed four, so the wave energy density exceeds the thermal level by not more than a factor of 50. Therefore, the Xmode has no noticeable effect on the relaxation of the electron distribution.
Figure 8 shows the average growth rates of the Z and Xmodes. One can see that the growth rate of the Xmode starts to decrease simultaneously with that of the Zmode. In general, the ratio of growth rates remains nearly the same throughout the relaxation process (despite a slight secondary increase in the average growth rate of the Xmode at the middle relaxation stage), so that the growth rate of the Xmode remains considerably lower than that of the Zmode. We found that the total energy density of the Xmode waves does not exceed 10^{13} erg cm^{3}, which is equivalent to 3 × 10^{11} of the energy density of the accelerated particles. The energy density of the Omode is much lower than that of the Xmode. Thus, in this simulation, the relaxation of the unstable electron distribution is caused entirely by the Zmode, and the intensities of the other wave modes are negligible.
Fig. 8 Time profiles of the average growth rates of the Z and Xmodes for the model 9 (Table 1). For the Xmode, the scale differs from that for the Zmode to enhance the growth rate variations. 

Open with DEXTER 
4.2. Comparison of the results for the different simulation models
The simulations of the coevolution of the electron distributions and plasma waves were performed for the various parameter sets (see Table 1). The cold plasma and vacuum dispersion relations were used for the models with n_{b}/n ≪ 1 and n_{b}/n = 1, respectively. In all cases, the extraordinary mode (or its slow branch, for the models including a thermal plasma) dominated strongly. Table 1 also contains the maximal initial growth rates, the timescales of the relaxation process, and the transformation coefficient of the particle energy into the waves. Figures 9–13 illustrate the effect of the various factors.
Fig. 9 Energy conversion efficiency and the characteristic timescales of the relaxation process for the different concentrations of the energetic electrons. The values correspond to the models 5, 9, and 19 (Table 1). 

Open with DEXTER 
4.2.1. Effect of varying the plasma parameters
Figure 9 shows the simulation results (relaxation timescales and conversion efficiency of the particle energy into the waves) for the different relative concentrations of the energetic electrons (n_{b}/n = 10^{4}, 10^{2}, and 1), while the total plasma density is assumed to be constant (n = 2 × 10^{7} cm^{3}, which corresponds to ω_{p}/ω_{B} = 10^{2}). The typical energy of the accelerated electrons E_{b} and the losscone boundary α_{c} equal 10 keV and 60°, respectively. One can see that the relaxation timescales are simply inversely proportional to the concentration of the energetic particles. The conversion efficiency of the particle energy into the waves increases slightly with increasing n_{b}/n varying from 11.2% to 13.6%.
Fig. 10 Same as in Fig. 9, for the different values of the total plasma density (concentration of the accelerated electrons n_{b} is constant). The values correspond to the models 15, 9, and 6 (Table 1). 

Open with DEXTER 
Figure 10 shows the simulation results for the case of when the concentration of the energetic particles is constant (n_{b} = 2 × 10^{5} cm^{3}), while the total plasma density varies (which results in different ratios of both n_{b}/n and ω_{p}/ω_{B}). With increasing ω_{p}/ω_{B} from 10^{3} to 10^{2}, both the growth rate of the extraordinary mode and the relaxation timescales remain nearly unchanged. At the same time, with increasing ω_{p}/ω_{B} from 10^{2} to 10^{1}, the growth rate decreases by about an order of magnitude (because of the changing dispersion parameters of the waves), while the relaxation timescales increase by the same factor. The conversion efficiency of particle energy into waves varies around 12% and slightly decreases with increasing ω_{p}/ω_{B}.
4.2.2. Effect of varying the distribution of the energetic electrons
We now consider different distributions of the energetic electrons. In all cases, the concentration of the energetic electrons is assumed to be the same (n_{b} = 2 × 10^{5} cm^{3}). Both the models without a thermal plasma component and the models including such a component (with n_{b}/n = 10^{2}) are considered; the corresponding points at Figs. 11, 12 are connected by the solid and dashed lines, respectively.
Fig. 11 Same as in Fig. 9, for the different energies of the accelerated electrons. Solid lines: models 13, 15, and 18 (without thermal component); dashed lines: models 7, 9, and 12 (with thermal component). 

Open with DEXTER 
Figure 11 shows the simulation results for the different typical energies of the accelerated particles (in all cases, the losscone boundary is α_{c} = 60°). For the models without a thermal plasma component, the conversion efficiency of the particle energy into the waves is almost constant (about 13.3%). For the models including a thermal plasma, the conversion efficiency increases with E_{b} (from 8% at E_{b} = 3 keV to 13% at E_{b} = 30 keV). An increase in the electron energy makes the relaxation process slower (both the timescales t_{ss} and τ_{sat} increase). This is because the “plateau” formation for the electron distributions with higher energy requires a larger displacement of particles in the momentum space and, consequently, a higher energy density of plasma waves (to ensure a stronger diffusion) and a longer time. For the models with a thermal component, the relaxation process is slower. With increasing electron energy, the relaxation timescales (as well as the energy conversion efficiency) for the models with a thermal component approach the corresponding values for the models without a thermal component.
Fig. 12 Same as in Fig. 9, for the different angular distributions of the accelerated electrons. Solid lines: models 14, 15, 16, and 17 (without thermal component); dashed lines: models 8, 9, 10, and 11 (with thermal component). 

Open with DEXTER 
Figure 12 shows the simulation results for the different losscone boundaries (in all cases, the beam energy is E_{b} = 10 keV). The models without a thermal plasma component always provide a higher conversion efficiency of the particle energy into the waves than the corresponding models with a thermal component. The highest conversion efficiency as well as the fastest relaxation occur for α_{c} = 60°. For the ringlike distribution (with α_{c} = 0), as well as for the losscone with α_{c} = 90°, the conversion efficiency and relaxation rate are slightly lower. The distribution with α_{c} = 120° is similar to those used in the paper of Bingham & Cairns (2000) and Bingham et al. (2001); it has no electrons with pitch angles around 90°. Therefore, for this distribution, the amount of free energy is relatively low and only about 5–6% of the initial electron energy can be transferred to waves. We highlight that for all considered distributions (including essentially anisotropic ones), the waves are generated preferably in the perpendicular direction to the magnetic field; for instance, for the distribution with α_{c} = 120°, the maximum of the wave intensity (at the late relaxation stage) is at θ ≃ 91°.
Fig. 13 Same as in Fig. 9, for the different strengths of the magnetic field. Solid lines: models 2, 4, and 15 (without thermal component); dashed lines: models 1, 3, and 9 (with thermal component). 

Open with DEXTER 
4.2.3. Effect of varying the magnetic field
We now investigate the influence of the magnetic field strength on the coevolution of the electron distributions and plasma waves. Figure 13 shows the simulation results for three values of the electron cyclotron frequency: 400 kHz (which is typical of the AKR sources in the terrestrial magnetosphere), 40 MHz (which is typical of the magnetosphere of Jupiter), and 4 GHz (which seems to be typical of the magnetospheres of UCDs). We assume that the plasma to cyclotron frequency ratio is constant and equals ω_{p}/ω_{B} = 10^{3} for the models without a thermal plasma component; thus the concentrations of the energetic electrons n_{b} equal 2 × 10^{3}, 2 × 10^{1}, and 2 × 10^{5} cm^{3}, respectively. In the models including a thermal plasma, we assume that ω_{p}/ω_{B} = 10^{2} and n_{b}/n = 10^{2}, which provides the same concentrations of the energetic electrons as above. The typical energy of the accelerated electrons E_{b} and the losscone boundary α_{c} equal 10 keV and 60°, respectively. Under the conditions considered, the growth rates of the plasma waves are proportional to the cyclotron frequency (see Table 1). On the other hand, the relaxation timescales decrease with increasing cyclotron frequency somewhat faster than the inverse proportionality law implies. The longest relaxation timescales (relative to the inverse cyclotron frequency) occur at the Earth since in a weaker magnetic field diffusion of particles on plasma waves also weakens, and the relaxation onset requires a higher wave amplification coefficient Λ (this effect is equivalent to that of reducing the initial temperature of plasma waves). The fastest relaxation (both in absolute and relative units) occurs for UCDs. A complete relaxation of the unstable electron distribution is achieved in less than 1 s at the Earth, less than 10 ms at Jupiter, and less than 0.1 ms at UCDs. The conversion efficiency of the particle energy into the waves slightly increases with increasing magnetic field strength; the models without a thermal plasma component always provide a higher conversion efficiency (13.1–13.3%) than the corresponding models with a thermal component (11.8–12.0%).
5. Discussion
Our simulations have shown that in a relatively lowdensity plasma with ω_{p}/ω_{B} = 10^{3} − 10^{1}, the ringlike or horseshoelike distributions of accelerated electrons mainly generate the extraordinary waves with a frequency slightly below the electron cyclotron frequency. This conclusion is consistent with the results of previous studies (e.g., Ergun et al. 2000). If the magnetoionic theory is applicable (i.e., the wave dispersion is dominated by a cold plasma component), the generated waves correspond to the Zmode branch of the extraordinary mode. Generation of other modes (such as the ordinary mode, the Xmode of the magnetoionic theory, and the waves near higher cyclotron harmonics) is also possible, but their final energy density is less than that of the fundamental extraordinary mode by many orders of magnitude, even if the difference in the initial growth rates is not so large. We have found that the dominant mode remains the same throughout the relaxation process of an unstable electron distribution. This differs from the results of Fleishman & Arzner (2000), where relaxation of a losscone distribution on the initially dominant mode (lowerhybrid waves) resulted in formation of a distribution that was stable with respect to excitation of lowerhybrid waves but still sufficiently anisotropic to amplify other types of waves (X and O). This is, most probably, because Fleishman & Arzner (2000) considered a different electron distribution (losscone vs. horseshoe) and different plasma parameters (ω_{p}/ω_{B} > 0.2). In contrast, in our simulations, the relaxation of an electron distribution on the waves of the dominating mode prevented the amplification of all other modes. We also note that for other types of unstable electron distributions, another wave mode can prevail; for example, a losscone distribution in a lowdensity plasma excites mainly the oblique Xmode (Aschwanden 1990).
For the ringlike or horseshoelike electron distributions, the transformation coefficient of the particle energy into waves can exceed 10%. This is considerably higher than for the losscone distributions (Aschwanden 1990). The reason is that for the losscone, only a small fraction of particles (with pitchangles around α_{c}) actually contributes to the generation of waves and provides them with energy; in contrast, for the horseshoelike distribution, the process of wave generation involves the majority of particles. We note that the conversion efficiency is very weakly dependent on the plasma parameters and is determined mainly by the shape of the nonthermal distribution. The presence of a thermal plasma component reduces the conversion efficiency.
We now estimate the required efficiency of the emission mechanism. As an example, we consider the M 9 dwarf TVLM 513–46546 located at a distance of d = 10.5 pc. This object is known to produce radio bursts with an intensity up to I = 6 mJy (at the frequency about 8 GHz) (Hallinan et al. 2007). The typical size of the AKR generation region (in the direction across the magnetic field) can be estimated to be 300 km. The radius of the UCD is expected to be comparable to that of Jupiter (i.e., ten times larger than that of the Earth), so we can assume that the radiation source has a size of about R_{ ⊥ } ≃ 3000 km. Thus the brightness temperature of the emission is about 10^{13} K. If the initial temperature of plasma waves equals T_{0} = 10^{6} K, then the logarithm of the amplification coefficient should be not less than lnΛ ≃ 16. For the growth rate γ ≃ 3 × 10^{6} s^{1} (see Sect. 4.1.1 for the beam and plasma parameters) and group velocity of the waves 3_{gr} ≤ c, this amplification can be achieved at a distance of no more than 1.6 km. Thus the estimated growth rates of the ECMI are well in excess of the required values, even for the relatively low concentrations of the accelerated electrons.
On the other hand, as stated above, estimates based only on the growth rate are insufficient and we need to check the conversion efficiency of the electron energy into radiation. The total power of radio emission from the UCD (for the case of isotropic radiation) can be estimated as F_{r} = 4πd^{2}IΔf, where Δf is the spectral bandwidth of the emission. For the above parameters, assuming that Δf ≃ f ≃ 8 GHz (this is an upper limit as the spectrum is unlikely to be this broad), we obtain F_{r} ≃ 6 × 10^{21} erg s^{1}. The energy flux of the accelerated electrons can be estimated to be , where n_{b}, 3_{b}, and E_{b} are the concentration, speed, and energy of particles, and is the crosssection area of the electron beam. We assume that the plasma to cyclotron frequency ratio is ω_{p}/ω_{B} = 10^{3} (which is smaller than that in the AKR sources) and all particles in the radiation source are accelerated (n_{b} = n). For the cyclotron frequencies f_{B} = 4 − 8 GHz, we then find that n_{b} ≃ (2 − 8) × 10^{5} cm^{3}. For the electron energy E_{b} = 10 keV, the corresponding speed 3_{b} = 5 × 10^{9} cm s^{1}, and the beam size R_{ ⊥ } = 3000 km, the energy flux of the particles is F_{b} ≃ (1.4 − 5.7) × 10^{24} erg s^{1}. Thus the conversion efficiency of the electron energy into the radio emission should be not less than F_{r}/F_{b} ≃ (1 − 4) × 10^{3}. According to our simulations, the ECMI efficiency ( ≳ 0.1) is much higher. The above estimates do not take into account the possible absorption of the radio emission during propagation, the radiation directivity, and (possibly) a narrow spectral band. Nevertheless, with reasonable assumptions about the source parameters, the ECMI is capable of providing the observed intensity of radio emission from UCDs.
The terrestrial AKR is generated in the auroral cavity where a cold plasma component is almost absent; this allows the waves produced below the electron cyclotron frequency to escape directly from the source region due to relativistic corrections to the dispersion relation. It is very likely that the radio emission of UCDs is generated under similar conditions. However, we discuss a possible case of when a cold plasma component dominates. Under these conditions, the waves (of the Zmode) excited by the ECMI cannot escape from the source because of a stop band at ω ≃ ω_{B}. Thus the Zmode has to be transformed into electromagnetic radiation, e.g., due to nonlinear processes. We can estimate the required efficiency of nonlinear conversion to be η_{NL} ≳ (1 − 4) × 10^{2} (when ω_{p}/ω_{B} = 10^{2}, n_{b}/n = 10^{2}, 10% of the particle energy goes into the Zmode waves, and all other parameters are the same as in the previous paragraph). Nonlinear processes require further investigation, but at present we cannot rule out the described twostage emission mechanism. We note that the growth rates of the X and Omodes (see Fig. 5), in a linear approximation, can also provide the required amplification coefficient at distances of from tens to hundreds of kilometers. However, the kinetic simulations have shown that the conversion coefficient of the particle energy into these waves is extremely small, so a direct generation of electromagnetic waves by the horseshoelike electron distribution in the presence of a cold plasma cannot be responsible for the observed emission.
According to the simulation results (see Table 1), the typical relaxation time of electron distributions with ω_{p}/ω_{B} = 10^{2} and n_{b}/n ≃ 1 in the terrestrial magnetosphere should be about τ_{diff} ≃ τ_{sat} ≃ 2 × 10^{3} s. On the other hand, the typical escape time of the waves from the generation region is τ_{esc} ≃ R_{ ⊥ }/3_{gr}. For R_{ ⊥ } ≃ 300 km and 3_{gr} ≃ c, we obtain τ_{esc} ≃ 10^{3} s. Thus τ_{diff} > τ_{esc}, and the diffusive limit is, in fact, not applicable. In the terrestrial magnetosphere, the escape of waves from the amplification region has a dominating effect on the relaxation of the electron beams. In particular, it allows the strongly unstable electron distributions with large values of ∂f/∂p_{ ⊥ } (such as those reported by Ergun et al. 2000) to exist in a quasistationary state. The escape of waves from the amplification region obviously reduces the relaxation rate, so that the processes of particle acceleration and magnetic mirroring are able to compensate for the changes in the distribution function caused by the diffusion on plasma waves. On the other hand, for UCDs we find that τ_{esc} ≃ 10^{2} s and τ_{diff} ≃ 10^{7} − 10^{5} s; thus τ_{diff} ≪ τ_{esc}, and the use of a diffusive limit for UCDs is clearly justified. We can expect the quasistationary electron distributions in the magnetospheres of UCDs to be very close to the relaxed state (as shown, e.g., at Figs. 2c,d). To derive more accurate results, one has to include the mechanism responsible for formation of an unstable electron distribution into the simulation model. The model that includes the escape of waves from the amplification region will be considered in a future paper.
The relaxation of unstable electron distributions in the magnetospheres of UCDs is very fast ( ≪ 1 s). On the other hand, the observed bursts of radio emission have a much longer duration (tens of seconds). It is very likely that the emission source is actually even more longlived, but we observe the bursts only when the narrowlydirected radio beam (whose direction changes owing to the rotation of the star) points towards the observer (Hallinan et al. 2006, 2007; Berger et al. 2009). This again means that the diffusive limit can be considered only as a very rough approximation. Actually, the emission properties are determined by the relatively stationary (but longliving) processes of particle acceleration; any model attempting to quantitatively interpret the observations has to include these processes. Nevertheless, we believe that the results obtained in this paper allow us to draw some conclusions about the processes in the magnetospheres of UCDs, and will be useful in developing more advanced models.
6. Conclusion
We have performed a numerical simulation of the electroncyclotron maser instability and investigated the coevolution of an unstable electron distribution and plasma waves in a lowdensity plasma (ω_{p}/ω_{B} ≪ 1). The electron distribution of the “horseshoe” type was considered as an energy source for the plasma oscillations. The simulation parameters were chosen according to the measurements in the terrestrial magnetosphere, and scaled to the conditions expected in the magnetospheres of UCDs (where the electron cyclotron frequency is a few GHz). The spatial movement of the waves and particles, as well as other processes modifying the electron distribution (except for the diffusion on plasma waves), were neglected. A number of simulations were performed for the different parameters of plasma and accelerated particles. We found that:

Under the conditions considered, the electron beams withringlike or horseshoelike distributions mainly generate theextraordinary waves that propagate across the magnetic field,with a frequency slightly below the electron cyclotron frequency.Other wave modes can also be amplified but their final intensity isless than that of the fundamental extraordinary mode by severalorders of magnitude.

The conversion efficiency of the energy of accelerated electrons into the extraordinary waves can exceed 10%, the parameter depending mainly on the shape of the nonthermal distribution.

A typical relaxation time of the unstable electron distribution is a few tens of inverse growth rates. In the magnetospheres of UCDs, a complete relaxation of the electron distribution should occur very rapidly (in much less time than one second); thus the observed emission requires a longlived source of accelerated particles.

Energy estimates show that the efficiency of the electroncyclotron maser instability is sufficient to provide the observed intensity of radio emission from UCDs under reasonable assumptions about the energies and concentrations of accelerated particles.
Acknowledgments
The Armagh Observatory is supported by a grant from the Northern Ireland Dept. of Culture Arts and Leisure. We also thank the Leverhulme Trust for financial support, without whose funding this work would not have been possible.
References
 Antonova, A., Doyle, J. G., Hallinan, G., Bourke, S., & Golden, A. 2008, A&A, 487, 317 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Aschwanden, M. J. 1990, A&AS, 85, 1141 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] (In the text)
 Aschwanden, M. J., & Benz, A. O. 1988, ApJ, 332, 447 [NASA ADS] [CrossRef] (In the text)
 Berger, E. 2005, ApJ, 648, 629 [NASA ADS] [CrossRef] (In the text)
 Berger, E., Rutledge, R. E., PhanBao, N., et al. 2009, ApJ, 695, 310 [NASA ADS] [CrossRef] (In the text)
 Bingham, R., & Cairns, R. A. 2000, Phys. Plasmas, 7, 3089 [NASA ADS] [CrossRef] (In the text)
 Bingham, R., Cairns, R. A., & Kellett, B. J. 2001, A&A, 370, 1000 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Delory, G. T., Ergun, R. E., Carlson, C. W., et al. 1998, Geophys. Res. Lett., 25, 2069 [NASA ADS] [CrossRef] (In the text)
 Ergun, R. E., Carlson, C. W., McFadden, J. P., et al. 2000, ApJ, 538, 456 [NASA ADS] [CrossRef] (In the text)
 Donati, J.F., Forveille, T., Cameron, A. C., et al. 2006, Science, 311, 633 [NASA ADS] [CrossRef] [PubMed] (In the text)
 Fleishman, G. D., & Melnikov, V. F. 1998, Physics–Uspekhi, 41, 1157 [NASA ADS] [CrossRef] (In the text)
 Fleishman, G., & Arzner, K. 2000, A&A, 358, 776 [NASA ADS] (In the text)
 Hallinan, G., Antonova, A., Doyle, J. G., et al. 2006, ApJ, 653, 690 [NASA ADS] [CrossRef] (In the text)
 Hallinan, G., Bourke, S., Lane, C., et al. 2007, ApJ, 663, 25 [NASA ADS] [CrossRef] (In the text)
 Hallinan, G., Antonova, A., Doyle, J. G., et al. 2008, ApJ, 684, 644 [NASA ADS] [CrossRef] (In the text)
 Le Quéau, D., & Louarn, P. 1989, J. Geophys. Res., 94, A2605 [NASA ADS] [CrossRef] (In the text)
 Louarn, P., & Le Quéau, D. 1996, Planet. Space Sci., 44, 211 [NASA ADS] [CrossRef] (In the text)
 Melrose, D. B., & Dulk, G. A. 1982, ApJ, 259, 844 [NASA ADS] [CrossRef] (In the text)
 Melrose, D. B., & Dulk, G. A. 1991, A&A, 249, 250 [NASA ADS] (In the text)
 Morin, J., Donati, J. F., Petit, P., et al. 2008, MNRAS, 390, 567 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
 Mutel, R. L., Menietti, J. D., Christopher, I. W., Gurnett, D. A., & Cook, J. M. 2006, J. Geophys. Res., 111, A10203 [NASA ADS] [CrossRef] (In the text)
 Mutel, R. L., Peterson, W. M., Jaeger, T. R., & Scudder, J. D. 2007, J. Geophys. Res., 112, A07211 [NASA ADS] [CrossRef] (In the text)
 Pottelette, R., Treumann, R. A., & Berthomier, M. 2001, J. Geophys. Res., 106, A8465 [NASA ADS] [CrossRef] (In the text)
 Reiners, A., & Basri, G. 2007, ApJ, 656, 1121 [NASA ADS] [CrossRef] (In the text)
 Robinson, P. A. 1986, J. Plasma Phys., 35, 187 [NASA ADS] [CrossRef] (In the text)
 Robinson, P. A. 1987, J. Plasma Phys., 37, 149 [NASA ADS] [CrossRef] (In the text)
 Strangeway, R. J. 1985, J. Geophys. Res., 90, A9675 [NASA ADS] [CrossRef] (In the text)
 Strangeway, R. J. 1986, J. Geophys. Res., 91, A3152 [NASA ADS] [CrossRef] (In the text)
 Treumann, R. A. 2006, A&AR, 13, 229 [NASA ADS] [CrossRef] (In the text)
 Willes, A. J., Melrose, D. B., & Robinson, P. A. 1994, J. Geophys. Res., 99, A21203 [NASA ADS] [CrossRef] (In the text)
 Winglee, R. M. 1983, J. Plasma Phys., 25, 217 [NASA ADS] [CrossRef] (In the text)
 Winglee, R. M. 1985, ApJ, 291, 160 [NASA ADS] [CrossRef] (In the text)
 Winglee, R. M., & Dulk, G. A. 1986, ApJ, 307, 808 [NASA ADS] [CrossRef] (In the text)
 Wu, C. S., & Lee, L. C. 1979, ApJ, 230, 621 [NASA ADS] [CrossRef] (In the text)
 Zarka, P. 1998, J. Geophys. Res., 103, E20159 [NASA ADS] [CrossRef] (In the text)
Appendix A: Dispersion of the magnetoionic modes
The refraction index of the electromagnetic waves (N_{σ}) in a cold magnetized plasma satisfies the dispersion equation (A.1)where (A.2)(A.3)ω and k are the wave frequency and wave vector, ω_{p} and ω_{B} are the electron plasma and cyclotron frequencies respectively, and θ is the angle between the wave vector and the magnetic field. In Eq. (A.1), σ = −1 for the Xmode (fast extraordinary) and the Zmode (slow extraordinary); σ = + 1 for the Omode (fast ordinary) and the Wmode (whistlers).
The polarization state of the magnetoionic modes is described by the parameters (A.4)(A.5)where T_{σ} is the axial ratio of the polarization ellipse and L_{σ} is the longitudinal part of the polarization.
The group velocity of the magnetoionic waves equals (A.6)where (A.7)The derivative of the refraction index with respect to the propagation direction can be calculated using the relation (A.8)The different magnetoionic modes exist in the following frequency ranges:

Xmode: at ω > ω_{c + };

Omode: at ω > ω_{p};

Zmode: at ω_{c − } < ω < ω_{r + };

whistlers: at ω < ω_{r − }.
The cutoff and resonance frequencies are given by (A.9)(A.10)In a vacuum, the dispersion parameters become (A.11)By expanding the polarization parameters (A.4, A.5) in the small parameters and and retaining only the zeroorder terms in the expansions, we obtain the polarization state of the magnetoionic modes near the cyclotron frequency in a lowdensity plasma (A.12)
Appendix B: Expressions for the growth rate
According to Melrose & Dulk (1982) and Aschwanden (1990), the growth rate of the magnetoionic oscillations equals (the wavemode index σ is hereafter omitted for brevity) (B.1)where v and p are the electron velocity and momentum, the indices z and ⊥ indicate the longitudinal and transversal components of the vectors with respect to the magnetic field (in particular, N_{z} = Ncosθ and N_{ ⊥ } = Nsinθ), β = 3/c, Γ = (1 − β^{2})^{ − 1/2} is the relativistic factor, J_{s}(λ) and are the Bessel function and its derivative over the argument λ, (B.2)and f(p) is the electron distribution function satisfying the normalization condition (B.3)where n is the total electron concentration. The growth rates of the vacuum modes can be obtained by substituting the corresponding dispersion parameters given by Eqs. (A.11, A.12).
If we introduce the dimensionless parameters (B.4)and, in addition, use the polar coordinates (u,α) for the distribution function (where α is the angle between the electron velocity and the magnetic field), then the expression for the growth rate takes the form (B.5)The dimensionless distribution function f(u) should satisfy the normalization condition (B.6)and λ = xN_{ ⊥ }u_{ ⊥ }.
Using the properties of the δfunction, we can reduce threedimensional integrals in Eq. (B.5) to onedimensional integrals over du_{z}(B.7)where the value u_{ ⊥ } at every point of the resonance curve is found from the resonance condition (B.8)and the integration limits u_{zmin} and u_{zmax} (which are dependent on the harmonic number) are the boundaries of the interval where Eq. (B.8) has a solution.
For  N_{z}  < 1 (which is always satisfied for X and Omodes), the resonance curve in the momentum space in Eq. (B.8) is an ellipse intersecting the axis u_{ ⊥ } = 0 at two points (B.9)where u_{z1} and u_{z2} correspond to signs “ − ” and “ + ”, respectively. In this case, u_{zmin} = u_{z1} and u_{zmax} = u_{z2}. For  N_{z}  > 1, Eq. (B.8) describes a hyperbola with only one branch being physical (for which, at u → ∞, the longitudinal component of the momentum u_{z} has the same sign as N_{z}). Therefore the integration limits in Eq. (B.7) equal (note that u_{z1} > u_{z2} in this case): u_{zmin} = −∞ and u_{zmax} = u_{z2} at N_{z} < 0; u_{zmin} = u_{z1} and u_{zmax} = ∞ at N_{z} > 0. Infinite integration limits can be avoided since the electron distribution function is usually defined only in a finite range of momenta, e.g., at u ≤ u_{high}. In this case, the infinite limit (upper or lower, depending on the sign of N_{z}) should be replaced by the coordinate of the intersection point of the resonance curve with the circle u = u_{high}, that is (B.10)where Γ_{high} is the relativistic factor of electrons with the momentum u_{high}.
The above formulae for the growth rate involve infinite sums over cyclotron harmonics s. However, only the harmonics in a certain range s_{min} ≤ s ≤ s_{max} contribute to the growth rate, since for the other harmonics either the resonance condition (for a given wave parameters) is never satisfied or the resonance curve lies outside the domain of the electron distribution function. In this work, the interval of harmonic numbers is taken with a large excess (say, from − 100 to 100), and then for each harmonic we verify whether it contributes to the growth rate; this method is simple and very fast, and allows us to take into account all harmonics having an effect on the growth rate.
Appendix C: Expressions for the quasilinear diffusion rate
The change in the electron distribution function due to diffusion on the magnetoionic waves is described by the last equation of the system of equations given in Eq. (1). In polar coordinates (p,α), this equation takes the form (Aschwanden & Benz 1988) (C.1)where D_{pp}, D_{pα}, D_{αp}, and D_{αα} are the components of the diffusion tensor (D_{pα} = D_{αp}). We note that the total diffusion tensor equals the sum of the diffusion tensors on separate modes (for simplicity, the formulae below refer only to one mode).
If we use the dimensionless momentum u and introduce, according to Aschwanden & Benz (1988), generalized diffusion coefficients (C.2)then Eq. (C.1) can be written in the form (C.3)where the expression elements are rearranged to make the numerical calculations more convenient and accurate.
The diffusion coefficients D_{r} (for an individual wave mode) are (Aschwanden & Benz 1988; Aschwanden 1990) (C.4)where W_{k} is the energy density of the considered mode in the space of wave vectors. The total energy density (the energy per unit volume) equals (C.5)The diffusion coefficients for the vacuum modes can be obtained by substituting the corresponding dispersion parameters (A.11, A.12).
If we use in Eq. (C.4) the dimensionless frequency x and the propagation direction θ as the wave characteristics instead of the wave vector, then the diffusion coefficients become (C.6)where is the relative (with respect to the level of thermal oscillations) energy density of plasma waves in the space of wave vectors (C.7)According to the initial conditions used, at the initial moment we have .
The integral over dθ in (C.6) can be calculated analytically by using the properties of the δfunction (see Aschwanden & Benz 1988; Aschwanden 1990). As a result, the expression for the diffusion coefficients takes the form (C.8)where the integration limits x_{min} and x_{max} correspond to the domain of the function .
The resonance curve (dependence of θ on x) is found using the following considerations. If the particle momentum u and wave frequency x are known, then the resonance condition results in (C.9)If the wave dispersion is assumed to be as in a vacuum, then N = 1 and we immediately obtain cosθ = N_{z} (the requirement  cosθ  ≤ 1 must be satisfied). If the cold plasma dispersion relation is used, then the solution becomes a bit more complicated. With a known N_{z} = Ncosθ, the dispersion Eq. (A.1) can be transformed to the biquadratic equation in the variable η = cosθ(C.10)where (C.11)(C.12)(C.13)This equation has a solution (C.14)To determine the correct sign (“ + ” or “ − ”) in the above formula for a given mode, one has to check whether the obtained values of the angle θ satisfy the dispersion Eq. (A.1). In practice, it is sufficient to check the validity of the inequality (C.15)In addition, the obvious requirements and η^{2} ≤ 1 must be satisfied. Depending on the conditions, Eq. (C.10) for a given magnetoionic mode can have up to two solutions. If a solution with respect to η^{2} exists then the sign of η coincides with that of N_{z} in Eq. (C.9).
The interval of harmonic numbers contributing to the diffusion coefficients (s_{min} ≤ s ≤ s_{max}) can be estimated analytically, by using the resonance condition. As a result, for X and Omodes (which have the refraction index N ≤ 1) we obtain (C.16)For the Zmode and whistlers, the refraction index is not limited from above (in the cold plasma approximation), and we have to introduce these limits artificially by assuming that N ≤ N_{max}. Then we obtain (C.17)In this work, we use N_{max} = 10, which exceeds the actual values for the waves in the region of positive growth rate. The above formulae, as a rule, estimate the range of harmonic numbers with a large excess.
Appendix D: Numeric code
In our simulations, the electron distribution function and energy density of plasma waves are represented as twodimensional arrays (D.1)where 0 ≤ i ≤ N_{x} − 1, 0 ≤ j ≤ N_{θ} − 1, 0 ≤ m ≤ N_{u} − 1, and 0 ≤ n ≤ N_{α} − 1. As a result, the growth rate and diffusion coefficients need to be calculated only in the nodes of the corresponding grids, and these values (for each considered mode) can also be represented as twodimensional arrays (D.2)Since the growth rate and diffusion coefficients are linearly dependent on the electron distribution function and wave energy density, respectively (see Eqs. (B.1) and (C.4)), they can be computed using tensor multiplication (Fleishman & Arzner 2000) (D.3)where the kernels R_{ijmn} and are constant throughout the simulation process. They can therefore be computed once, which reduces the computation time considerably.
Since the resonance curves in the spaces of particle momentums and wave vectors, as a rule, do not pass through the grid points, the calculation of the growth rate and diffusion coefficients in the discretized case requires interpolating. In this work, we use a bilinear interpolation for the energy density of plasma waves. For the electron distribution function, we need to calculate not the function itself but its partial derivatives; this is done using a linear interpolation in one variable and a Lagrange interpolation on three closest points in another variable (where the derivative is needed).
For the electron distribution function, we use the grid where the momentum nodes are evenly distributed in the open interval from u_{low} to u_{high}, and the pitchangle nodes are evenly distributed in the open interval from 0 to π, so that (D.4)For this grid, the best results are achieved (i.e., the conservation laws for the total energy and particle number are fulfilled with the highest accuracy) with the following boundary conditions which are numerically implemented as an extrapolation of the arrays f_{mn} and beyond the grid (D.5)For the energy density of plasma waves , boundary conditions are not needed since the equations used do not contain derivatives of this value with respect to the wave parameters.
All Tables
All Figures
Fig. 1 Initial growth rates of the extraordinary a) and ordinary b) waves. Simulation parameters correspond to the model 15 (Table 1), and the wave dispersion is assumed to be similar to that in a vacuum. 

Open with DEXTER  
In the text 
Fig. 2 Time evolution of the electron distribution for the model 15 (Table 1). 

Open with DEXTER  
In the text 
Fig. 3 Time evolution of the extraordinary waves for the model 15 (Table 1). The contour levels are 0.0075, 0.015, 0.03, 0.06, 0.125, 0.25, 0.5, and 0.9 of the maximum wave energy density. 

Open with DEXTER  
In the text 
Fig. 4 Time profiles of the total wave energy a) and average growth rate b). In panel a), the solid line is the result of numerical simulations and dashed line is a simplified functional fit. Simulation parameters correspond to the model 15 (Table 1). 

Open with DEXTER  
In the text 
Fig. 5 Initial growth rates of the extraordinary a) and ordinary b) waves in the presence of a thermal electron component. Simulation parameters correspond to the model 9 (Table 1), and the cold plasma dispersion relation is used. 

Open with DEXTER  
In the text 
Fig. 6 Electron distribution function a) and the Zmode energy density b) at the late relaxation stage. The thermal component of the distribution function is shown only partially. Simulation parameters correspond to the model 9 (Table 1). 

Open with DEXTER  
In the text 
Fig. 7 Time evolution of the Xmode waves for the model 9 (Table 1). The contour levels are 0.125, 0.25, 0.5, and 0.9 of the maximum wave energy density. 

Open with DEXTER  
In the text 
Fig. 8 Time profiles of the average growth rates of the Z and Xmodes for the model 9 (Table 1). For the Xmode, the scale differs from that for the Zmode to enhance the growth rate variations. 

Open with DEXTER  
In the text 
Fig. 9 Energy conversion efficiency and the characteristic timescales of the relaxation process for the different concentrations of the energetic electrons. The values correspond to the models 5, 9, and 19 (Table 1). 

Open with DEXTER  
In the text 
Fig. 10 Same as in Fig. 9, for the different values of the total plasma density (concentration of the accelerated electrons n_{b} is constant). The values correspond to the models 15, 9, and 6 (Table 1). 

Open with DEXTER  
In the text 
Fig. 11 Same as in Fig. 9, for the different energies of the accelerated electrons. Solid lines: models 13, 15, and 18 (without thermal component); dashed lines: models 7, 9, and 12 (with thermal component). 

Open with DEXTER  
In the text 
Fig. 12 Same as in Fig. 9, for the different angular distributions of the accelerated electrons. Solid lines: models 14, 15, 16, and 17 (without thermal component); dashed lines: models 8, 9, 10, and 11 (with thermal component). 

Open with DEXTER  
In the text 
Fig. 13 Same as in Fig. 9, for the different strengths of the magnetic field. Solid lines: models 2, 4, and 15 (without thermal component); dashed lines: models 1, 3, and 9 (with thermal component). 

Open with DEXTER  
In the text 