O + and H + ion heat fluxes at high altitudes and high latitudes

Higher order moments, e.g., perpendicular and parallel heat fluxes, are related to non-Maxwellian plasma distributions. Such distributions are common when the plasma environment is not collision dominated. In the polar wind and auroral regions, the ion outflow is collisionless at altitudes above about 1.2 RE geocentric. In these regions wave–particle interaction is the primary acceleration mechanism of outflowing ionospheric origin ions. We present the altitude profiles of actual and “thermalized” heat fluxes for major ion species in the collisionless region by using the Barghouthi model. By comparing the actual and “thermalized” heat fluxes, we can see whether the heat flux corresponds to a small perturbation of an approximately biMaxwellian distribution (actual heat flux is small compared to “thermalized” heat flux), or whether it represents a significant deviation (actual heat flux equal or larger than “thermalized” heat flux). The model takes into account ion heating due to wave–particle interactions as well as the effects of gravity, ambipolar electric field, and divergence of geomagnetic field lines. In the discussion of the ion heat fluxes, we find that (1) the role of the ions located in the energetic tail of the ion velocity distribution function is very significant and has to be taken into consideration when modeling the ion heat flux at high altitudes and high latitudes; (2) at times the parallel and perpendicular heat fluxes have different signs at the same altitude. This indicates that the parallel and perpendicular parts of the ion energy are being transported in opposite directions. This behavior is the result of many competing processes; (3) we identify altitude regions where the actual heat flux is small as compared to the “thermalized” heat flux. In such regions we expect transport equation solutions based on perturbations of bi-Maxwellian distributions to be applicable. This is true for large altitude intervals for protons, but only the lowest altitudes for oxygen.


Introduction
Non-Maxwellian ion velocity distributions have been observed at different high altitudes and high latitudes (e.g., Slapak et al., 2011;Waara et al., 2010;Huddleston et al., 2000;Winningham and Burch, 1984).The features of these velocity distributions are related to higher order moments like heat flux.For example, Biddle et al. (1985) have observed ion heat flux in light ion polar wind.Analyses of satellite measurements have indicated that polar wind can contain suprathermal components of both light and heavy ions (e.g., Moore et al., 1985Moore et al., , 1986)).These observations are characteristic of any plasma flow where collisions are not dominated.Different theoretical studies have discussed the ion outflows in the polar wind and auroral regions either by using kinetic models (e.g., Lemaire and Scherer, 1973, and for more details see the review by Tam et al., 2007), transport theory approach (e.g., Schunk and Watkins, 1982;Demars and Schunk, 1989, 1992, 1994;Khazanov et al., 1984;Ganguli et al., 1987, Ganguli andPalmadesso, 1987) or the Monte Carlo method (e.g., Barakat and Schunk, 1983;Barghouthi et al., 1993Barghouthi et al., , 2001;;Barakat et al., 1995).These models have investigated the ion outflows either in the collision-dominated region or collisionless regions or in both regions.These studies did not include the effect of wave-particle interactions.Other studies have considered the collisionless plasma outflows in the polar wind and auroral regions and have included the effect of wave-particle interactions.These studies have replaced the Fokker-Planck collision term in the Boltzmann equation by the effect of wave-particle interactions (e.g., Published by Copernicus Publications on behalf of the European Geosciences Union.I. A. Barghouthi et al.: O + and H + ion heat fluxes Retterer et al., 1987Retterer et al., , 1994;;Crew et al., 1990;Barghouthi et al., 1998Barghouthi et al., , 2007Barghouthi et al., , 2008; Barghouthi and Atout, 2006;Barghouthi, 1997Barghouthi, , 2008;;Bouhram et al., 2003aBouhram et al., , b, 2004)).In these studies the ion motion was followed along geomagnetic field lines under the effects of external forces (gravity and polarization electric field) and the interactions between the ion and the electromagnetic waves observed in those regions.They obtained altitude profiles for ions density, drift velocity, and parallel and perpendicular temperatures and the ion velocity distribution functions.These and other similar studies have not discussed the higher order moments of the ion velocity distribution such as ion heat flux and viscous stress.However, Biddle et al. (1985) observed ion heat flux in light polar wind and found that collisionless effects are stronger than predicted and dominate the heat transport throughout the entire transition from supersonic to subsonic flow.Their observations demonstrate a capability to observe transport effects in plasmas as higher order moments of the ion velocity distribution function; also, they provide evidence that there is a significant difference between the flow velocity and temperature of the ion velocity distribution core and the energetic tail.
Two published models by Brown et al. (1995) and Wu et al. (1999) discussed ion outflow.These models are collisional at lower altitudes and collisionless at higher altitudes, which extend up to geocentric distances of 3 R E .The Barghouthi model (Barghouthi, 2008) is similar to the collisionless region of their models, but extend the flux tube to 10 R E .The Barghouthi model used altitude-and velocity-dependent diffusion coefficient, which are important for the much larger altitude interval spanned by this model as shown in Barghouthi (1997Barghouthi ( , 2008)).
There is a series of papers about a model which has sought to include wave-particle interactions, as well as the effect of field line divergence, within a multi-moment fluid context (Jasperse et al., 2006a(Jasperse et al., , b, 2010a, b), b).The model includes wave-particle interactions through a Fokker-Planck formalism, as do most of the above references.In particular, Jasperse et al. (2010a) give an approximate method for evaluating wave-particle interactions based on satellite measurements for a particular wave mode.There is also at least one published paper of a kinetic model which includes both wave-particle interactions and diverging magnetic field (Jasperse, 1998).Although the Jasperse et al. work focuses on calculating the large-scale potential on auroral field lines, they also calculate how electron and ion moments evolve with altitude.Note that parallel electric fields, both quasistatic and Alfvénic, should significantly affect the dynamics in the auroral zone.
Currently, we and Barghouthi et al. (2011) are searching through a significant amount of observational literature on ion outflow at high altitude and high latitude (e.g., Yau and Andre, 1997;Andre and Yau, 1997;Paschmann et al., 2002;and Moore and Horwitz, 2007), and according to our knowledge, there is no observations related to O + and H + ion heat fluxes in the altitude range (1.7-10 R E ), i.e., the range of our of study.
The objectives of this paper are the following: (1) to investigate the general trends of the H + and O + ion outflows in the polar wind and auroral zone by using Barghouthi model,(2) to explain the ion heat flux in terms of the ion velocity distribution function and to single out the contributions of the ions in the energetic tail of the velocity distribution, (3) to provide typical altitude profiles for H + and O + ion heat fluxes in the collisionless polar wind region and auroral zone -these numerical values are important for fluid model closure, since most fluid models assume a particular form of the heat flux to produce a solvable system of equations, and (4) to provide a recommendation on the usefulness of the 16-moment equations for discussing H + and O + ion outflows in the polar wind and auroral regions.
This paper is organized as follows: following a description of the theoretical formulations (Sect.2), we present H + and O + ion velocity distributions (Sect.3), compute and discuss H + and O + ion heat fluxes (Sect.4), compare the actual H + and O + ion heat fluxes and H + and O + ion "thermalized" heat fluxes (Sect.5), discuss the applicability of our results (Sect.6), and draw some general conclusions based on our results (Sect.7).

Theoretical formulations
The ion heat flux at high altitudes and high latitudes under the effects of gravity, polarization electric field, diverging geomagnetic field lines, and the effect of wave-particle interactions (WPIs) can be obtained by integrating the product (v s − u s ) 2 (v s − u s )f s over all velocities and then dividing by the ion density n s (Schunk and Nagy, 2000).It is given mathematically by the following expression: where c s is the thermal velocity, and it is equal to v s − u s , where u s is the drift velocity, and v s is the ion velocity, f s is the ion velocity distribution and is obtained by solving the following Boltzmann equation: where g is the acceleration of gravity, E is the polarization electrostatic field, B is the geomagnetic field, e s and m s are the charge and the mass of the ion (i.e., either H + or O + ) respectively, c is the speed of light, ∂/∂t is the time derivatives, ∇ is the coordinate space gradient, and ∇ vs is the velocity space gradient.The right-hand side of the Boltzmann equation δf s /δt represents the rate of change of f s (v s , r s , t) in a given region of phase space (v s , r s ) as a result of collisions, i.e., the Fokker-Planck collision term.Because we are in a collisionless region and the effect of wave-particle interactions is significant, we have replaced this Fokker-Planck collision term by the following diffusion equation as shown in Retterer et al. (1987): where (D ⊥ ) is the quasi-linear velocity diffusion coefficient.
The influence of wave-particle interactions on the ion during t is taken into consideration by incrementing the ion's perpendicular velocity by a random increment v ⊥ such that where t is the time interval chosen randomly, and D ⊥ (r, v ⊥ ) is the quasi-linear velocity diffusion rate perpendicular to the geomagnetic field B, and is given by the following expression (Barghouthi, 2008): where i is the ion gyrofrequency, and k ⊥ is perpendicular wave number and related to the characteristic perpendicular wavelength of the electromagnetic turbulence λ ⊥ .
In the polar wind region, D ⊥ (r) is given by Barghouthi et al. (1998) as follows: And in the auroral region D ⊥ (r) is given by Barghouthi (1997) as follows: These altitude-dependent diffusion coefficients given in Eqs. ( 6) and ( 7) are derived empirically from fits to DE-1 data.Due to the shortage of data for electric wave activity at high altitudes, we have extrapolated the low-altitude diffusion coefficients to high altitude.Recently, Slapak et al. (2011) and Waara et al. (2011) have provided updated diffusion coefficients for the high-altitude polar cap.These show that the altitude dependence of the electric field wave activity as well as the absolute level is lower than implied by our extrapolation.On the other hand, the fraction of the wave activity which is efficient in heating the ions would have to be higher as compared to our study in order to explain their observations.We have therefore chosen to use our old values in this study.
Different models have been used to solve the Boltzmann equation in the collisionless region, such as kinetic, semikinetic, hydromagnetic, generalized transport, and Monte Carlo models.In general, these models were used to find the ion velocity distributions, density, drift velocity, and parallel and perpendicular temperatures.According to our knowledge, none of the related studies (i.e., the studies of the ion outflows at high altitudes and high latitudes) have investigated the effects of gravity, polarization electric field, diverging geomagnetic field lines, and wave-particle interactions on the ion heat flux in collisionless polar wind or auroral regions.
According to Schunk and Nagy (2000), the heat flux has been given in terms of the parallel and perpendicular heat fluxes as follows: The magnitudes of the total, parallel and perpendicular heat fluxes are given by q s = 1 2 (q s + 2q s⊥ ), q s = q s • b = n s m s (v s − u s ) 3 , and with b defined as the unit vector in the direction of the geomagnetic field B.
The ion motion under the effect of gravity, polarization electric field, geomagnetic field, and wave-particle interactions has been described in detail by Barghouthi (2008).In essence, we will inject a test ion (H + or O + ) at the lower boundary and follow the motion of the test ion under the effects of the above external forces and wave-particle interactions.We keep updating the velocities and the position of the test ion at the end of each randomly chosen period of time.At different altitudes in the simulation tube, we accumulate data about the test ions, and consequently the analysis of this data at the given altitude produces the ion velocity distribution and the velocity moments.Here we are interested in the velocity moments that give the above total, parallel and perpendicular heat fluxes.The Monte Carlo code, i.e., the Barghouthi model (Barghouthi, 2008), has been used to compute the actual heat fluxes in the polar wind and auroral regions and for H + and O + ions.

I. A. Barghouthi et al.: O + and H + ion heat fluxes
The boundary conditions selected for the polar wind region are similar to those of Barghouthi et al. (1998).At the lower boundary (i.e., 1.7 R E ), we set the O + ion drift velocity at 0 cm s −1 , the oxygen ion density at 100 cm −3 , and the O + ion temperature at 3000 K.For H + ions, we set the H + ion drift velocity at 11 km s −1 , the hydrogen ion density at 200 cm −3 , and the H + ion temperature at 3000 K.The electron temperature was kept constant at 1000 K along the entire simulation tube (1.7-10 R E ).We also assumed the velocity distribution to be a drifting Maxwellian for H + ions and the up-going half of non-drifting Maxwellian for O + ions at the lower boundary.
The boundary conditions selected for the auroral region are similar to those of Barghouthi (2008).At the lower boundary (1.2 R E ), we set the H + ion drift velocity at 16 km s −1 , the H + ion density at 100 cm −3 , and the H + ion temperature at the lower boundary at 0.2 eV (2320 K).We set the O + ion drift velocity at 0 km s −1 , the oxygen ion density at 5000 cm −3 , and the O + ion temperature at 0.2 eV (2320 K).The electron temperature is kept constant at 1000 K along the entire simulation tube (1.2-10 R E ).We also assumed the velocity distribution to be a drifting Maxwellian for H + ions and the up-going half of non-drifting Maxwellian for O + ions at the lower boundary.The geomagnetic field B was taken to be proportional to r −3 , where r is the geocentric distance.
The potential energy due to the body forces (i.e., gravitational force and electrostatic polarization field) is given by Barakat and Schunk (1983) as follows: where k is Boltzmann's constant; T e is the electron temperature which is kept constant in the simulation tube; and n e and n e0 are the electron densities at r and r 0 (i.e., the lower boundary), respectively.These can be calculated from the quasi-neutrality condition [n e = n(O + ) + n(H + )].G is the gravitational constant, M E the mass of the earth, and m the ion's mass (i.e., H + or O + ).The altitude profiles of the potential energy due to the body forces given by Eq. ( 11) are obtained self-consistently: running the Barghouthi model for the above conditions produces ion densities at different altitudes, and these densities will be inserted into Eq.( 11) to calculate the potential energy.We keep repeating this till the results converge.Different altitude profiles for ion potential energy in both regions are given in Barghouthi et al. (2011) and Barghouthi (2008).

H + and O + ion velocity distributions
The main objective of this paper is not to present and discuss, in detail, the H + and O + ion velocity distributions in the polar wind and auroral regions.It has been discussed by Barghouthi et al. (2011) and Barghouthi (2008).However, one of the main objectives of this study is to discuss H + and O + ion heat fluxes in polar wind and auroral regions and to investigate the relationship between the shape of the ion velocity distribution function and the ion heat flux.According to Eqs. ( 8)-( 10), the H + and O + ion heat fluxes in the polar wind and auroral regions are determined by the ion velocity distribution functions.Therefore, Fig. 1 presents contour plots for H + and O + ion velocity distributions in both regions and at different altitudes.Briefly, in the polar wind region f (H + ) moves from Maxwellian at lower altitude (i.e., it is consistent with the initial conditions) to bi-Maxwellian at 2 R E ; there is a temperature anisotropy.The width of the distribution function in the perpendicular direction is greater than the width in the parallel direction.This means that the ion perpendicular temperature is greater than the ion parallel temperature (see Fig. 2c, and d, solid lines).At altitude 5 R E this anisotropy has been reversed.At around 7 R E , the anisotropy is inverted again and the lower half of the ion velocity distribution is larger than the upper half.At higher altitude the velocity distribution function shows conic features due to mirror force, and at very high altitudes (not shown here) the toroidal features are well established.
The O + ion velocity distribution function in the polar wind moves directly from a Maxwellian distribution at 2 R E to a conic distribution at 4 R E , and then to a toroid at 6 R E and finally to a well-pronounced toroid at higher altitudes.
In the auroral region the effect of wave-particle interactions is stronger than in the polar wind region.Therefore, the role of wave-particle interactions in heating the ions in the perpendicular direction produces an increase in the ions perpendicular temperature (energy) as shown in Fig. 2g.Part of this gained energy will be converted to the parallel direction by the mirror force in combination with energy conservation, resulting in an increase of the parallel temperature (energy).The H + ion velocity distribution in the auroral region displays a bi-Maxwellian shape at 2 R E , conic at 3 R E , and toroids at 8 R E and above.The toroidal features are well established and appear at low altitudes; this is due to the effect of wave-particle interactions.O + ions in the auroral region move from Maxwellian at 1.2 R E , to conic at 2 R E , and to well-pronounced toroid at 4 R E and altitudes above.It is obvious that the influence of wave-particle interactions on O + ion velocity distributions is stronger than on H + ion velocity distributions, and in the auroral region stronger than in the polar wind region.

Heat flux simulations
In this section, we present the perpendicular, parallel, and total actual heat fluxes (solid lines) for H + and O + ions in the polar wind and auroral regions (Figs. 3,and 4,respectively).In Sect.5, we present the corresponding "thermalized" heat fluxes (Figs. 3 and 4, dashed lines).

Polar Wind Region
Aurora Region H 1/2 , where T is the ion temperature at the given geocentric altitude.The contour levels decrease successively by a factor e 1/2 from the maximum.
We will discuss the behavior of the actual heat fluxes in light of the following arguments: 1.The effect of wave-particle interactions.This effect contributes to heat (energize) the ions in the perpendicular direction, and consequently part of this energy transfers to the parallel direction due to mirror force.
2. The effect of perpendicular adiabatic cooling.This effect contributes, through the mirror force, to transfer energy from perpendicular direction to parallel direction, as ions drift upward along geomagnetic field lines, in order to keep the first adiabatic invariant and energy constant.
3. The effect of ion potential energy.This potential energy is due to the gravity and polarization electric field.
Previously, Barghouthi and Barakat (1995), Barghouthi et al. (1998), and Barghouthi (1997Barghouthi ( , 2008) ) showed that O + ions are heated more than H + ions in the polar wind and auroral regions as shown in Fig. 2. Similarly, Chang et al. (1986) gave an explicit mass dependence for ion heating due to a power law spectrum.The heating process in the auroral region is higher than that in the polar wind region, as the perpendicular diffusion coefficient, which is directly proportional to the power spectral density for the electric field of the waves, in the auroral region is much higher than that in the polar wind region.O + ions at the lower boundary are gravitationally bound (i.e., they move in a monotonically increasing potential, the potential that is due to the gravity and polarization electrostatic field (Eq.11) and need more energy to overcome the potential barrier).These ions are energized due to wave-particle interactions and then move upward along geomagnetic field lines, due to the mirror force.Ions located in the energetic tail may overcome the potential barrier.Ions that do not overcome the potential barrier will be reflected downward and get more energy due to wave-particle interactions and then reflected upward due to the mirror force.They will keep moving up and down till they can overcome the potential barrier and finally escape from the top of the simulation tube.This is the pressure cooker effect (Barakat and Barghouthi, 1994;Gorney et al., 1985), and this explains why O + ions are heated more than H + ions at low altitude.

H + and O + heat fluxes in the polar wind
Based on the definitions of the vector and magnitude of the parallel, perpendicular and total heat fluxes, given in Sect.2, and before presenting their altitude profiles, it is worthwhile to mention the following general characteristics of the ion heat flux: 1. Total heat flux represents the rate of the total energy transfer per unit area per unit time along geomagnetic field lines.Positive heat flux means that the energy is moving upward along geomagnetic field lines.That is, energy moves from low altitude to higher altitude.In other words, the ion velocity distribution function has an upward energetic tail.
2. Negative heat flux means that the energy is moving downward along geomagnetic field lines.That is, energy is moving from higher altitude to lower altitude.The ion velocity distribution has a lower tail.
3. When the heat flux decreases (increases), this means that the rate of energy transfer along geomagnetic field lines decreases (increases).
4. When the parallel heat flux increases (decreases) and is positive (negative), this indicates that the rate of the parallel part of the ion energy transported increases (decreases) in the upward (downward) direction along geomagnetic field lines.
5. When the perpendicular heat flux increases (decreases) and is positive (negative), this indicates that the rate of the perpendicular part of the ion energy transported increases (decreases) in the upward (downward) direction along geomagnetic field lines.
It is important to note that, in Figs. 3 and 4, the graphs are given in logarithmic scale for the heat flux.The solid line shows the magnitude of the actual heat flux, with red color indicating positive values and blue color indicating negative values.
Figure 3 (solid lines) presents the magnitudes of the perpendicular (top panels), parallel (middle panels), and total (lower panels) actual heat fluxes in the polar wind for H + ions (left panels) and O + ions (right panels).
The behavior of the actual H + perpendicular heat flux (Fig. 3a) is controlled by the competition between the effects of wave-particle interactions and perpendicular adiabatic cooling.At altitudes below 7.34 R E it decreases, becomes negative at 6 R E , and reaches its minimum value of about −1×10 −11 erg cm −2 s −1 at 7.34 R E .At altitudes above 7.34 R E it starts to increase, and becomes positive above 8 R E .To understand the behavior (decreases, increases, positive, and negative) of the actual H + perpendicular heat flux, we need to consider the definition of that perpendicular heat flux, i.e., Eq. (10).At low altitudes, perpendicular adiabatic cooling is significant and dominates the effect of waveparticle interactions in the polar wind region.Therefore, as altitude increases, perpendicular velocity decreases, and this decreases the perpendicular temperature (as shown in Figs. 1  and 2c).In other words, the width of the ion velocity distribution is decreasing along the perpendicular direction.This decreases the actual H + perpendicular heat flux.At altitudes above 7.34 R E , the effect of wave-particle interactions dominates the effect of perpendicular adiabatic cooling.Thus perpendicular velocity starts to increase, and this increases corresponding temperature and the actual H + perpendicular heat flux.As altitude increases, energy transfer from perpendicular direction to parallel direction enhances the bulk drift velocity and makes it larger than H + parallel velocities for H + ions located in the lower part (tail) of the distribution function.This produces a negative sign in Eq. ( 10), which means that the contributions of the ions in the lower tail of the distribution are very important.In other words, when the bulk drift velocity of the whole H + population becomes larger and the number of ions in the lower tail becomes large, this will enhance the negative values of the heat flux.This explains the behavior of the actual H + perpendicular heat flux in the altitude range 6 to 8 R E .These findings are consistent with the H + ion velocity distributions (Fig. 1, left panel).There is almost symmetry between upper and lower parts of the distribution functions at about 6 and 8 R E .This symmetry cancels the contribution that is due to the difference between ion parallel velocity and the bulk drift velocity.Therefore, we have zero heat flux at these altitudes.However at around 7 R E the contribution of the lower part is bigger.In other words, we have more ions in the lower part of the distribution function.Therefore, we have negative heat flux.The actual H + parallel heat flux as shown in Fig. 3b (solid line) decreases with altitude, and slightly increases at very high altitudes.As mentioned before, the effect of perpendicular adiabatic cooling is, through conservation of energy, to enhance the bulk drift velocity of the H + ions in the upward direction.The increase in the parallel energy at the expense of the parallel thermal velocity is known as the parallel adiabatic cooling (Barakat and Lemaire, 1990).The actual H + parallel heat flux depends on the cubic difference between ion parallel velocity and bulk drift velocity.When the drift velocity increases, the difference between the parallel velocity of particles in the energetic tail and the drift velocity decreases.Thus the temperature related to the parallel velocity decreases, and the parallel heat flux will be decreased.At very low altitudes the difference between the parallel velocity for ions located in the forward energetic tail and the bulk drift velocity is large.However as the ions drift upward this difference turns to decrease to reach a minimum value at about 9 R E .At very high altitudes it increases because the difference between drift velocity and the thermal parallel velocity starts to increase due to a combination of energy input to the ions due to wave-particle interactions and transfer from perpendicular direction to the parallel direction due to the mirror force.
It is important to outline that the contributions of the particles located in the tail of the distribution are very significant to the magnitude of the heat flux because the difference between the bulk drift velocity and the parallel velocities of these ions is very large, and consequently the heat flux becomes very large.Therefore it is crucial to take into consideration the contributions of the particles in the tail of the velocity distribution and not to concentrate on the contribution of the bulk ions.The difference between parallel velocities of the bulk ions and the drift velocity is small, and the net heat flux is minimal.For example the heat flux for a Maxwellian distribution is zero.
The lower panel of Fig. 3 shows the altitude profile of the actual total H + heat flux.Below about 8 R E the total heat flux decreases because both parallel and perpendicular heat fluxes are decreasing.Above 8 R E it increases because perpendicular heat flux was increasing very rapidly, and the parallel heat flux decreases very slowly.In general, the total heat flux is always positive; this means that the energy is moving in the upward direction along the geomagnetic field lines.Figure 3 (right panels, solid lines) presents the behavior of O + actual perpendicular (panel d), parallel (panel e), and total (panel f) heat fluxes.The differences between the behavior of O + ions and H + ions are due to the effect of wave-particle interactions which dominates the effect of perpendicular adiabatic cooling at low altitudes and to the potential energy barrier which is monotonically decreasing for H + ions and monotonically increasing for O + ions (Barghouthi et al., 2011).O + actual perpendicular heat flux (Fig. 3d) decreases at very low altitudes, at 2.5 R E reaches its minimum negative value of −1.11 × 10 −11 erg cm −2 s −1 , and then starts to increase and at about 3 R E becomes zero.In the range 3 to 5.2 R E it increases very rapidly, and at altitudes higher than 5.2 R E it decreases very slowly.At very low altitudes, the perpendicular temperature is almost constant due to the competition between perpendicular adiabatic cooling and ion heating due to wave-particle interactions; therefore the contribution from the square term in Eq. ( 10) is constant.The main contribution to the actual perpendicular heat flux, at low altitudes, is controlled by two contributions: (1) the ions in the tail of the distribution which overcome the potential barrier produce positive value, and (2) the ions in the bulk of the distribution and reflected downward due to the potential barrier produce negative value.It is clear that the role of O + potential barrier is to create a downward tail; this is very obvious in the shape of the ion velocity distribution as shown in Fig. heat flux at low altitudes by the following argument: as altitude increases from 1.2 to 3 R E , the number of O + ions in the lower part of the distribution increases, and this increases the negative contribution to the perpendicular heat flux.On the contrary, the number of ions that overcome the potential barrier and contribute positively to the perpendicular heat flux decreases.The actual perpendicular heat flux increases very quickly in the range 2.8 to 5.2 R E because wave-particle interactions turns out to be very effective; consequently most of O + ions overcome the potential barrier and move upwardthis will contribute positively to the perpendicular heat flux.However at altitudes above 5.2 R E it decreases slowly because of the effect of finite gyroradius (Bouhram et al., 2004;Barghouthi and Atout, 2006) which makes the ion heating process become saturated.These results are consistent with the contour plots of the O + ion velocity distribution in the polar wind (Fig. 1).It is obvious that at 3 R E the lower tail is broader than the upper part, and the distribution function saturated above 6 R E .O + actual parallel heat flux (Fig. 3e, solid line) is always positive, decreases at altitudes below 4 R E , increases very rapidly in the range 4 to 5.5 R E , and above 5.5 R E it decreases very slowly.The behavior of O + actual parallel heat flux can be explained as follows: (1) at low altitudes, the ions that are located in the tail of the velocity distribution need less additional energy to overcome the O + potential barrier.Now it looks like we have a group of O + ions with high parallel velocity which are above the O + potential barrier, and the rest are coupled in the bulk of the ion velocity distribution and are below the O + potential barrier.According to Eq. (9) these highly, but few, energized ions contribute more and positively to the parallel heat flux, while ions in the bulk contribute negatively to the heat flux.The combination of these two contributions produces the altitude profile at low altitude.(2) In the altitude range 4 to 5.5 R E , the majority of O + ions drift upward and have overcome the O + potential barrier; this means that the contribution of the difference between parallel and drift velocities to the parallel heat flux (Eq.9) is increasing.In other words the waveparticle interactions heat the ions.The faster ions overcome the O + potential barrier and produce a large difference between its parallel velocities and the drift velocity of the whole distribution.This process keeps going: faster ions overcome the barriers before the bulk ions.Then these bulk ions keep moving upward and reflected downward, become energized due to wave-particle interactions and then drift upward till they overcome the O + potential barrier.This is the pressure cooker effect (Barakat and Barghouthi, 1994).(3) At 5.5 R E and above, the heating process due to wave-particle interactions is saturated.This is because at a specific altitude and when the perpendicular wavelength of the electromagnetic turbulence, which is responsible for ion heating, is comparable to or greater than the ion gyroradius, the heating process saturates.The O + velocity distribution function (Fig. 1) does not change, and the changes in the ion outflow characteristics are very minimal (Bouhram et al., 2004;Barghouthi, 2008).
The lower panel of Fig. 3 gives the total actual O + heat flux (panel f); it is positive and increases, very rapidly, in the altitude range from 1.7 to 5 R E due to the above mentioned behaviors of parallel and perpendicular heat fluxes.Above 5 R E it remains almost constant and positive.The total actual O + heat flux is always positive.This means that the energy has been transferred from low altitudes to higher altitudes along geomagnetic field lines.Above 5 R E the rate of energy transfer is almost constant and we have almost a saturated process: the shape of the ion velocity distribution does not change.The shape is expanded equally, but slowly, in all directions.As mentioned before, the plots of the ion velocity distribution are normalized with respect to the total temperature at that altitude and shifted backward along the parallel direction by the amount of the drift velocity at that altitude; otherwise, it goes out of the page.

H + and O + heat fluxes in the auroral region
Before discussing the ion heat flux in the auroral region, it is helpful to emphasize that the ion heating process in the auroral region is more efficient than in the polar wind region.The strength of the electric spectral densities of the electromagnetic waves, and consequently the diffusion coefficients are much higher in the auroral region than in the polar wind region.Also, O + potential energy due to the polarization electric field and gravity is much higher in the auroral region than in the polar wind region.This creates a barrier in the upward direction, and thus O + ions need more energy to overcome this barrier.Finally, O + ions with enough energy will be able to overcome this barrier and therefore appear with a higher energy (temperature) than the corresponding ions in the polar wind region.In light of these arguments, we will discuss the ion heat fluxes in the auroral region.Figure 4 (solid lines) presents the magnitudes of the perpendicular (top panels), parallel (middle panels), and total (lower panel) actual heat fluxes in the auroral region for H + ions (left panels) and O + ions (right panels).H + actual perpendicular heat flux (Fig. 4a) increases at very low altitude due to the assumed boundary conditions.At about 1.4 R E it becomes negative and reaches its minimum value of −1.88 × 10 −8 erg cm −2 s −1 at 1.8 R E .Between 2 and 4 R E it increases very rapidly and reaches its maximum positive value, because the effect of wave-particle interactions becomes very significant.Consequently, this increases the perpendicular velocity, and this increases the corresponding perpendicular heat flux.Above 4 R E it decreases very slowly because the ion heating turned out to be saturated due to the effect of finite gyroradius (mentioned above).However, the effect of perpendicular adiabatic cooling comes into play, and part of the energy in the perpendicular direction transfers into the parallel direction to keep the first adiabatic invariant and energy constant.This decreases the perpendicular velocity and consequently the associated perpendicular heat flux.
The actual parallel heat flux of H + ions is almost constant in the altitude range 1.2 to 2.5 R E , increasing very rapidly in the altitude range from 2.5 to 4 R E , reaches its maximum positive value at 4 R E , and then reverses and decreases slowly above 4 R E .It is important to know that, in the case of H + ions, the potential energy is always negative and decreasing as altitude increases.Therefore, it will not act as a potential barrier (Barghouthi, 2008).In the altitude range 1.2 to 2.5 R E , the effects of wave-particle interaction and perpendicular adiabatic cooling are competing with each other, this makes the difference between parallel and drift velocities to be constant.Consequently, according to Eq. ( 9) the parallel heat flux remains constant.In the altitude range 2.5 to 4 R E , the heating process due to wave-particle interactions increases H + ions perpendicular temperature (energy), and part of this gained energy will be converted into the parallel direction in order to keep the first adiabatic invariant and energy constant.This results in an increase in the ion parallel velocity and drift velocity.The H + ions located in the tail of the velocity distribution are already highly energized, and with little heating due to wave-particle interactions they will move faster than the rest of the bulk ions.Therefore, their contribution, even if there were few ions at the beginning, will be high because it is given in terms of the cubic difference between parallel velocity and drift velocity as shown in Eq. ( 9).This explains the rapid increasing behavior of the actual parallel heat flux in that range.The actual H + parallel heat flux decreases slowly above 4 R E .This is due to the effect of parallel adiabatic cooling (mentioned before).The total actual H + heat flux is shown in the lower panel of Fig. 4 (panel c, solid line).The behavior of the total heat flux is controlled by the contributions of both parallel and perpendicular heat fluxes and because, almost, both of them are increasing at altitudes below 4 R E .Therefore, it increases at low altitudes.Above 4 R E it decreases, very slowly, similar to the behaviors of parallel and perpendicular heat fluxes.
Figure 4 (right panels, solid lines) presents actual perpendicular (top panel), actual parallel (middle panel), and actual total (lower panel) heat fluxes for O + ions in the auroral region.In the case of O + ions, the potential energy due to gravity and the polarization electric field is always positive and monotonically increasing (Barghouthi, 2008).Therefore, when O + ions drift upward due to the mirror force, they will face that potential barrier.Ions with high parallel velocity (i.e., those located in the upward tail of the ion velocity distribution) will be able to overcome that potential barrier and escape to higher altitudes and leave the bulk of the ions.
At altitudes below 2 R E , both actual O + perpendicular (Fig. 4d) and parallel (Fig. 4e) heat fluxes are increasing rapidly and positively.This is due to the energetic O + ions that are located in the tail of the velocity distribution.These ions are able to overcome the O + potential barrier.These highly energized escaping ions contribute much to the parallel heat flux in terms of the cubic difference between parallel and drift velocities.They also contribute to the perpendicular heat flux but with fewer amounts, i.e., in terms of the difference between parallel and drift velocities.Also, for the case of O + ions in the auroral region the effect of waveparticle interaction dominates, contributing to increasing perpendicular temperature (Fig. 2g, dashed line) and drift velocity (Fig. 2f).These contribute to increase parallel and perpendicular heat fluxes.Both fluxes are decreasing, slowly, above 2 R E , because the majority of O + ions are drifting upward and the difference between parallel and drift velocities is decreasing slowly.Consequently, the contribution of this difference to Eqs. ( 9) and ( 10) is decreasing.Also, this behavior is consistent with the saturation of the O + velocity distribution function as shown in Fig. 1.The behavior of the actual O + total heat flux is related to the behaviors of both parallel and perpendicular heat fluxes; it increases very rapidly at very low altitudes, and above 2 R E the heat flux slowly increases with altitude.

O + and H + "thermalized" ion heat fluxes
Based on Schunk (1977, and references therein), Schunk and Watkins (1982) developed a 13-moment system of transport equations and described collisional and collisionless heat flows.This set of transport equations has been used in modeling the polar wind, and they obtained altitude profiles for higher order moments such as ion heat flux for subsonic and supersonic flows.Demars and Schunk (e.g., 1989, 1992, 1994) developed this set of transport equations and obtained a new set of generalized transport equations based on bi-Maxwellian ion velocity distribution.In their 1992 paper, they replaced the actual velocity in the definition of the higher order moments by the thermal velocity, and defined the "thermalized" ion heat flux as follows: In these equations, n and m are the ion number density and mass, and k is the Boltzmann constant.Demars and Schunk (1992) compared the values of the "thermalized" ion heat flux obtained by using the above Eqs.( 12)-( 17 ion heat flux obtained from the solutions to the 16-moment equations, in order to test the convergence of their transport solutions and to investigate the applicability of the 16moment expansion, and, hence, the 16-moment equations, to the polar wind flows.They indicated that if the values of the actual heat flow are smaller than the corresponding values of the "thermalized" heat flow, then the series expansion converges for that specific ion outflow, and, then, the 16-moment set of equations is a useful tool for studying that ion outflow. In this section, we calculate altitude profiles of the "thermalized" ion heat fluxes (presented in Figs. 3 and 4, dashed lines) for both H + and O + ions in polar wind and auroral regions by using Eqs.( 12)-( 17) in order to compare it with the corresponding profiles of the actual ion heat fluxes obtained by using the Barghouthi model (presented in Sect. 4,Figs. 3 and 4,solid lines).We conduct this comparison, similarly to Demars and Schunk (1992), in order to provide a recommendation about the suitability of using 16-moment equations to investigate ion outflows in the polar wind and auroral regions.Also, according to Eqs. ( 12)-( 17), it is clear that the "thermalized" ion heat flux is related to the properties of the core of the ion velocity distribution.In other words, in this section, we want to emphasize the report of Biddle et al. (1985) about the significant difference between the flow velocity and temperature of the velocity distribution core and the energetic tail.
Figure 3 presents a comparison between the actual heat fluxes (solid lines) and the "thermalized" heat fluxes (dashed lines) for both H + and O + ions in the polar wind region.Similarly Fig. 4 presents a comparison between the ion actual heat fluxes (solid lines) and the "thermalized" heat fluxes (dashed lines) for both H + and O + ions in the auroral region.The "thermalized" heat flux is an estimate of the relative significance of the bulk as compared to the tail of the ion distribution.When the "thermalized" heat flux is large compared to the actual heat flux, the actual heat flux can be considered to result from a small perturbation of the thermal population (Demars and Schunk, 1992), and for that case the 16-moment set of equations is likely a useful tool for studying "thermal" space plasmas that develop non-Maxwellian features.On the contrary, when the actual heat flux is significant or larger than the "thermalized" heat flux, the 16-moment approximation may not be adequate to investigate non-Maxwellian space plasmas.Based on this and the results in Figs. 3, and 4, the 16-moment equations may be appropriate to study the H + ion outflow in the polar wind in the altitude range 1.7-10 R E , O + ion outflow in the polar wind in the altitude range 1.7-4 R E , and H + ion outflow in the auroral region in the altitude range 1.2 R E -3 R E .It may not be adequate to use the 16-moment equations to discuss O + ion outflow in the polar wind at altitudes above 4 R E , H + ion outflow in the auroral region at altitudes above 3 R E , and O + ion outflow in the auroral region at altitudes above 1.2 R E .

Applicability of our results
We have used the Barghouthi model which is based on Monte Carlo simulation, to study ion outflow in a polar cap and auroral region environment.The model includes some significant simplifications, as do all models.What has primarily been left out in our work is the dynamics the electrons may have in certain situations.We do include the effect of electron pressure gradients, and therefore the ambipolar electric field, although for a rather simple description of the background electrons as a fluid at constant temperature.In the model, we have chosen to represent the electrons as having a temperature of 1000 K at all locations along the modeled field line.There is experimental support for such a behavior in the perpendicular electron temperature (Boehm et al., 1995), but in the auroral region electrons are significantly heated in the parallel direction (Carlson et al., 1998;Chaston, 2003).Conservation laws predict that, at least in the auroral zone, parallel electron temperature should decrease as perpendicular ion temperature increases (Jasperse et al., 2010b).In the auroral region, electron pressure gradients significantly modify the parallel electric field (Jasperse et al., 2010b).Our model is not able to handle such effects as double layers, quasistatic field-aligned electric fields and field-aligned electric fields associated with Alfvén waves, which may affect the ion outflow.There are codes, such as Jasperse et al. (2010b) as well as Vlasov codes (e.g., Watt et al., 2004) which can be used to study the evolution of such electric fields selfconsistently, but it is clearly outside the scope of our study.Handling these self-consistently is not possible in a particle simulation over the time and spatial scales we are studying.A Vlasov simulation of ion outflow in the downward current region in the vicinity of auroral arcs has been published by Hwang et al. (2009).They found that the interaction between electric field double layers and the ions would affect the ion outflow, though only in the very limited altitude region close to the average auroral acceleration region.We want to study the evolution of the ions over a broader altitude interval.
Furthermore, there is good reason to believe that there are large regions, in particular the magnetospheric cusp, where there is strong wave activity but no field-aligned electric fields.This wave activity is typically sporadic, but it has been shown that despite this the average properties of the waves well explain the average properties of the ions as they evolve along the field line (Waara et al., 2011).The large variability of observed wave electric fields and associated heating still means that we have a large uncertainty, and at this stage it is not meaningful to put too much emphasis on details.It is therefore useful to study a smoothed out, steady-state average field line.The net effect of variable heating and indeed sporadic field-aligned electric fields in limited altitude regions can be expected to have a relatively good similarity to this average result.In particular our simulation provides a good similarity to many observations (Barghouthi, 2008;Barghouthi et al., 2011).Another strong point with using our steady-state solution is that if we did not run our model until steady state is achieved, we would have to initialize the density profile with a realistic density profile, and this choice of density profile would strongly affect our results.
We ignored the effect of centrifugal acceleration.This effect is distributed along the entire flight path, as could be in principle long-range potential differences also (Jasperse et al., 2006a(Jasperse et al., , b, 2010a, b), b).For centrifugal acceleration, at least we know that in the high-altitude cusp and mantle it is significant, but the net effect is much smaller than the wave-particle interaction (Nilsson et al., 2008).For the magnetotail lobes (i.e., polar wind), it has been shown that the effect of centrifugal acceleration is significant but still very small (Nilsson et al., 2010).
By using a steady-state solution and assuming that the magnetic field falls as a dipole with distance, we also ignore the effect of convection (E×B drift).One effect of the E×B drift would be to mix the auroral plasma with the polar wind plasma; we intend to provide a clear, separate picture for ion outflow in both regions.Therefore, we have considered the case of small drift which occurs at quiet geomagnetic activity (low Kp).This agrees well with our steady-state solution.Another effect of the E ×B drift would be to change how the magnetic field falls off with altitude.If convection dominates over parallel velocity, this could provide an important modification.Once again we have chosen to keep the situation simple here.Future studies could be more elaborate.
In the case of high geomagnetic activity, we believe a three-dimensional model similar to Barakat and Schunk (2006) would be needed for an accurate description.This is beyond the scope of this study.A 3-D model which can reproduce reliable ion heat flow estimates would demand the use of a very large number of particles.Possibly a single field-line model subject to temporal variations associated with convection could be used instead.Our work provides a reference point for such future more advanced models.
We neglected the effect of Coulomb collisions of the ions.The wave-particle interactions are strong in the auroral region, making the effect of Coulomb collisions relatively unimportant.However, the same may not be true in the polar wind region.In the model, the effect of ion heating for H + in the polar wind is 4 orders of magnitudes smaller than the corresponding values for the auroral region.Moreover, the effect of Coulomb collisions is stronger at lower altitudes while the modeled ion heating based on Eqs. ( 6) and ( 7) is the weakest there.The effect of Coulomb collisions may change the shape of the ion velocity distribution (Pierrard and Barghouthi, 2006), which, in turn, affects the heat flux.
Therefore we present representative average outflow of ions in wave environments typical of the polar wind and auroral zone (where the latter is even more applicable to the cusp region than the main auroral oval).These can be considered to be somewhat smoothed out as compared to spatially/temporally transient regions.It still describes the typical evolution of ion distributions over large distances along the field line, with good similarity to many reported observations.

Discussions and conclusions
The obtained behavior of the ion heat flux has been discussed in terms of the ion potential energy due to gravity and the polarization electric field, perpendicular adiabatic cooling, and wave-particle interactions.Our main findings are the following: 1.In the polar wind region and for H + ions, the total heat flux is always positive.This means that the flow of energy is upward along the geomagnetic field lines, and energy is going from low altitudes to higher altitudes.from about 2 up to 8 R E , which means that the flow rate of energy into the upward direction is decreasing.However, above 8 R E it is increasing; that is, more energy has been converted into the upward direction.This is consistent with the H + ion velocity distributions in the polar wind.At low altitudes the width of the velocity distribution function is decreasing in the perpendicular direction and increasing in the parallel directions.At about 8 R E and above, this situation is inverted.This has been explained in terms of perpendicular adiabatic cooling and wave-particle interactions, where at low altitudes the effect of perpendicular adiabatic cooling was dominant.There is a competition between these two effects at intermediate altitudes, and at higher altitudes the effect of wave-particle interactions is dominant.
2. In the polar wind region and for O + ions, the actual total heat flux is always positive and increasing.At altitudes below 5 R E , the heat flux is increasing very rapidly, and the width of the velocity distributions is well pronounced in the perpendicular direction.This is due to the significant role of the effect of wave-particle interactions which dominated the effect of perpendicular adiabatic cooling.At altitudes above 5 R E , the heat flux is increasing very slowly; this is because the heating process turned out to be saturated.
3. According to the magnitudes of the H + and O + ion total heat fluxes in the polar wind, the amount of energy transferred from low altitudes to higher altitudes by O + ions is more than that transferred by H + ions.
4. In the auroral region, H + ion actual total heat flux is negative at very low altitudes (i.e., around 1.8 R E ).This means that energy is going downward along the geomagnetic field lines: energy is moving from high altitudes to lower altitudes.The actual total heat flux was increasing very rapidly in the altitude range (2-4 R E ) due to the effect of wave-particle interactions, which means that the amount of energy that has been transferred along geomagnetic field lines is increasing.At altitudes above 4 R E the heat flux is positive and decreasing very slowly due to the competition between perpendicular adiabatic cooling and wave-particle interactions.
5. The behavior of O + ion total actual heat flux in the auroral region is similar to its behavior in the polar wind except that the magnitudes are higher in the auroral region and tend to increase at a slow rate at altitudes above 5 R E in the polar wind region and at altitudes above 2 R E in the auroral region.
6.In the auroral region the magnitudes of H + and O + ion total heat fluxes are on the same order.That is, the amount of energy transferred from low altitude to higher altitude by both ions is similar.This is due to the role of O + potential barrier which blocks many O + ions from moving in the upward direction to reach higher altitudes.Specifically, the higher O + potential energy barrier compensates for the higher O + density at the lower boundary.
7. Heat flux magnitudes are small near the bottom (lower boundary) of the simulation region because the ion velocity distribution at those altitudes is either Maxwellian or very close to Maxwellian.
8. When the parallel and perpendicular heat fluxes have different signs at the same altitude (e.g., H + perpendicular and parallel heat fluxes in the altitude range (6-8 R E ) in the polar wind), this indicates that the parallel and perpendicular parts of the ion energy are being transported in opposite directions.This behavior is the result of many competing process between perpendicular adiabatic cooling, parallel adiabatic cooling, wave-particle interactions, finite gyroradius effect, and potential energy.9.By comparing with "thermalized" heat flux, a measure of the relative weight of the thermal populations as compared to the net heat flux, we investigated regions where the heat flux could be considered to be a perturbation of a basically bi-Maxwellian core distribution.In the H + polar wind this was the case in the altitude range 1.7-10 R E , for O + ion outflow in the polar wind in the altitude range 1.7-4 R E , and H + ion outflow in the auroral region in the altitude range 1.2-3 R E .
10.We believe that this work is potentially highly significant because most fluid models close on the heat flux moment: they assume a particular form of the heat flux in order to produce a complete set of equations (see Schunk, 1977, for a review).The practical reason for doing so is that there are large statistical errors on highorder distribution function moments derived from particle detectors.Jasperse et al. (2006bJasperse et al. ( , 2010b) devoted a significant amount of discussion to the choice of closure relations for their model.The current results imply that certain choices of heat flux closure cannot capture the physics of ionospheric ion outflow.As noted earlier, when the ion heat flux components are small relative to the values of the "thermalized" ion heat flux, the 16-moment approximation may be adequate to investigate the ion outflow.It is not recommended to use 16-moment equations for cases when the ion heat flux is not small and larger than the "thermalized" heat flux (i.e., for conic and toroidal velocity distribution functions).
According to our knowledge, this is the first study that discusses H + and O + ion heat fluxes in the collisionless polar wind and auroral regions, that takes into consideration the effects of gravity, polarization electrostatic field, the divergence of the geomagnetic field lines, and wave-particle interactions, provides altitude profiles for ion heat fluxes in polar wind and auroral regions, and sheds light on the contributions of the energetic tail to the ion heat flux.In addition, it reports that on modeling the collisionless outflow in terms of the thermal properties of the ion velocity distribution leads to a significant difference between the modeled values and the actual ones.It is recommended to keep considering the shape of the distribution function and the energetic tail, especially when dealing with higher velocity moments like ion heat flux.

Figure 1 .
Figure 1.Ion velocity distribution functions at different geocentric altitudes in the polar wind (two left panels) and auroral (two right panels) regions, for H + and O + ions.The ion velocity distribution is represented by equal values contours in the normalized velocity ( c , c⊥ ) plane, where c = [v − u]/ 2kT /m 1/2 , where T is the ion

Figure 3 .
Figure3:PW Figure 3. Altitude profiles for actual (solid lines, the line shows the magnitude of the actual heat flux, with red color indicating positive values and blue color indicating negative values) and "thermalized" heat fluxes (dashed lines) of H + ions (left panels) and O + ions (right panels) in the polar wind region.The heat fluxes considered here are perpendicular heat fluxes (top panels), parallel heat fluxes (middle panels) and total heat flux (bottom panels).The unit of the heat flux is erg cm −2 s −1 .The arrow and the symbol next to it represent a label for the corresponding line.