A linear auroral current-voltage relation in fluid theory

Progress in our understanding of auroral currents and auroral electron acceleration has for decades been hampered by an apparent incompatibility between kinetic and fluid models of the physics invo ...


Introduction
Numerical simulations are an important tool in studies of auroral phenomena and magnetosphere-ionosphere coupling.Simulations that treat the electrons as particles can describe fast, small-scale phenomena, but fluid models are required for slow (several seconds), large-scale (several R E ) processes that determine the global dynamics.When simulating the large-scale behavior of an auroral flux tube, a major difficulty has been the inability of fluid models to properly describe the generation of field-aligned electric fields and the associated electron acceleration.
Correspondence to: J. Vedin (jorgen.vedin@space.umu.se)Observations from rockets and satellites (Hoffman, 1993;Evans, 1968Evans, , 1974;;Mizera and Fennel, 1977) indicate that precipitating auroral electrons are freely accelerated by a field-aligned potential drop at altitudes around 1 R E .Many features that are observed in the particle distributions are predictable from a one-dimensional description of the kinematics of collisionless electrons in static magnetic and electric fields (Whipple, 1977;Chiu and Schulz, 1978;Lundin and Eliasson, 1991).Kinetic models also show that the current is linearly related to the potential drop if an isotropic Maxwellian equatorial source distribution is assumed (Knight, 1973;Fridman and Lemaire, 1980).Fluid models used in simulations of the upward auroral current region have not been consistent with this kinetic description.While it is well known that shear Alfvén waves with short perpendicular wavelength will produce a field-aligned electric field, the connection between the quasi-static structures considered in kinetic models and the Alfvén waves remains unclear.In particular, it has not been possible to derive a linear current-voltage (C-V) relation from the collisionless fluid models describing Alfvén waves.To obtain a linear C-V characteristic within a fluid model, it has been necessary to introduce anomalous resistivity (Lysak and Dum, 1983;Streltsov et al., 2002), which seems incompatible with the collisionless theory that fits particle observations.The hierarchy of fluid equations is usually closed by assuming a local relation between the density and temperature in the form of an equation of state, but for field-aligned flows in a collisionless plasma far from thermal equilibrium it is difficult to justify such an equation.It seems that all fluid models of auroral electron acceleration, in the absence of a vindicable equation of state, have been based on the assumption that temperature variations can be neglected.This applies to the classical studies by Goertz and Boswell (e.g. 1979) and Lysak and Dum (1983), as well as more recent work (e.g.Rönnmark and Hamrin, 2000;Streltsov et al., 2002).However, in this study we will show that the electron temperature variations, caused by the auroral current, have profound effects on the dynamics of the auroral acceleration region.Although there is no local equation of state, from a kinetic model we can determine how the electron temperatures depend on the field-aligned current.The chain of fluid equations may then be closed by introducing nonlocal equations of state in the form of approximants that describe the temperature variations.We find that the temperature gradients play a decisive role in the electron momentum equation, and that the proper inclusion of these gradients is essential to the derivation of a linear C-V relation from fluid theory.

Theory
Let z be a coordinate along the magnetic field line, with z=0 at the equatorial plane.The boundary of the generator region is at z G , the bottom of the acceleration region at z a and the ionospheric boundary at z I .Introducing the field-aligned velocity v z =ż and the perpendicular velocity v ⊥ , the Vlasov equation can be written (1) The Vlasov equation determines the evolution of the phase space density f =f (v z , v ⊥ , z, t).
Integrating Eq. ( 1) over velocity space, we obtain a continuity equation for the electrons.Introducing the charge density ρ=e(n i −n), where e is the proton charge, n i is the ion density and n is the electron density, in the presence of an inhomogeneous magnetic field B we can write the equation of current continuity as if we neglect for simplicity the ion motion.In this equation j z is the field-aligned electron current density.Multiplying the Vlasov Eq. ( 1) by v z before integrating over velocity space we find the momentum equation, or an equation for the evolution of the field-aligned current density j z .Including an electrostatic potential φ and kinetic temperatures T z and T ⊥ , we find (Rönnmark, 2002) (3) As before we neglect the ion contributions to this equation, since they are small by at least a factor of √ m/m i , where m is the electron mass and m i is the ion mass.A more complete version of Eq. (3) was derived by Mitchell and Palmadesso (1983), who included ions and gravitation as well.
When the phase space density is known, we can calculate fluid quantities such as the electron density the field-aligned current density the perpendicular temperature and the parallel temperature Introducing new independent variables µ and H , where µ is the magnetic moment and H is the total energy we can describe the phase space density by the functions g ± defined by or equivalently Since they describe the phase space density, the functions g ± must satisfy the Vlasov equation, which in these variables takes the form since μ= Ḣ =0.Assuming a stationary state with ∂ t g ± =0 it follows from Eq. ( 11) that g ± is independent of z, and that the phase-space density will be constant along the trajectories defined by Eqs. ( 8) and ( 9).Hence, if we specify the velocity distribution function F G =F G (v z , v ⊥ ) at the generator boundary (z=z G ), where φ=0 and B=B G , this defines along all trajectories that pass through this boundary with v z >0.Similarly, we specify an ionospheric distribution F I that defines f on trajectories that pass the ionospheric boundary at z I with v z <0.
The current density is calculated by determining which particles that can reach the ionosphere.Provided the potential, for all z, satisfies known as the Fridman-Lemaire (F-L) condition (Fridman and Lemaire, 1980), the current depends only on the total potential drop φ along the field lines and is independent of the shape of the potential φ(z).In this equation we have introduced B I ≡B(z I ).The F-L condition is, however, not sufficient to make the density and temperatures independent of how the potential varies along the field line.For this the much more stringent condition (Janhunen, 1999) must be satisfied.Here, we will determine the limits of integration in velocity space by a method that takes the shape of the potential fully into account, and consequently, is independent of conditions ( 13) and ( 14).

Method
It is convenient to consider the limits of integration for the integrals ( 4)-( 7) in the µ-H plane (Whipple, 1977).An electron at the generator boundary, where B=B G and φ=0 must have µ≤H /B G according to Eq. ( 9).The line µ=H /B G , corresponding to v z =0, is the turning point line at the generator (z=z G ).At any altitude we introduce the turning point line μ, defined by Electrons coming from the generator region will always be within the region µ≤ μ and H ≥0, but parts of this region may be inaccessible to downgoing electrons.In order to reach a level z, the electron must have µ≤ μ at all levels between the generator and z.Defining the turning point boundary µ tpb by we find that electrons can reach z, if and only if they have µ≤µ tpb (z, H ). We determine the turning point boundary numerically by the method illustrated in Fig. 1.Starting from µ tpb (z G , H )= μ(z G , H ) we take a small step to z 1 =z G + z, calculate a new turning point line and record the point where µ tpb (z G , H ) and μ(z 1 , H ) intersect.The smaller of µ tpb (z G , H ) and μ(z 1 , H ) defines µ tpb (z 1 , H ). A typical pattern of intersections is shown in Fig. 1, where the shading indicates the accessible region.Continuing this process recursively towards the ionosphere in about one hundred small steps, we obtain a table of intersection points that allows us to determine µ tpb accurately for any z and H . Magnetospheric electrons will reach the ionosphere if they have µ<µ tpb (z I , H ), and we assume that these electrons are lost.In the interior of the flux tube, in addition to the downgoing electrons, there will be reflected magnetospheric electrons with µ tpb (z I , H )<µ<µ tpb (z, H ).
Electrons originating at the ionospheric boundary (with v z ≤0) must have µ≤(H + eφ I )/B I = μ(z I , H ) and H ≥−eφ I .Reaching the altitude z≤z I the ionospheric electrons must have µ≤µ lcb (z, H ), where Numerically we determine the loss cone boundary µ lcb by starting from μ(z I , H ). Applying the method outlined above for µ tpb in the reverse direction we then recursively go from a level z k to z k−1 =z k − z, to build up a table of intersections between µ lcb (z k , H ) and μ(z k−1 , H ). From this table we find µ lcb (z, H ) by interpolation.
Assuming an isotropic velocity distribution at the generator and ionospheric boundaries the distribution function g ± can be written independent of µ as g ± (H ) and the integrals over µ can then easily be evaluated analytically.For example, the density integral (4) can be expressed as a sum of components of the form where g ± represent the distributions of up-or downgoing magnetospheric or ionospheric electrons, and µ max equals µ tpb or µ lcb .In these integrals H min is the H -value at µ=0 for the desired boundary.The remaining integrals over H are evaluated numerically.

Auroral flux tube model
In this study we use an isotropic Maxwellian velocity distribution at the generator boundary.In Eq. ( 19) N G is the density and T G is the temperature in energy units.In the results we have used the numerical values N G =1 cm −3 and T G =1 keV.Changing variables to H and µ the distribution (19) corresponds to and the volume element in velocity space is expressed as Similarly at the ionospheric boundary we specify an isotropic Maxwellian with temperature T I =10 −3 T G and density All the results in this study are based on an electrostatic potential φ(z) of the form where φ is the main potential drop.How the parameter p is chosen is discussed more thoroughly in the next section.This form of the potential is well suited to model the observed maximum in the field-aligned electric field at altitudes around 1 R E (Reiff et al., 1993;Hull et al., 2003).
The bottom of the acceleration region is often located at altitudes around z a =0.7 R E (McFadden et al., 1999), which corresponds to a magnetic field strength B a =B(z a )=200B G at the lower boundary of the main potential drop.At lower altitudes we have only an ambipolar potential φ amb to ensure that the ionospheric electrons are confined to altitudes below z a also when φ=0.The numerical values of the potential parameters used in the results are φ amb =9 V and φ=10 kV unless otherwise stated in a specific figure.The variation of the magnetic field is for z<z I =9 R E given by B(z By choosing isotropic distributions at the boundaries we have also specified the electron density along the field lines when φ=0.We assume that the ion density is equal to the electron density when φ=0.When the current is increased a new equilibrium between the electrons and the potential will be established on a time-scale characterized by the electron transit time, which is 1-10 s.The ions are too heavy to move significantly during this time, and the electron density must remain roughly constant to keep the plasma quasineutral.After a few minutes the ion density may be affected by processes related to the presence of the field-aligned current.However, a self-consistent treatment of the ions is beyond the scope of this study, and we will assume that the ion density remains fixed.To maintain quasi-neutrality we must then demand that the electron density is essentially independent of the magnitude of φ. (2002) argued that it was necessary to abandon the assumption that the source distribution of electrons should be independent of the current, in order to compensate for the electron density increase near z a caused by the potential drop.However, the results in Rönnmark (2002) were obtained with a potential that violates assumption ( 14), using the classical integration boundaries μ(z, H ) derived from this assumption.The difference between the locally determined boundary μ(z a , H ) and the true boundary µ tpb (z a , H ) is illustrated in Fig. 2. Electrons with (H, µ) within the shaded region are reflected at higher altitudes, although they would have reached z a if the potential had satisfied (14).Figure 3 shows the corresponding region in velocity coordinates.The thinner, solid lines show contours of the distribution function and the dashed, bold line at the innermost contour is the integration boundary corresponding to μ(z a , H ). Considering that the volume element increases in proportion to v ⊥ , it is clear that the shaded area corresponds to a significant fraction of velocity space.Since this volume is mainly at v ⊥ v z , the use of μ(z a , H ) instead of µ tpb (z a , H ) will underestimate the parallel temperature and overestimate the perpendicular temperature and density.When the correct boundaries are used for the velocity space integrals it becomes possible to keep the electron density independent of the field-aligned current without modifying the source distribution.

Rönnmark
Observations strongly indicate that most of the potential drop occurs at altitudes around 1 R E , but the constraint (14) requires a sufficient potential drop at higher altitudes.It is easy to see that Eq. ( 14) is marginally satisfied if φ is proportional to B (Janhunen, 1999).Adopting the model potential Eq. ( 23) we find that Eq. ( 14) is violated when p>1.The dependence of the density on the value of p for a potential drop φ=10 kV is illustrated in Fig. 4. For p=1, when condition ( 14) is satisfied, we find at altitudes around 1 R E a substantial electron density enhancement that would violate quasi-neutrality.When p is increased the electron density is reduced until for p=2 it is comparable to the ion density.Figure 5 shows the electron density for p=2 and different values of φ, and demonstrates that the choice p=2 makes the electron density rather insensitive to the magnitude of the total potential drop.Considering, for example, the uncertainty of the assumed source distribution and background ion density we refrain from searching for an optimal potential shape that would satisfy quasi-neutrality even better, and use the model potential ( 23) with p=2 for this study.
Figure 3 also indicates why the current is much less sensitive to the shape of the potential than the density.The current is carried by electrons in the loss cone, and the shaded area between the p=1 and p=2 contours is almost entirely outside the loss cone.This illustrates that the F-L condition ( 13) is essentially satisfied even when p=2.Since the F-L condition is satisfied, the current-voltage relation is close to linear and can be approximated by (Fridman and Lemaire, 1980) where j z I is the field-aligned current density at the ionospheric boundary.Figure 6 presents the calculated currentvoltage relation showing that the linear approximation is satisfactory for φ, up to at least 20 kV.By evaluating integrals ( 6) and ( 7) we find the temperatures T ⊥ and T z .From a distribution function that satisfies the time independent Vlasov equation we have then calculated all the fluid quantities in the momentum Eq. ( 3).If we then rewrite the stationary version of this equation using the equation of current continuity (2), we obtain In a collision dominated medium the temperatures are locally related to the density by an equation of state, and Eq. ( 25) then gives a local relation between the electric field and the current density.In a collisionless auroral flux tube the plasma properties at different altitudes are tightly coupled by electrons moving along the field lines, and there will be no local equation of state.However, if we globally can approximate the temperatures by functions of B and φ, in a fluid model these approximations can replace the equation of state and provide a local current-voltage relation.
We introduce normalized variables B=B/B I , and =e φ/T G .To find approximations for the temperatures we first consider the region from the generator boundary to the bottom of the acceleration region.Since B is a monotonous function of z we can use B instead of z as the independent variable, and assuming n≈N G in the region B<B a ≡B a /B I we rewrite Eq. ( 25) with the pressures P z and P ⊥ instead of the temperatures to find To obtain this equation we have substituted j z = j z I B, where j z I is expressed in terms of using the linear approximation Eq. ( 24).We have also inserted φ from ( 23) with p=2.The coefficients α and β are , and β = 2 B −2 a + α .
Expressing P ⊥ and P z as power series in B≤B a <1 with coefficients depending on we obtain according to the derivation presented in Appendix A. In this derivation we have also used analytical expressions for the pressure at =0 to determine some of the The remaining coefficients q i are calculated using a least-square fit to the pressure profiles at different .The coefficients are thus determined to be q 11 = 0.138 B −1 a q 12 = −0.0191B −1 a q 21 = 0.556 B −2 a q 22 = 0.0368 B −2 a q 31 = 0.335 B −3 a In the region B∈[B a , B I ], where B I =1, we approximate the pressure with a second degree polynomial in 1−B.The polynomial is chosen to have certain values at B=B a , B=0.5>BI and B=B I , according to Appendix B. These choices are made merely to make the approximation as good as possible.
For both the parallel and the perpendicular pressures we now have two polynomials valid as approximations of the pressure in two different regions Here the coefficients c i are chosen equal to those in the first or second equation in Eq. ( 27), depending on which of the pressure components we want to approximate.The coefficients d i are presented in Appendix B for both the parallel and the perpendicular pressures.Two polynomial approximations valid in two different regions may cause problems when constructing algorithms for solving the fluid equations.In numerical simulations it is more convenient to have a uniform, continuously differentiable approximant for the pressure.To achieve this we continue by finding a Padé approximant approximating the pressure profile in the whole region B ∈ [B G , B I ].In the Padé approximant, Q(B, ) and R(B, ) are polynomials with as low an order as possible but still fulfilling P=P in a satisfactory way.For the perpendicular pressure we choose the polynomials Q(B, ) and R(B, ) to be of order three, and for the parallel pressure we choose fourth order polynomials.The expressions for the coefficients in the Padé approximants are given in Appendix B.
We now have approximations for the pressure in the entire region B∈[B G , B I ] and an approximation for the plasma temperature can then be calculated as In Figs. and 8 we present the parallel and perpendicular temperature profiles, respectively.The profiles obtained by numerical integration of Eqs. ( 6) and ( 7) are in these figures compared with the corresponding approximants for two different φ.
Notice that although we have written the approximants in terms of the potential drop, by using Eq. ( 24) we can replace φ by j zI or the local current density j z =j zI B. Using the approximants (30) we can rewrite the momentum Eq. ( 25) as This equation allows us to calculate a local electric field and a total potential drop when the current density is known.As an example we show in Fig. 9 comparisons between the potentials derived from Eq. ( 31) and the model potentials at the three current densities j z I =−1.7, −5.1 and −9.5 µA, corresponding to φ≈1, 5 and 10 kV.We see that the linear C-V relation is accurately reproduced by Eq. ( 31), and that large electric fields appear at the expected altitudes.Within a thin layer near the bottom of the acceleration region the main contribution to the parallel electric field is from the parallel density gradient, T z ∂ z ln n.This is consistent with the observations of McFadden et al. (1999) and Hull et al. (2003).At higher altitudes, where the density is rather constant, the electric field is supported mainly by the parallel temperature gradient.In general, the second term, involving the temperature anisotropy and magnetic field gradient, gives a smaller but still significant contribution to the electric field.In comparison to these terms, the electron inertia term is almost negligible.
Figure 9 reveals that higher order approximants are needed to obtain an accurate fit to the sharp edge of the model potential at z a .However, even with a higher order model we cannot expect a fully realistic description of details in the acceleration region.Observations indicate that there may be structures with strong parallel electric fields that are highly localized in altitude (McFadden et al., 1999).Even if quasineutrality locally may be violated within such regions, these strong fields must still add up to a potential that maintains quasi-neutrality on a global scale.Nonlocal equations of state should then still be able to describe large-scale temperature variations and represent a substantial improvement to the momentum equation in fluid simulations of electrodynamics of auroral flux tubes.
The electron temperatures in a current carrying auroral flux tubes are in this study approximated by the functions T ⊥ and T z .In a fluid model, if we use these approximations as equations of state, the electron momentum equation provides a local relation between the field-aligned electric field and the current.When applied to the entire flux tube, this momentum equation will yield a linear C-V relation.The mechanism behind this linear relation is the fluid analogue of the classical kinetic result (Knight, 1973).Most of the energy gained by the electrons falling through the field-aligned potential drop is converted into an increased kinetic temperature.The resulting temperature gradients are large, and it is essential to include them in the fluid momentum equation.Since the approximants T ⊥ and T z only depend on parameters that are available within a fluid model, they can be used as nonlocal equations of state to close the set of fluid equations.
Earlier calculations of fluid quantities (e.g.Knight, 1973;Chiu and Schulz, 1978;Fridman and Lemaire, 1980;Janhunen, 1999) have all been based on assumptions, such as Eqs.13) or ( 14), that allow the limits of integration at z to be determined from the local turning point line μ(z, H ). The new method we have devised takes the shape of the potential fully into account.This is important, since condition ( 14) is violated when the potential increases faster than B, which is needed to maintain quasi-neutrality and is also suggested by observations.The density and temperatures are sensitive to the shape of the potential, even if the F-L condition is satisfied so that the current is independent of the potential shape.
Models of quasi-neutral parallel electric fields, including a kinetic, self-consistent description of the ions, have been developed by Chiu and Schulz (1978), Stern (1981), and Ergun et al. (2000).While this approach is formally elegant, it is not obviously superior in practice.Due to the slow motion of the ions, such a model can only be safely applied after the potential has remained stationary for several minutes.During this time, quasi-neutrality must be maintained by the electrons, and a model with a given, possibly slowly time dependent ion density profile is more useful.
In this study we have sometimes sacrificed realism in details to maintain simplicity and a close connection with the classical kinetic model.As seen from Fig. 5, the electron densities in our model depend to some extent on the magnitude of the potential drop.It would be possible to adjust the electron density further by allowing the boundary distribution functions or the shape of the potential to depend on φ.While it is likely that such dependences exist, we prefer to keep the model simple by ignoring them.The consequences of small errors in the density are not too severe for our pressure and temperature estimates.We have also made calculations with the exponent p=1.5 and p=2.5 in the model potential (23), and the calculated temperature profiles agree to within about 15% with the results for p=2.
There are some conceptual similarities between this study and the work of Janhunen (1999), who also used the density and pressure calculated from a kinetic model to obtain an equation of state.However, there are also a number of important differences.Janhunen does not require that the electron density should be independent of the potential drop, and he considers only a scalar pressure.Janhunen also assumes a local equation of state of the form P =C>n γ , which makes it hard to obtain a good fit to the pressure even in the case φ=0, where he derives simple formulas for the altitude variation of P and n.He finds that γ =3 should produce a linear C-V relation, but also makes several reservations in his conclusions.Still, some of the important questions answered in this study were first asked by Janhunen (1999).
The temperature approximants T ⊥ and T z apply to quasistationary situations, when the electrons have reached an equilibrium with the potential.Considering that the potential variations are concentrated to the lower end of the field line and that the velocity of a 1 keV electron is about 3 R E /s, we still expect our approximations to be better than a conventional equation of state (e.g.isothermal) when the potential varies on time scales slower than a few seconds.The utility of our particular formulas for the temperatures may also be limited by the assumptions that the source distributions are isotropic and that the potential variation is proportional to B 2 .We have also neglected backscattered and secondary electrons, as well as magnetospheric electrons that may become trapped below the potential drop during its buildup (Eliasson et al., 1979).Such electrons may contribute significantly to the density and pressure in the lower part of the acceleration region.Although specific problems may require adjustments of the model, it seems clear from our results that nonlocal equations of state are needed in any fluid model of active auroral flux tubes.

Conclusions
Non-local equations of state are required to describe the temperature variations in an auroral flux tube with upward current.Uniformly valid approximations for the temperature can be derived from kinetic theory, and expressed in terms of fluid variables.When supplemented with the nonlocal equations of state the electron momentum equation will determine the local parallel electric field from the current density.Integration of this parallel electric field will yield a potential drop φ that is linearly related to the current in the flux tube.These equations of state will allow fluid simulations of magnetosphere-ionosphere coupling that include a realistic description of the parallel electric field.P ⊥ (B I , ) = 0.50 + 1.5 + N I T I N G T G , where the values at =0 are taken from P 0 ⊥ in (A4).The other numerical values in these expression are chosen to make the Padé approximant fit the calculated pressure profiles as well as possible.The Padé approximant is given by where the coefficients Q i are The coefficient R 0 =1, while R 1 , R 2 , and R 3 are obtained as the solution of Eq. (B3).
For the parallel pressure the Padé approximant is given by where the coefficients Q i are The coefficient R 0 =1, while R 1 , R 2 , R 3 , and R 4 are obtained as the solution of Eq. (B5). (B5)

Fig. 1 .
Fig. 1.Illustration to the first steps in the calculation of the turning point boundary µ tpb .The shaded area marks the region in the µ-H plane accessible to downgoing electrons.

Fig. 2 .
Fig. 2. The area accessible to downgoing electrons is the part in the µ-H plane above the turning point boundary µ tpb if the shape of the electrostatic potential is taken fully into account.The difference between this boundary and the local turning point line is shaded.

Fig. 3 .
Fig. 3. Solid, bold lines show the integration boundaries in the velocity space.The dashed, bold line is the integration boundary corresponding to the local value μ(z, H ) valid for potentials (23) with p≤1.The thinner lines show some contours of the distribution function corresponding to the values 1, 10 −2 and 10 −4 when normalized with the maximum value of the distribution function.The figure applies to z=z a .

Fig. 4 .
Fig. 4. Electron density profiles for φ=10 kV and different p in the model (23), together with the ion density.The point z a marks the bottom of the acceleration region.

Fig. 5 .
Fig. 5. Dependence of the electron density profile on the potential drop φ for p=2.

Fig. 6 .
Fig.6.Current-Voltage relation computed with the described method compared with the linear approximation from Eq. (24).

Fig. 7 .
Fig. 7. Parallel temperature profiles compared with the approximations T z (dashed line) for different φ.

Fig. 9 .
Fig.9.Comparison between the model potential and the potential derived from the approximants to the temperatures at different current densities.