Field-aligned particle acceleration on auroral field lines by interaction with transient density cavities stimulated by kinetic Alfvén waves

. We consider the ﬁeld-aligned acceleration of energetic ions and electrons which takes place on auroral ﬁeld lines due to their interaction with time-varying density cavities stimulated by the strong oscillating ﬁeld-aligned currents of kinetic Alfv´en waves. It is shown that when the ﬁeld-aligned current density of these waves increases, such that the electron drift speed exceeds the electron thermal speed, ion acoustic perturbations cease to propagate along the ﬁeld lines and instead form purely-growing density perturbations. The rarefactions in these perturbations are found to grow rapidly to form density cavities, limited by the pressure of the bipolar electric ﬁelds which occur within them. The time scale for growth and decay of the cavities is much shorter than the period of the kinetic Alfv´en waves. Energetic particles traversing these growing and decaying cavities will be accelerated by their time-varying ﬁeld-aligned electric ﬁelds in a process that is modelled as a series of discrete random perturbations. The evolution of the particle distribution function is thus determined by the Fokker-Planck equation, with an energy diffusion coefﬁcient that is proportional to the square of the particle charge, but is independent of the mass and energy. Steady-state solutions for the distribution functions of the accelerated particles are obtained for the case of an arbitrary energetic particle population incident on a scattering layer of ﬁnite length along the ﬁeld lines, show-ing how the reﬂected and transmitted distributions depend on the typical “random walk” energy change of the particles within the layer compared to their initial energy. When this typical energy change is large compared to the initial energy, the reﬂected population is broadly spread in energy about a mean which is comparable with the initial energy, while the transmitted population has the form of a strongly accelerated ﬁeld-aligned beam. We suggest that these processes are responsible for the occurrence of accelerated ﬁeld-aligned beams of ions and electrons that are commonly observed on auroral ﬁeld lines in planetary magnetospheres.


Introduction
It is now well established that charged particles are accelerated by a number of specific processes in planetary magnetospheres, contributing to the hot plasma populations that they contain, as well as to auroral precipitation and the generation of electromagnetic emissions. Such processes include the conversion of magnetic energy in the current sheets downstream from reconnection sites on the magnetopause and in the magnetic tail, field-aligned acceleration by large-scale parallel electric fields in regions of strong field-aligned current, resonant interactions with electromagnetic waves excited in the plasma, and stochastic acceleration due to plasma turbulence. Regions of particular importance for particle acceleration are those containing strong background fieldaligned electric currents, which map into the auroral ionosphere at the Earth. Observations by spacecraft, such as Freja, FAST, and Polar, at altitudes ranging from 2000 km to 12 000 km, have shown that these regions are characterized by specific features, in particular (a) low background plasma density such that the condition ω pe <ω Be is satisfied, where ω pe is the electron plasma frequency and ω Be the electron cyclotron frequency, (b) non-isothermal plasma satisfying T e >T i , as indicated by the existence of ion acoustic waves which only propagate under this condition (Jayasree et al., 2003), and (c) the presence of kinetic Alfvén waves with wave vectors nearly perpendicular to the ambient magnetic field and with small spatial scales across the field, typically k ⊥ c/ω pe ∼1 (Genot et al., 1999(Genot et al., , 2000. Kinetic Alfvén Published by Copernicus GmbH on behalf of the European Geosciences Union.  1. Sketch of the energetic particle acceleration layer on auroral field lines, associated with growing and decaying density cavities along the field direction that are stimulated by the varying field-aligned currents of kinetic Alfvén waves propagating nearperpendicular to the magnetic field B 0 . The latter waves are indicated by the propagation vector k. waves are associated with instabilities in regions where fieldaligned currents vary strongly across the field lines (Kozlovsky and Lyatsky, 1997;Wu and Seyler, 2003), and have been studied intensively, both experimentally and theoretically (Louarn et al., 1994;Carlson et al., 1998;Chaston et al., 1999Chaston et al., , 2003Chaston et al., , 2004Stasiewicz et al., 2000;Wygant et al., 2000;Paschmann et al., 2002;Wu and Chao, 2004). These studies have shown that the kinetic Alfvén waves are accompanied by spikes in the electric field directed both parallel and perpendicular to the ambient magnetic field and associated decreases in plasma density with n/n∼50% and transverse scales of the order of ∼c/ω pe (Louarn et al., 1994;Chaston et al., 1999;Genot et al., 1999;Wu and Chao, 2004).
The region on auroral field lines containing kinetic Alfvén waves and density cavities with strong electric fields evidently plays a major role in the acceleration of charged particles. The physical nature of the acceleration processes are of great importance in magnetosphere physics and have been intensively studied in the literature (Wu and Chao, 2004;Genot et al., 2004;Chaston et al., 2004). Observations, e.g. from the Freja and FAST spacecraft, have shown that electrons and ions are accelerated along the field lines to energies of ∼1-10 keV (e.g. Khotyaintsev et al., 2000;Chaston et al., 2000;Wygant et al., 2002). The accelerated electrons form counter-streaming (bi-directional) beams directed almost along the auroral field lines, while the accelerated ions form broader distributions (Chaston et al., 2004). Evidence exists that similar wave processes may also operate in Jupiter's middle magnetosphere (Saur et al., 2003), which is connected magnetically to the Jovian main auroral oval. Ex-tended regions containing bi-directional distributions of accelerated ions and electrons have also been observed in the outer parts of Jupiter's magnetosphere, in this case at energies from a few 10 s of keV to a few MeV (Lanzerotti et al., 1993;Staines et al., 1996).
In this paper we thus consider the processes that are responsible for the field-aligned acceleration of charged particles in association with kinetic Alfvén waves and density cavities, according to the following scheme. First, the field-aligned currents that occur in kinetic Alfvén waves are shown to provide conditions for the instability of small-scale density variations along the magnetic field lines. Development of the instability leads to the formation of small-scale transient density cavities with intense time-varying electric fields. Charged particle interaction with these cavities then changes the particle energy, and can provide an effective diffusive acceleration of the particles along the field lines. In this paper each of these topics is discussed in turn. In Sect. 2 it is shown that the field-aligned currents of kinetic Alfvén waves can become very large if they propagate nearly perpendicular to the field lines, and that under such circumstances the plasma will become unstable to the formation of density cavities which arise out of acoustic perturbations propagating along the field lines, as sketched in Fig. 1. Following the earlier works of Misonova (2001, 2002), in Sect. 3 we then consider the acceleration of charged particles through their interaction with the time-varying electric fields of such growing and decaying cavities, modelling the interaction as a series of discrete random perturbations that can be described by the Fokker-Planck equation. Solutions for the particle distribution function are obtained for the case in which a steady source of particles, ions or electrons, is incident upon a scattering layer of finite width, it being shown that the transmitted population forms an accelerated field-aligned beam under appropriate circumstances. We thus suggest that these processes are responsible for the accelerated field-aligned populations of ions and electrons which have been observed on auroral field lines in both the terrestrial and Jovian magnetospheres.

Field-aligned currents in kinetic Alfvén waves
The essential initial feature of the system proposed here is the existence of large-amplitude kinetic Alfvén waves on high-latitude auroral field lines and their associated oscillating field-aligned electric currents. As we will see in the sections below, if these currents are of sufficient amplitude, they lead to the growth of density cavities along the auroral field lines. We thus begin by examining the properties of kinetic Alfvén waves, and consider a plane sinusoidal wave ∼ exp i (ωt−k.r), which propagates at angle θ with respect Ann. Geophys., 24, 2313-2329, 2006 www.ann-geophys.net/24/2313/2006/ to a background magnetic field B= (0, 0, B 0 ), such that k=k (sin θ, 0, cos θ) = (k x , 0, k z ). The wave fields (indicated in the following by a tilde) are correspondingly given byẼ= Ẽ x , 0,Ẽ z andB= 0,B y , 0 , related by Faraday's law. The current density in the wave is given in terms of the wave magnetic field by Ampère's law, since the displacement current density may be neglected when the Alfvén speed V A is much less than the speed of light, as assumed here. The field-aligned current densityj z is thus given by where we note that the x component of the current density is simply related through divj =−i k·j =0. The value of k x is given by the dispersion relation for kinetic Alfvén waves (Streltsov and Lotko, 1995;Lysak and Lotko, 1996) where v T e is the electron thermal speed, assumed small compared with V A , and ω BH is the hybrid gyrofrequency equal to the geometric mean of the ion and electron gyrofrequencies ω Bi and ω Be , respectively, i.e.
Setting k z =k x cot θ , Eq.
(3) becomes a quadratic equation for k 2 x for a wave of given angular frequency ω and propagation angle θ. The root of physical interest has k 2 x ≥0 for all θ , with |k x | increasing monotonically from zero at θ =0 to infinity at θ =π/2. In relation to Eq. (2), this implies that the field-aligned current densityj z increases monotonically with θ for a wave of given magnetic amplitude, such that the waves of interest here with large field-aligned current density, are those propagating near-perpendicularly to the background magnetic field. Of course in practice, |k x | cannot increase indefinitely, but is limited to a value of the order of ∼ω pe /c (Genot et al., 1999(Genot et al., , 2000, where ω pe = e 2 N/ε o m e is the electron plasma frequency. For large |k x | Eq. (3) gives k z ≈ω/v T e , such that the corresponding limitation on the propagation angle is tan θ≤ ω pe v T e /ω c . From Eq. (2), the limiting field-aligned current density for waves of given amplitude B y is then given by µ ojz ≤ B y ω pe /c.

Condition for instability of the field-aligned current leading to density cavity formation
We now consider the stability of the regions of strong fieldaligned current in the kinetic Alfvén waves to the formation of small-scale density cavities. We note that although the current densityj z in principle varies along the field lines as cos (ωt−k z z+ϕ), where ϕ is the phase on a particular field line, here we consider acoustic perturbations whose spatial and temporal scales are both small in comparison with the kinetic Alfvén wave. In this case we can treat the kinetic Alfvén wave current densityj z as an effective constant in the analysis. The quasi-hydrodynamic equations governing the acoustic perturbations are as follows (e.g. Dungey, 1958). The continuity equation is where ρ=m i N i +m e N e is the mass density of the plasma, and u is the centre-of-mass velocity given by (m i +m e ) u=m i v i +m e v e . In these expressions the ion and electron masses are m i and m e , respectively, their number densities are N i and N e , and their individual bulk speeds are v i and v e . The momentum equation can then be written as where F is the force per unit volume on the plasma, the dyadic tensor (uu) kl =u k u l , ρ q =e (N i −N e ) is the charge density (assuming singly-charged ions), where e is the magnitude of the electron charge, j =e (N i v i −N e v e ) is the current density, and the elements of the mechanical stress tensor S are given by where P is the sum of the ion and electron pressures, each assumed isotropic in their individual rest frames. Assuming that m i m e , N i m i N e m e , and |j | ρ q u , we then have the approximate form that will be used in subsequent calculations S kl ≈ P δ kl + m e e 2 N e j k j l .
We now consider one-dimensional acoustic perturbations that vary in the z direction along the background field, with particle velocity perturbations that occur also only in the z direction. From Eqs. (6) and (7) the force per unit volume on the plasma then only has a z component, given by where we have used the second approximate form of Eq. (7) and Gauss's law (divE=ρ q /ε o ) in Eq. (6). We further assume that in the initial stages of development the plasma remains quasi-neutral with N i ≈N e =N , such that the electric field term can be neglected, and that j z remains constant at the valuej z determined by the kinetic Alfvén wave. We then find where we have also assumed that the plasma pressure obeys the polytropic law (P /N γ ) =const in the perturbations. It can be seen that if the current density is small such that the pressure term dominates, then F z and (∂N /∂z) will be of opposite sign. If we then consider some perturbation of the density away from its initial value, the direction of the force will be away from the peaks of the perturbation, and towards the troughs, thus returning the plasma to the initial condition. If, however, the current density is sufficiently large that the first term dominates in Eq. (9), i.e. if then F z and (∂N/∂z) will be of the same sign, such that an initial perturbation will grow, leading to the formation of density cavities in the plasma. The instability condition given in Eq. (10) is essentially such that the electron speed along the field lines associated with the field-aligned current v e should exceed the electron thermal speed v T e . For a given field-aligned current densityj z , Eq. (10) can also be regarded as a condition that the background plasma density is sufficiently low. Using the cold ion approximation T e T i , we have P ≈N κ B T e , where κ B is Boltzmann's constant, such that the condition given by Eq. (10) becomes N <N th , where the threshold density is The corresponding condition on the amplitude of the kinetic Alfvén wave which drives the current is given by Eq. (2) We pointed out in Sect. 2.1 above that for kinetic Alfvén waves of given angular frequency ω, |k x | monotonically increases from zero for field-aligned propagation, to a maximum value of ∼ω pe /c for nearly field-perpendicular propagation (Genot et al., 1999(Genot et al., , 2000. Thus, waves whose amplitude satisfies B 2 y >µ o γ P will exceed the threshold field strength B th if they propagate sufficiently perpendicular to the magnetic field lines. It can be seen from Eq. (12) that instability is favoured in regions of low background plasma density. In fact, the processes considered here preferentially take place in regions of comparatively rarefied plasma, where ω pe <ω Be .

Time scale for density cavity growth
We now consider the development of the ion acoustic (ion sound wave) instability in more detail. Putting ρ≈N m i in Eqs. (5) and (6), the one-dimensional continuity and momentum equations become and where we have again neglected the electric field term in Eq. (14). We have also put where subscript "o" indicates initial unperturbed quantities, again assuming that the pressure perturbations obey the polytropic law (P /N γ ) =const (noting that Eq. 15 becomes P =κ B T eo N for γ =1), and that T eo T io , such that the ion sound waves are not strongly Landau damped by the ions. From Eq. (11) we also have wherej z is the current density of the kinetic Alfvén wave, taken to be essentially a constant on the time scale of the ion sound waves, as indicated above, and is the ion sound (acoustic) speed whenj z =0 (see later). A typical numerically computed solution of Eqs. (13) and (14) is shown in Fig. 2. In this case we have taken (N th /N o ) =1.05, such that the instability threshold Eq. (10) is slightly exceeded, and consider initial density and velocity perturbations given by N (z, 0)/N o =1−0.063 cos (πz/δ S ) and u z (z, 0)/V So =0.001 cos (π z/δ S ), such that one wavelength of the perturbation corresponds to (z/δ S ) =2. For demonstration purposes we have also assumed for simplicity that γ =1, corresponding to the isothermal approximation, though the results are not sensitively dependent on this choice (the actual value of γ is determined by the relation between the characteristic time scale of the process considered and the time scale for heat conduction in the medium). In Figs. 2a and b we show the normalised density and velocity plotted versus (z/δ S ) at the initial time (blue curves), and also at normalised times (V So t/δ S ) =0.6 (green), and 0.9 (red). It can be seen that the initial negative perturbation in the density grows strongly to produce a density cavity on the order of unit normalised time, associated with modest plasma velocities away from the cavity centre.
We now examine Eqs. (13) and (14) to elucidate this behaviour. If we differentiate Eq. (14) with respect to z, neglect the second term on the LHS containing u 2 z , and employ Eq. (13), we obtain the following equation governing the density Assuming plane sinusoidal waves ∼ exp i (ω S t−k S z), we thus have the ion sound wave dispersion equation It can thus be seen that as the current density in the kinetic Alfvén wave grows, such that N th increases slowly towards N o , the phase speed of the ion sound waves decreases from V So towards zero. More specifically, since N th is a slowlyvarying function of time, but is almost independent of position z along the field lines, the solutions of Eq. (19) are such that the angular frequency ω S of a wave decreases with time as the current increases, while the wave-number k S remains nearly constant. For that part of the kinetic Alfvén wave cycle in which N th >N o (as in Fig. 2), however, Eq. (19) then implies instability, with purely growing density perturbations whose growth rate is given by It can be seen that the largest growth rate corresponds to the largest value of k S , or equivalently, the smallest wavelength. Prior to instability, it can be seen from Eq. (20) that the largest value of k S corresponds to the maximum angular frequency ω S of the ion sound waves, which is limited to ω S <ω pi in the linear regime, where ω pi is the ion plasma frequency (ω pi = e 2 N /ε o m i ). The limit on the wave vector is thus given by which also limits the growth rate given by Eq. (21), since as we have just indicated, the value of k S remains constant during the evolution of the wave as the current density changes. The time scale for growth of the density cavities near threshold can be estimated directly from Eq. (18). We first assume that N th grows only slowly compared with the time scale for growth of the perturbations, such that we may put N th ≈N o . Then writing ∂/∂t→1/τ S and ∂/∂z→1/δ S , where δ S is the spatial scale of the perturbations along the field as in Fig. 2, and retaining only the dominant second term on the RHS of Eq. (18) for N<N o , we obtain This approximation is in accord with the results of numerical investigations, such as those shown in Fig. 2. Again, Eq. (23)

Fig. 2.
Formation of a density cavity as determined by numerical integration of Eqs. (13) and (14), with (N th /N o ) =1.05 and γ =1. The initial conditions are N (z, 0)/N o =1−0.063 cos (πz/δ S ) and u z (z, 0)/V So =0.001 cos (πz/δ S ). Each parameter is plotted versus (z/δ S ) at three times, corresponding to (V So t/δ S ) equal to zero (blue lines), 0.6 (green), and 0.9 (red). The plots show (a) the normalised plasma density (N/N o ), (b) the normalised plasma velocity (u z /V So ), and (c) the normalised electric field (E z /E zo ), where E zo = m i V 2 So /2eδ S . The latter parameter has been determined from Eq. (28) with the neglect of the first (time-derivative) term on the RHS, and where the last term on the RHS becomes logarithmic when γ =1.
implies that the fastest growing perturbations are those with the shortest scale lengths δ S along the field. From Eq. (22) this is limited by giving The growth of the density cavities is eventually terminated by the growth of the electric field term in Eq. (6), which has been neglected above in Eqs. (14) and (18). The electric field can be determined by consideration of the electron fluid equations equivalent to Eqs. (5) and (6). For one-dimensional perturbations along the field, the momentum equation is where v ez is the electron bulk velocity along the field, and the electron pressure is P e =κ B T eo N γ e /N γ −1 o (or P e =κ B T eo N e for γ =1) similarly to Eq. (15) above. We also recall that the current carried by the electrons along the field lines associated with the kinetic Alfvén waves is constant on the time scale of the acoustic instability, so that the continuity condition becomesj z =−eN e v ez ≈const. Employing this in the second term on the LHS, we have where we have used Eq. (16), assuming that N th ≈N o near the threshold of instability, and also Eq. (17). Substituting Eq. (27) and the expression for P e into Eq. (26) then gives where we note that the second term on the RHS effectively corresponds to the gradient of the scalar potential of the electric field. The electric field corresponding to the numerically computed solution in Fig. 2a is shown in Fig. 2c. Here the electric field has been normalised to m i V 2 So /2eδ S , the first term on the RHS of Eq. (28) has been neglected, and the second term within the bracket on the RHS becomes 2 log (N e /N o ) when γ =1, as assumed in this case. It can be seen that the electric field is bipolar in form, accelerating the electrons on entry to the cavity, and decelerating them on exit. Since Eq. (16) at the threshold of instability N th ≈N o implies m i V 2 So ≈m e v 2 ezo , where v ezo is the unperturbed electron bulk velocity associated with the field-aligned current, it can also be seen that the change in the electron energy in the cavity ∼eE z δ S can be by a significant factor, depending on the depth of the density depression in the cavity.
The ultimate depth of the density cavity is limited by the growth of the electric field pressure term in Eqs. (6) and (8), as earlier neglected. From Eq. (6), cavity growth ceases when where N e min is the equilibrium electron density within the cavity. Setting ∂/∂z→1/δ S in Eq. (28) and retaining only the leading (first) term in the potential gradient we find Substituting Eq. (30) into Eq. (29), and using Eq. (16) forj z with N th ≈N o near threshold then yields The minimum electron density is thus about one-half of N o for the minimum cavity size given by Eq. (24), and decreases for cavities of larger size as δ −2/3 S . Neglect of the electric field pressure in Eq. (14), as assumed earlier in this section, is thus a valid approximation for densities somewhat larger than this in the cavity-formation process.
The results derived in this section thus show that the fieldaligned currents of kinetic Alfvén waves of sufficient amplitude propagating nearly perpendicular to the magnetic field can provide the conditions for a density instability to take place in the plasma. This instability then leads to the formation of small-scale density cavities along the field lines, which possess strong field-aligned electric fields. Charged particles which pass through these cavities then undergo changes in energy that can lead to significant overall accelerations of the particles along the field lines. Examination of these energy changes is the goal of next section of the paper.

Energetic particle interaction with an individual density cavity
We now consider the interaction of high-energy particles with the density cavities analysed in Sect. 2, and the nature of the accelerated distributions thereby produced. In this subsection we first estimate the change in the particle kinetic energy along the field lines δε which occurs when such a particle interacts with one growing cavity as the current rises. In general this is given by Ann. Geophys., 24, 2313-2329, 2006 www.ann-geophys.net/24/2313/2006/ where q is the particle charge (positive or negative), E z is given by Eq. (28), t (z) is determined by the particle equation of motion, and the integral is taken over the density cavity in the direction of particle motion. Here we consider particles whose speed along the field lines is much larger than those in the background plasma, such that t can be taken to be constant in the integral. In this case integration of the second (spatial gradient) term in Eq. (28) over the whole cavity clearly goes to zero, such that the only term in the electric field to produce a net effect is the first, involving the partial time derivative of the electron bulk velocity. Thus, the acceleration process discussed here requires time-dependency of the cavities produced through the time-varying currents of the kinetic Alfvén waves, and will not occur for quasisteady cavities driven by quasi-steady currents. Substituting the time-derivative term from Eq. (28) into Eq. (32), and writing ∂/∂t→1/τ S then yields where the sign of δε depends on the sign of the charge q and the direction of the energetic particle motion across the cavity. For a growing cavity in which the current-carrying electron speed increases as the density falls, the electric field given by the first term in Eq. (28) is in the same direction as the field-aligned current, accelerating energetic electrons that move across the layer in same sense as the current-carrying electrons, and retarding those moving in the opposite direction. For ions, the sense is reversed. For decaying cavities, for which the electric field is directed opposite to the fieldaligned current, the opposite signs are appropriate. Setting which is expected to be comparable with the kinetic energy associated with the bulk motion of the plasma electrons, and significantly smaller than the energy of the particle considered. We note that the magnitude of these energy exchanges are independent of the mass and sign of charge of the energetic particles, assuming singly-charged particles, and also of their energy.

Determination of the energetic particle distribution function
Having estimated the change in energetic particle kinetic energy along the field which occurs in the interaction with a single density cavity, δε, we now consider the result of random particle interactions with a set of cavities along the magnetic field lines, as illustrated schematically in Fig. 1. As we have shown in Sect. 2, the time scale for growth of individual cavities (Eq. 23) is much shorter than the period of the kinetic Alfvén waves, and the time scale for their subsequent decay is probably comparable. Considering, therefore, that energetic particles can interact with either growing or decaying cavities when the field-aligned current exceeds the threshold, we assume that the change in particle kinetic energy which occurs in the interaction with one cavity is +δε or −δε with equal probability, depending on whether the density cavities are growing or decaying and the direction of particle motion along the field. As above, we assume that the particles' energy ε significantly exceeds the change in energy in a single cavity interaction, such that ε δε. We also assume that the velocity of the cavities can be neglected, and that the length of the density cavities along the field, δ S , is much less than the distance between them, d, such that the particles can be taken to experience a series of discrete random perturbations. In this case the interaction of the particles with the density cavities can be described by the Fokker-Planck equation. For simplicity we restrict our attention to the steady-state case only, such that the Fokker-Planck equation can be written as where f is the particle distribution function, and is the energy diffusion coefficient. From Eq. (34) this may be estimated as independent of the particle energy, where the angle bracket indicates an average over the ensemble of density cavities. Then introducing the spatial variable ξ =D ε z, Eq. (35) can be written as where f + and f − are the distribution functions for particles moving in the positive and negative z directions, respectively. Continuity at ε=0 then requires that and We further assume that the density cavities are confined to a layer of length L along the field lines in the region 0≤z≤L, or equivalently 0≤ξ ≤h, where h=D ε L, and that magnetospheric particles with distribution function f + o (ε) are incident on the layer at ξ =0. Steady-state solutions as determined here are then appropriate, provided that the incident distribution is essentially unchanging on the time scale required for particles to pass through the scattering layer. It is readily shown that a solution of Eq. (38) for f + (ξ, ε) is valid for any values of the constants a, ξ , and ε . For ε >0 it is evident that this represents a steady-state solution in which a source of mono-energetic particles of initial energy ε=ε moving in the +ξ direction is present at ξ =ξ , and spreads increasingly in energy with increasing distance from the source for ξ >ξ . It is also evident that for ε <0, the solution in Eq. (40) goes to zero at all energies at ξ =ξ . Here we require a solution of Eq. (38) in which f + reduces to f o (ε) at ξ = 0, and hence take where the first integral on the RHS gives f o (ε) at ξ =0, and the second, containing an arbitrary function of energy g + (ε), gives zero in this limit. The latter function represents particles moving in the positive ξ direction that have undergone previous reflections. As can be seen from Eq. (41), Eq. (40) is the Green function for the diffusion equation obeyed by f + . Similarly, for the diffusion equation obeyed by f − (Eq. 38), the corresponding solution is valid for ξ <ξ . In this case, the required solutions are those for which f − =0 at ξ =h, since we assume for simplicity that no flux moving in the −ξ direction is incident on the layer from the region ξ >h. We therefore set for some arbitrary function of energy g − (ε). If we now apply the boundary conditions at ε=0 given by Eqs. (39a) and (39b) we find and which are the equations from which g + and g − may be determined for given f o . The calculations of g + and g − are somewhat lengthy, such that the details are given in Sect. A1 of the Appendix. The results are and given by Eqs. (A12b) and (A19b), respectively, where the functions C (z) and S (z) are the cosine and sine Fresnel integrals defined by Eq. (A18). The general solution of our problem for an arbitrary incident particle distribution is thus given by Eqs. (41) and (43), taken together with Eqs. (45) and (46). To determine the distribution functions of the particles that are transmitted through (f T (ε)) and reflected from (f R (ε)) the layer, we thus substitute Eqs. (45) and (46) into Eqs. (41) and (43), respectively. The details of the calculations are given in Sect. A2 of the Appendix, with the results and  (A25b) and (A29b). The functions erf (z) and erfi (z) are the error function and the error function of imaginary argument, respectively, defined by Eq. (A24).
To illustrate these results we choose for simplicity an initial distribution function f o (ε) in which the parallel energy of the particles is mono-energetic at energy ε o , i.e.
Substituting Eq. (49) into Eqs. (47) and (48), the transmitted and reflected distribution functions can be written as and From Eqs. (50) and (51) it can thus be seen that the form of the parallel energy distribution, relative to the initial energy ε o , depends on the dimensionless ratio √ 2h/ε o , where as above h=D ε L=(δε) 2 L/2d. The interpretation of this parameter becomes clear if we consider that the "random walk" change in parallel energy of the particle after n cavity interactions is ε≈ √ n δε. For a particle traversing the scattering layer of width L with cavities separated by distance d along the field, as above, the number of interactions is n≈L/d, so that a measure of the change in particle energy is ε≈ √ L/d δε= √ 2h. Thus, parameter √ 2h/ε o is the ratio of the typical energy change due to cavity interactions for particles traversing the layer, compared with the initial energy. Some results are shown in Fig. 3, where we plot the normalised distribution functions of transmitted (solid lines) and reflected (dashed lines) particles versus (ε/ε o ) for various values of √ 2h/ε o . The vertical dotted lines show the parallel energy of the incident particles (ε/ε o ) =1. In Fig. 3a we show results for a "thin" layer in which The transmitted distribution takes the form of a beam which is somewhat spread about the initial energy, with "tails" that extend to higher and lower energies. The reflected distribution is relatively weak and broadly spread at low energies. Results for an intermediate layer with Fig. 3b. The distributions have now broadened in energy, and the flux of reflected particles has significantly increased. In Fig. 3c we then show results for a "thick" layer with ε/ε o = √ 2h/ε o =100. Here the reflected distribution peaks near (ε/ε o ) ≈1, but is very broadly spread in energy. The transmitted distribution is also Plot (a) is for a "thin" layer in which the typical change in particle energy is small compared with the initial energy, specifically for √ 2h/ε o =0.3. Plot (b) is for a layer of intermediate thickness in which the typical change in particle energy is comparable with the initial energy, specifically for √ 2h/ε o =1. Plot (c) is for a "thick" layer in which the typical change in particle energy is large compared with the initial energy, specifically for www.ann-geophys.net/24/2313/2006/ Ann. Geophys., 24, 2313-2329, 2006 broadly spread, but peaks near ε≈ √ 2h, such that it forms a broad and strongly accelerated particle beam along the field lines. The number of cavity interactions required to produce these distributions is readily estimated from the above argument as n≈ √ 2h/ε o 2 (ε o /δε) 2 . If we take, for example, δε/ε o ≈0.1, then n is approximately ∼10, ∼100 and ∼10 6 for ε/ε o = √ 2h/ε o =0.3, 1, and 100 in Figs. 3a, b, and c, respectively. Rough estimates indicate that values up to ∼10 5 may be realistic in the terrestrial magnetosphere, with even larger values being possible at Jupiter and Saturn.
We finally calculate the average parallel energy of the particles in the transmitted and reflected distributions for the "thick" layer case, for which ε/ε o = √ 2h/ε o 1, whereε o is the average energy of the incident population. In general, the average energy is given bȳ From Eq. (47) the transmitted distribution can be written as where for small argument Substituting Eq. (54) into Eq. (53) we thus find and so substituting into Eq. (52) we havē which is thus in conformity with the simple "random walk" estimate ε≈ √ 2h made above within a factor of order unity. Setting h=D ε L=(δε) 2 L/2d in Eq. (56) thus yields where δε is given by Eq. (34). Since the field-perpendicular energy of the particles remains ∼ε o during the interaction with the scattering layer, Eq. (57) provides a measure of the anisotropy of the particle distribution along and across the field. This ratio increases in proportion to the square root of the width of the layer L, and in direct proportion to the typical energy exchange with a single cavity δε. The latter is proportional to the modulus of the particle charge (Eq. 34), but otherwise depends only on the properties of the cavities themselves. Thus, singly-charged ions and electrons are accelerated along the field to the same energies as each other within the acceleration layer.
For the reflected distribution we substitute Eq. (51) into Eq. (52) and reverse the order of integration to find For ε 2 /4h 1 we then have so that substituting into Eq. (58) we havē The average energy of the reflected population is thus generally of the same order as the energy of the initial distribution.

Discussion and summary
In this paper we have shown how charged particles can be significantly accelerated along magnetic field lines due to their interaction with transient density cavities that are formed in regions of strong varying field-aligned currents associated with kinetic Alfvén waves. It has been shown that if the field-aligned currents associated with these waves become sufficiently large, such that the electron drift along the field lines exceeds the electron thermal speed, then the Ann. Geophys., 24, 2313-2329, 2006 www.ann-geophys.net/24/2313/2006/ plasma becomes unstable to the formation of density cavities. For kinetic Alfvén waves of a given amplitude, the field-aligned current density increases monotonically with the angle of propagation relative to the background field, so that the instability will generally be associated with waves propagating nearly perpendicular to the background field. It has been shown that as the field-aligned current density in the wave increases towards the threshold, the phase speed of ion sound waves propagating along the field direction decreases towards zero. When the instability threshold is exceeded, however, these waves become purely growing density perturbations along the field, the rarefactions of which develop rapidly into density cavities. The initial time scale for cavity growth is comparable to δ S /V So (Eq. 23), where δ S is the perturbation scale length along the field lines and V So is the sound speed for zero field-aligned current density, and is thus shortest for the smallest cavities along the field. However, δ S is limited to ∼V So /ω pi , where ω pi is the ion plasma frequency, so the minimum time scale is limited to ∼ω −1 pi (Eq. 25). The cavities have an associated parallel electric field which is approximately bipolar in form (Eq. 28), and such that the current-carrying electrons are accelerated along the field as they enter the cavity and retarded as they leave. The existence of such cavities on auroral field lines has been reported, for example, in data from the FAST spacecraft (Carlson et al., 1998;Chaston et al., 1999). The growth of the cavities is eventually limited by the pressure of this electric field, such that the density depletion is limited to a factor of about one-half for cavities with the minimum scale length of ∼V So /ω pi , the depletion factor then decreasing as δ −2/3 S for larger spatial scales (Eq. 31).
We then considered the interaction of individual energetic particles with these cavities. Clearly, the "electrostatic" electric field associated with the cavities produces no net effect on the particle energy, such that there is no net acceleration if such particles interact with quasi-steady cavities formed in regions of strong but steady currents. However, when the cavity properties are changing with time, a net acceleration is found to be present which is proportional to the rate of change of the velocity of the current-carrying electrons (Eq. 28). This effect gives rise to increases or decreases in energetic particle energy, depending upon whether the cavities are growing or decaying, and on the direction in which the energetic particles traverse the cavities. The change in energy is proportional to the charge of the energetic particle, but is independent of its mass and energy, such that both ions and electrons are equally affected (Eq. 33).
Overall, these interactions can be treated as a series of discrete random perturbations that can be described by the Fokker-Planck equation (Eqs. 35 and 37). Appropriate solutions have been derived for the case in which a steady source of particles is incident upon a scattering layer of length L along the field, within which the parallel energy of the particles is changed randomly by ±δε at discrete events sep-arated by distance d along the field lines (where d δ S ). General expressions for the distributions of the transmitted and reflected particles have been derived (Eqs. 50 and 51), and have been illustrated for the case of a mono-energetic particle source. The typical change in energy of a particle traversing the layer is given by the "random walk" formula ε∼ √ L/d δε, and the nature of the results depends on whether this is large or small compared with the initial energy of the particles ε o . When ε ε o , the transmitted distribution is weakly-scattered about the initial energy, and the reflected population is weak and broadly-spread at low energies. When ε ε o , however, the transmitted population is strongly accelerated along the field to peak at an energy of ∼ ε, which thus increases as the square root of the length L of the scattering layer and in proportion to the change in parallel energy |δε| for the interaction with an individual cavity. We recall that the latter is proportional to the particle charge, but otherwise depends only on the cavity properties (Eq. 33), such that singly-charged ions and electrons are accelerated along the field to the same energies. Under these circumstances, then, the transmitted particles form a highly anisotropic field-aligned beam along the field lines. The reflected particles are also broadly spread in energy, but have an average energy which is comparable with the initial energy. Of course, the geometry of the calculation presented here is simplistic. If the initial population corresponds to magnetospheric particles incident on the top of a scattering layer, as illustrated in Fig. 1, then the transmitted population is that which emerges underneath, directed towards the ionosphere, which then contributes to the auroral particle precipitation. However, the particles that are mirrored by the strongly increasing magnetic field underneath the layer will then pass back through the layer and be further accelerated along the field, emerging into the magnetosphere as a strongly fieldaligned beam. We suggest that such particles correspond to the field-aligned beams of accelerated particles that are commonly observed on auroral field lines in planetary magnetospheres.
We note that while for simplicity the present work has been based on analytic solutions of a simplified Fokker-Planck equation, important results on related problems have been obtained by other authors using numerical methods (e.g., Isliker et al., 2000Isliker et al., , 2001Arzner et al., 2006). In future work we intend to employ such methods to extend the study presented here to cases where a Fokker-Planck approach is not applicable, and to non-steady situations.
In summary, therefore, the results obtained in the present study are as follows. We have (a) shown that the strong field-aligned currents associated with kinetic Alfvén waves can stimulate the formation of density cavities on auroral field lines; (b) evaluated the main characteristics of the varying density cavities, namely their characteristic time scale, their spatial scale along the field lines, the electric field that occurs within them, and the limiting depletion within the cavities; (c) determined the change in energy that takes place when an energetic ion or electron interacts with an individual time-varying cavity; (d) solved the Fokker-Planck equation for the steady-state distribution functions of energetic particles interacting with an ensemble of growing and decaying cavities within a layer of finite length along the field lines, and have shown that under appropriate circumstances the transmitted population will form a field-aligned beam.
We finally emphasise that the model of auroral fieldaligned particle acceleration proposed here conforms with the principal features observed in regions containing strong field-aligned electric currents, that is the existence of density cavities with strong electric fields in association with kinetic Alfvén waves (Louarn et al., 1994;Chaston et al., 1999), and the correlation of the latter phenomena with accelerated particle fluxes (Wu and Chao, 2004;Genot et al., 2004). However, existing data do not as yet permit a detailed comparison to be undertaken of the actual characteristics of density cavities and accelerated particle fluxes with analytical results such as those derived here. The physical nature of transient density cavities and their role in the formation of accelerated electron and ion fluxes in real magnetospheric conditions is a topic of continuing interest for future studies.
Appendix A A1 Determination of the functions g + (ε) and g − (ε) In this Appendix we determine the functions g + (ε) and g − (ε) in terms of f o (ε) from the boundary conditions at ε=0 given by Eq. (44) which we reproduce here as and To do this, we first writẽ so that Eq. (A1) becomes and We now use the integral expressions and Substituting these into the RHS of Eq. (A3) and reversing the order of integration yields and If we then substitute these expressions into the RHS of Eq. (A3), and compare the integrands with those on the LHS, we find the two relations Subtracting these expressions gives We now use the identities (e.g. Gradshteyn and Ryzhik, 1980, their Eq. 3.751) and and substituting into Eq. (A8), reversing the order of integration, and comparing integrands, gives Thus, the Fourier sine transform of √ ε g + ε is equal to minus the Fourier cosine transform of √ ε f o ε . Performing the inverse sine transform gives Reversing the order of integration, and integrating over the trigonometric functions we finally find or using Eq. (A2) To determine g − (ε) we add together Eqs. (A7a) and (A7b) to find In the second integral we write (e.g. Gradshteyn and Ryzhik, 1980, their Eq. 3.691) and also introducing Eq. (A12a) and rearranging the order of integration yields Now (e.g. Gradshteyn and Ryzhik, 1980, their  where C (z) and S (z) are the cosine and sine Fresnel integrals defined by C (z) = 2 π z 0 dt cos t 2 and S (z) = 2 π z 0 dt sin t 2 .
Substituting Eqs. (A16) and (A17) into Eq. (A15), and the latter into Eq. (A13) then yields or using Eq. (A2) we have finally A2 Determination of the transmitted and reflected distribution functions We now wish to use the results of Sect. A1 to calculate the distribution functions of the particles which are transmitted through and reflected by the diffusive layer. From Eq. (41) the transmitted distribution function is Then substituting forg + (ε) from Eq. (A12a) and reversing the order of integration we have We now use (e.g. Gradshteyn and Ryzhik, 1980, their Eq. 3.896) exp substitute into Eq. (A21) and rearrange the order of integration to find where erf (z) is the error function, and erfi (z) is the error function of imaginary argument (erfi (z) =erf (iz)/i), given by erf (z) = 2 √ π z 0 dt exp −t 2 and erfi (z) = 2 √ π z 0 dt exp t 2 .
Substituting Eq. (A23) into Eq. (A20) then gives the transmitted distribution function as or using Eq. (A2) Similarly, from Eq. (43) the reflected distribution function is Substituting from Eq. (A7a) and reversing the order of integration, we find We now substitute Eq. (A12a) forg + (ε) and reverse the order of integration to find Substituting Eq. (A28) into Eq. (A27), and then into Eq. (A26) we thus find the reflected distribution function to be given by or using Eq. (A2)