A&A 372, 686-701 (2001)
DOI: 10.1051/0004-6361:20010552

On the temperature profile and heat flux in the solar corona:
Kinetic simulations

S. Landi - F. G. E. Pantellini

Département de Recherche Spatiale, Observatoire de Paris-Meudon, 92195 Meudon Cedex, France

Received 29 January 2001 / Accepted 26 March 2001

Abstract
In the solar corona the collisional mean free path $\lambda$for a thermal particle (electrons or protons) is of the order of 10^-2 to 10^-4 times the typical scale of variation H of macroscopic quantities like the density or the temperature. Despite the relative smallness of the ratio $\lambda/H$, an increasingly large number of authors have become convinced that the heat flux in such a plasma cannot be described satisfactorily by theories which suppose that the local particle velocity distribution functions are close to Maxwellian. We address this question through kinetic simulations of the low solar corona by assuming that non thermal velocity distribution functions are present at the base of the corona. In particular, we show that if one assumes that the electron velocity distribution functions at the base of the corona have sufficiently strong suprathermal power law tails, the heat flux may flow upwards, i.e. in the direction of increasing temperature. Using kappa velocity distribution functions as prototypes for non thermal velocity distributions, we find that the heat conduction can be properly described by the classical Spitzer & Härm (1953) law provided the kappa index is $\rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$>$ }\;5$. This value is much smaller than the value previously found by Dorelli & Scudder (1999). In addition we show that, unless extremely strong power law tails are assumed near the base of the corona (i.e. $\kappa < 4$), a local heating mechanism (e.g. waves) is needed to sustain the temperature gradient between the base of the corona and the coronal temperature maximum.

Key words: Sun: corona - methods: numerical - plasmas - conduction

1 Introduction

In this paper we present results from a one dimensional kinetic model of a semicollisional electron-proton plasma plunged in a gravitational field. The model is especially suited for stationary (not necessarily static) flows and for Knudsen numbers $K=\lambda/H\rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$>$ }\;10^{-4}$, where H is a typical scale of variation of a macroscopic quantity, such as the density or the temperature and $\lambda$ the distance between two successive collisions of a typical particle in the system. A simplified version of the model has previously been used by Pantellini (2000) to simulate a one-species atmosphere in a constant gravitational field. As expected, the result was the formation of a stratified isothermal atmosphere with an exponentially decreasing density known as the barometric law. The fact that the barometric law could be recovered was the first confirmation of the fact that, despite being one dimensional, the model could correctly reproduce known results. More recently, we implemented a more sophisticated version of the model to simulate an electron-proton plasma confined to the space between two conducting plates held at different temperatures and not subject to any external force (Pantellini & Landi 2000). We could show that the thermoelectric field needed to ensure quasi-neutrality in our simulation compares quite well with results of Fokker-Planck calculations with all possible interspecies collisions included (Spitzer & Härm 1953). These encouraging results motivated us to use the model to address the question of the heat flux in the solar corona.

The motivation for applying the model to the solar corona stems from the fact that observations suggest that above the transition region the typical thermal Knudsen number $K_{\rm T}\equiv \lambda~\partial \ln T/\partial z$ is of the order 10^-3 or larger (e.g. Dupree 1972; Ko et al. 1997; David et al. 1998; Fludra et al. 1999). It has been demonstrated that such a value, despite being much smaller than unity, is large enough for the classical transport coefficients (obtained by applying the Chapmann-Enskog formalism to the Fokker-Planck equation) to become substantially modified because of the presence of high non thermal energy tails in the electron velocity distribution functions (e.g. Shoub 1983; Scudder 1992b). Whence the necessity of using a numerical model appropriate for the solar atmosphere above the chromosphere-corona transition region where high values of the thermal Knudsen number $\rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$>$ }\;$10^-3 are commonplace. The difficulty with the coronal plasma (and even for the solar wind plasma out to distances of the order of astronomical units) is that neither a collisionless model based on the Vlasov equation nor a fluid model with Spitzer-Härm transport coefficients provide a convenient framework for investigation. Unfortunately, acceptance of the postulate of the existence of non thermal electron velocity distributions in the solar corona introduces an infinite number of additional free parameters required to define the distributions at the boundaries of the system. However, it is common practice to generalize the standard Maxwell-Boltzmann velocity distribution function using kappa distributions (see Eq. (18) below) which have the substantial advantage of requiring one additional free parameter only (the $\kappa$ index). Depending on the value of the parameter $\kappa$ the distribution departs more or less significantly from a Maxwell-Boltzmann distribution due to the presence of a more or less large excess of high energy particles. At least two studies have already discussed the fate of electron kappa distributions in the solar corona under the action of collisions. Anderson (1994) shows that collisions do strongly affect density and temperature profiles obtained using Scudder's (1992b) collisionless approach. After assuming that collisions do merely introduce first order perturbations to the collisionless distribution function he finds that the actual perturbations are of order unity or larger showing that collisions need to be treated self-consistently. To a certain extent this has been done by Dorelli & Scudder (1999) who let collisions affect the first order term of the Legendre polynomial expansion of the electron distribution function without assuming that this term was small but with the assumption of all higher order terms being zero. However, as we shall demonstrate below, higher order Legendre terms cannot be neglected. For example we find that collisions substantially modify the collisionless temperature profile (this has been observed by Anderson 1994, as well) indicating that at least the second order Legendre term must be retained in the expansion. The approach of Lie-Svendsen et al. (1999) is not substantially different from that of Dorelli & Scudder (1999) since they also use a first order truncated Legendre expansions for the electron distribution function. According to Chapman & Cowling (1970), such a truncation is valid for $K_{\rm T} \ll 1$ (weak inhomogeneity assumption) but, as pointed out by Shoub (1983) and Anderson (1994), $K_{\rm T}~\rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$<$ }\;10^{-3}$ is probably a more appropriate condition for the first order truncated Legendre expansion to remain a justified approximation, especially in the case of non-thermal boundary conditions. Unfortunately, 10^-3 is a typical value for the thermal Knudsen number $K_{\rm T}$in the corona and the weak inhomogeneity assumption may be regarded as questionable. Assuming much stronger temperature gradients than the ones assumed in both the Dorelli & Scudder (1999) paper and in the present work, Lie-Svendsen et al. (1999) argue that the classical Spitzer & Härm (1953) heat flux adequately describes the heat flux in the lower solar corona if Maxwellian boundary conditions are chosen at both ends of the simulated plasma slab. All these Fokker-Planck based models eventually are affected by additional limitations. For example, Dorelli & Scudder (1999) use the standard hydrostatic equilibrium equation as a closure whereas, following Shoub (1983), Lie-Svendsen et al. (1999) use a contestable zero-gravity pressure equilibrium condition. We do not need such a fluid closure equation nor do we require the velocity distribution functions to be of any particular form. However, the principal advantage of our model stems from the fact that collisions are included self-consistently, even though their treatment is strongly simplified with respect to the complexity of collisions in a real plasma. We are unable to evaluate precisely the importance of the simplified treatment of the collisions on our results. However, the very fact that the transport properties measured in test simulations compare well with those predicted by Fokker-Planck calculations suggests that our simplified way of handling collisions allows us to retain most of the essential physics occurring in a non-magnetized plasma.

Since the simulation model we use has never been described in full, we shall devote the next section to doing so. Non-essential details of the algorithm are presented in Appendix A. A brief discussion of the differences and similarities between our model and conventional Fokker-Planck models is given in Appendix B. The derivation of some relevant quantities (density, temperature and heat flux) for a collisionless plasma is given in Appendix C.

2 The model

A qualitative sketch of the model is shown in Fig. 1.

Figure 1: Schematic illustration of the model for a proton-electron plasma plunged in a z aligned constant gravitational field. The particles' velocities are 3D even though the model is spatially 1D. E0 is a z aligned external electric field which is generally needed to ensure local quasi-neutrality. In some cases a constant electric field does not suffice for the system to be neutral everywhere. In this case an additional small inter-particle electric field is introduced (as shown on the left hand side of the figure) to compensate for these polarization effects. When two particles encounter each other, they may collide as described in Sect. 2.3. If a particle of species $\alpha$ hits one of the two boundaries at z=0,L it is injected back into the system according to prescribed velocity distributions $f_\alpha ^{0,{\rm L}}(v_{z},v_\perp )$ (cf. Sect. 2.4) where vz is the particle's velocity along the z axis and $v_\perp$ the absolute value of the particle's velocity in a plane perpendicular to z.

Open with DEXTER

The model is based on the numerical integration of the one dimensional motion (one space and three velocity components) of N protons and N electrons plunged in a constant gravitational field and an electric field which is generally needed to ensure quasi-neutrality everywhere in the system. Particles are confined to the interval $z\in [0,L]$ by two conducting plates. Each time a particle hits one of the plates, it is instantly reinjected into the system according to a user defined prescription (e.g. elastic reflection, contact with a heat bath, etc.). Whenever the world lines of any two particles meet they may (or not) make an elastic collision depending on the magnitude of their relative velocity. The functional form of the velocity-dependent collision probability strongly influences the macroscopic behavior of the plasma (cf. 2.3).

2.1 Equation of motion

Both the gravitational field acceleration $\vec{g}$ and the electric field $\vec{E}(z)=[E_0+\delta E(z)]\vec{\hat z}$ are directed along the z axis so that the equation of motion for a particle of species $\alpha$ can be written as

  

$${\rm d}v_{z}/{\rm d}t = -g + q_\alpha E(z)/m_\alpha$$ (1)

$$v_{z} = {\rm d}z/{\rm d}t, \;\;\; v_{x} = {\rm const.}, \;\;\; v_{y} = {\rm const.}$$ (2)

where $q_\alpha$ and $m_\alpha$ are the charge and the mass of the particles of species $\alpha$ (electrons or protons). As can be seen on the left hand side of Fig. 1, we assume that the electric field between adjacent particles is constant, eventually increasing (decreasing) discontinuously by an arbitrary amount $\epsilon$ at the position of each proton (electron) in the system as if the particles where a succession of condensator plates attracting or repulsing each other depending on the charges on each plate. We then take the electrostatic field felt by a given particle (or condensator plate) to be the mean of the electrostatic field on either side of the particle. Thus, if we number all the particles in the system from 1 (bottom particle) to 2N (top particle) it follows that particle i feels the electrostatic field $E=\bar{E_i}\equiv(E_i+E_{i+1})/2$ during the whole time interval $\delta t_i$ until the earliest of the three possible collisions between i and $i\pm 1$ and between i-1 and i-2. Note that only a collision between particles of different charges may modify the topology of the electrostatic field in the system, provided the interacting particles exchange their relative position during the collision. Since the electrostatic field is piecewise constant, we always integrate Eq. (1) using a constant E field.

  
2.2 The electric field

It is often possible to ensure a satisfactory charge neutrality in the system without the need for an interparticle electrostatic field $\epsilon$. In these cases, the same constant external electric field E₀ can be used for all particles in the system. Particles do not feel each other. However, in the general case, a constant electric field is not good enough, as the plasma may behave like a dielectric medium, where polarization effects are no longer negligible. In that case, the Poisson equation must be written in the form

$${{\partial \vec{D}}\over{\partial z}}\equiv {\partial\over{\partial z}} (\varepsilon_0 \vec{E}+\vec{P}) = \rho$$ (3)

where $\rho$ is the charge density, $\varepsilon_0$ the permittivity of the vacuum, $\vec{D}$ the electric displacement and $\vec{P}$ the polarization. If $\vec{P}$ is negligible or (and) independent of z Eq. (3) implies that the electric field $\vec{E}$ needed to ensure charge neutrality ($\rho=0$) is a constant $\vec{E}=E_0\vec{\hat z}$. This may be the case when the density or temperature gradients in the system are small (Pantellini & Landi 2000). In cases where polarization effects are not negligible, a better charge neutrality can be obtained by introducing a moderately strong interparticle electric field $\epsilon$ as shown in Fig. 1. Let E_i denote the electric field between particles i and i-1 and let E₀ be an external electric field, for example the gravitoelectric field (Rosseland 1924). Let us suppose that the boundaries at z=0 and z=L are conducting. It is well known that in this case a charged particle is attracted by the boundary as if there was a particle of opposite charge in the symmetric position behind the conductor. In this case, it is easy to see that the electric field E₁between the conducting wall at z=0 and particle number 1 as well as the field E_i+1 above particle i must be given by

  

$$E_1\,\,\,\,\, E_0 + (q_{2\rm {N}} - q_1)\epsilon/2$$ = (4)

$$E_i + q_i\epsilon, \;\;\; i\in\{1,2,...,2N\}.$$ E_i+1 = (5)

We emphasize that to simulate a plasma, $\epsilon$ should be taken to be small enough for the electrostatic energy between neighboring particles $\vert q_i q_{i-1}\,\epsilon \delta z\vert$ ($\delta z$ is a typical interparticle distance) to be much smaller than the typical kinetic energy of the particles. Otherwise particles of opposite charge become bounded and form atoms.

  
2.3 Collisions

A particle of the system shown in Fig. 1 can either collide with another particle or with one of the two conducting walls at z=0,L. The former is elastic, i.e. both total momentum and total energy of the colliding particles are conserved while the latter is not, as the walls reflect particles according to a prescribed velocity distribution function regardless of the particles' velocities before the collision. Let us discuss the case of a particle-particle collision first; we shall come back to the particle-wall collision in Sect. 2.4.

Two particles may collide if they simultaneously occupy the same position along the z-axis. During an elastic collision between two particles (labeled 1 and 2) the velocity changes according to the well-known rules (e.g. Landau & Lifshitz 1960)

  

$$\vec{v}_1^\prime {{m_1\vec{v}_1+m_2\vec{v}_2}\over{m_1+m_2}} + {{m_2}\over{m_1+m_2}} u\, \vec{n}(u)$$ = (6)

$$\vec{v}_2^\prime {{m_1\vec{v}_1+m_2\vec{v}_2}\over{m_1+m_2}} - {{m_1}\over{m_1+m_2}} u\, \vec{n}(u)$$ = (7)

where primed and unprimed velocities represent pre- and post-collision velocities, respectively and $u\equiv \vert\vec{v}_2-\vec{v}_1\vert$. We note that the orientation of the unity vector $\vec{n}$in velocity space is not specified by the requirement of the total energy (1 equation) and total momentum (3 equations) to be conserved during the collision, as the number of unknowns is 6 (the three velocity components for both particles). Let us use spherical coordinates to define the orientation of $\vec{n}$ in velocity space with $\theta^\prime$ being the angle between $\vec{n}$ and the z-axis and $\phi^\prime$ the angle between the projection of $\vec{n}$in the (x,y) plane and the x axis. The question is: how shall we choose the angles $\theta^\prime$ and $\phi^\prime$ for any given collision in the system? Given the rotational symmetry of our system around the z axis, it is quite natural to choose $\phi^\prime$ according to a uniform probability distribution in the interval $[0,2\pi[$. The answer is not as obvious concerning the angle $\theta^\prime$ since the system is not spherically symmetric. In the case of particles of equal mass (Pantellini 2000) the probability distribution for $\theta^\prime$ is specified by the requirement that the system (in its most simple configuration, e.g. without external forces and with elastic boundary conditions) must relax towards a stationary state where the particle velocity distribution function is isotropic. This condition implies the post-collision angle $\theta^\prime$ should be chosen according to (see Sect. IIC in Pantellini 2000)

$$\theta^\prime = \arccos(\sqrt{P}),~ {\rm with}~ P={\rm random~number}\in[-1,1].$$ (8)

It is straightforward to convince oneself that the same result holds in the case of particles of unequal mass.

Now, even though the orientation of the $\vec{n}$ in Eqs. (6) and (7) must be chosen according to the above probability distributions if one requires the relaxed particle velocities to become distribute isotropically, one is still free to decide whether or not two particles which encounter each other effectively make a collision. If one decides that there isn't a collision, the two particles just go through each other without changing their velocities. If one decides that there is a collision, we compute the new velocities of the particles using Eqs. (6)-(8) to determine $\vec{n}$. In general, one is allowed to decide if two encountering particles collide depending on the magnitude u of their relative velocity only. The collision probability cannot depend on the orientation of $\vec{n}$ as the relaxed state would no longer be characterized by an isotropic velocity distribution.

	Figure 2: Collision probability R for hard sphere type collisions (solid line) and for Coulomb type collisions (dashed line) as a function of the relative velocity u.
Open with DEXTER

Figure 2 shows two choices for the velocity dependence of the collision probability R on the relative velocity u. One may interpret R as the collisional cross-section. Accordingly, we call the case R=1, where particles collide at each encounter, the "hard spheres'' case and

 

$$R_{\alpha\beta}(u) = \left\{ \begin{array}{ll} 1 & {\rm if}~ u < u_{\alpha\beta}\\ (u_{\alpha\beta}/u)^4 & {\rm otherwise} \end{array}\right.$$ (9)

the "Coulomb'' case. In Eq. (9) the indices $\alpha$ and $\beta$ indicate that one is considering collisions between a particle of species $\alpha$ and a particle of species $\beta$ (e.g. electron-electron, proton-electron, proton-proton).

A detailed comparison of the present model with other numerical models based on the Fokker-Planck or the Boltzmann equation (e.g. Shoub 1992) is beyond the scope of the present paper. Qualitatively speaking, the justification of the model stems from the fact that the scattering cross-section due to the cumulated effect of distant encounters in a near-equilibrium plasma is proportional to $\propto u^{-4}$ (e.g. Chandrasekhar 1943). Accordingly, one may interpret one collision in our model as representing the cumulated effect of a large number of distant encounters in a real plasma. Given that the most widely used (and best justified) numerical model to simulate distant encounter dominated plasma are Fokker-Planck models, we discuss more thoroughly the relation between our model and Fokker-Planck models in Appendix B. In the paper by Pantellini & Landi (2000) we show that using cut-off velocities $u_{\alpha\beta}$ of the order of, or smaller than, the typical relative velocity between $\alpha$ particles and $\beta$ particles, our model gives results that compare well with results from Fokker-Planck calculations.

  
2.4 Boundary conditions

There are several ways of treating the problem of a particle hitting one of the boundaries at z=0 and z=L. One may, for example, let the particle rebound elastically by simply changing the sign of the z component of its velocity. In this case, total energy is exactly conserved. However, neither temperature gradients nor non-Max- wellian distribution functions can be simulated in such a system. Given that we are interested in situations where both temperature gradients and non-Maxwellian distributions are present, we shall use more sophisticated boundary conditions allowing the injection of an arbitrary velocity distribution function. The prescription is as follows. Each time a particle of species $\alpha$ hits one of the boundaries, it is reinjected into the system following a specified velocity distribution function $f_\alpha^{0,{\rm L}}(\vec{v})$. This implies that in a stationary state, and apart from statistical fluctuations, the bulk velocity along z must be zero everywhere. We further assume that particles are injected following isotropic velocity distribution functions, i.e. $f_\alpha^{0,{\rm L}}(\vec{v}) = f_\alpha^{0,{\rm L}}(v)$, where $v=\vert\vec{v}\vert$ is the magnitude of the velocity of the injected particle. Accordingly, the theoretical flux of particles coming from the boundary with velocity $\vec{v}$ in the magnitude interval $[v,v+{\rm d}v]$ and orientation with respect to the z axis in the range $[\theta,\theta+{\rm d}\theta]$ is given by

$${\rm d}F_\alpha^{0,{\rm L}}(v,\theta)\, =\, v\cos\theta\:f_\alpha^{0,{\rm L}}(v)\:2\pi\:v^2\: \sin\theta {\rm d}\theta\:{\rm d}v.$$

This expression can be integrated separately for both the velocity v and the angle $\theta$, leading to the probability distributions P_v and $P_\theta$ of observing a particle entering the system from the boundary at z=0 with velocity V<v and $\theta$ in in the range $[0,\theta]$

  

$$P_\theta(\theta) = A_\vartheta \int_0^\theta \cos\vartheta\,\sin\vartheta\,{\rm d}\vartheta = \sin^2\theta$$ (10)

$$= A_{v} \int_0^v f(V)\,V^3\,{\rm d}V$$ P_v(v) (11)

where we have suppressed the species index $\alpha$ for readability and where A_v and $A_\theta$ are normalization constants such that $P_\theta(\pi/2) = P_{v}(\infty) =1$. Thus, each time a particle hits the boundary at z=0 it is reinjected following the probability distributions (10) and (11). This requires that the expressions (10) and (11) be solved for $\theta$ and v . For example, the angle $\theta$ is obtained by computing $\theta=\arccos(\sqrt{P})$ where P is a random number in the range [0,1]. Similarly v=P^-1_v(P) where P is again a random number in the range [0,1] and where the function P^-1_v is obtained by inverting (in most cases numerically) Eq. (11). Of course this procedure applies for both boundaries and all species.

  
3 Results

In the present simulations we consider a thin layer of a fully ionized electron-proton plasma plunged in a uniform gravitational field $g = GM_{\odot}/R_{\odot}^2$ where $M_{\odot}$ and $R_{\odot}$ are the solar mass and the solar radius, respectively while G is the universal constant of gravitation.

Following the reference paper by Dorelli & Scudder (1999), we assume typical temperatures and densities at the z=0 boundary to be $T_0^{\rm M} \equiv 5\times10^5~{\rm K}$ and $n_{\rm e}(0)=n_0^{\rm M} \equiv 10^8 ~{\rm cm^{-3}}$ (we shall use these quantities for normalization in the remaining of the paper). The typical temperature gradient between the two boundaries at z=0 and z=L is of the order of $1.4\times10^6~{\rm K}/R_{\odot}$. For a system length $L=0.1~ R_{\odot}$ this leads to an upper boundary temperature $T_{\rm L} = 6.4\times10^5 ~{\rm K}$. These parameters correspond to a thermal Knundsen number $K_{\rm T}\equiv \lambda_{\rm ee} (\partial T/\partial z)/T$ (T is the temperature and $\lambda_{\rm ee}$ the mean free path for electron-electron collisions) of the order 10^-4 to 10^-3, which is typical for the low solar corona in coronal holes (e.g. Ko et al. 1997; David et al. 1998; Fludra et al. 1999). The Fokker-Planck electron-proton collision frequency for such a plasma is given by Eq. (B.11). The same collision frequency is obtained in our simulation model if the number of electrons (or protons) N is of the order $N\simeq 1000$ (cf. Appendix B), which is therefore a typical value for all the simulations presented in the paper.

In order to reduce computational time a proton-to-electron mass ratio $m_{\rm p}/m_{\rm e}=100$ has been chosen for all simulations. This can be done provided the relevant dimensionless parameter

$$\gamma \equiv {{g (m_{\rm p}+m_{\rm e})R_{\odot}}\over{2k_{\r... ... \approx{{g m_{\rm p} R_{\odot}}\over{2k_{\rm B}T_0^{\rm M}}}$$ (12)

is chosen to be the same as in the real world, i.e. $\gamma\approx 23$. We note that $\gamma$ is the length $R_{\odot}$ expressed in units of the isothermal scale height of the atmosphere. Thus if one assumes protons to be less massive than in the real world, one has to assume gravity g to be stronger than in the real world (i.e. a fictitious Sun more massive then the real Sun), so as to ensure $\gamma$ remains unchanged. On the other hand, as we shall see below, a thermoelectric field $E_{\rm T}$ is needed to ensure quasi-neutrality in a plasma with an imposed temperature gradient (cf. Eq. (17)). Since $E_{\rm T}$ does not depend on the mass of the particles, at least as long as $m_{\rm e}/m_{\rm p}\ll 1$, there is no reason to use the real mass ratio in the simulation, the only requirement being the condition $m_{\rm e}/m_{\rm p}\ll 1$.

In simulations it is generally convenient to suitably normalize all physical quantities. Thus, throughout the remainder of the paper, we shall assume that velocities are normalized to $v_0^{\rm M}\equiv\sqrt{2k_{\rm B}T_0^{\rm M}/m_{\rm e}}$, distances to the slab thickness L, time intervals to $t_0^{\rm M}\equiv L/v_0^{\rm M}$, electric fields to $E_0^{\rm M}\equiv m_{\rm e} (v_0^{\rm M})^2/(eL)$ and heat fluxes to $q_0^{\rm M}\equiv m_{\rm e} n_0^{\rm M} (v_0^{\rm M})^3$.

Distribution functions and moments are constructed by regularly sampling positions and velocities of the particles in the system. In practice, we sample positions in bins of width 0.03125 L and velocities in bins of width of the order of 0.4 times the thermal velocity of the given population. In a typical simulation, 10³ particles encounter some 10¹⁰times and the distribution functions are obtained by sampling positions and velocities every 10⁴ encounters. This sampling interval is roughly the time it takes for a thermal proton to cross the plasma slab, which is also an estimate of the time memory of the system.

The just described procedure allows the construction of density or heat flux profiles which are not yet normalized. In order to do so one has to determine the "real" number density (in ${\rm cm}^{-3}$) somewhere in the system, for example at z=0. This is impossible in a collisionless stationnary and quasi-neutral system where the absolute density is an arbitrary parameter which can be eliminated from the equations (e.g. from the Vlasov equation). In a collisional system, however, the number density is no longer an arbitrary parameter given its intimate (roughly linear) connection with, for example, the electron-proton collision frequency. Thus, by recording the electron-proton collision frequency somewhere in the system (in units of $1/t_0^{\rm M}$ and thus in ${\rm s}^{-1}$), say at z=0, one can determine the absolute density there, provided a relationship between density and collision frequency has been previously established in some way. Such a relationship may have been established experimentally by measuring the collision frequency in a real Maxwellian plasma as a function of temperature and density. As we shall discuss below, and in Appendix B, we much more pragmatically adopt the relationship provided by a Fokker-Planck model. In brief, our strategy goes as follows: we choose a number of simulation particles N such that the recorded collision frequency near z=0 corresponds to a typical Fokker-Planck collision frequency for a plasma with an electron (or proton) number density n(0) of about $n_0^{\rm M}=10^8 ~{\rm cm}^{-3}$. In practice we just ensure that the Knudsen number in our simulation and in the solar corona are the same, despite the fact that the number of particles N in our system is ridiculously small compared to the number of particles which populate the solar corona. Fortunately, only 10³ to 10⁴ particles are required to simulate the corona. A number $N \sim 10^5$ would already require a computational power well beyond present day computer capabilities.

In the following subsections we shall discuss the behavior of a slab of solar corona for three different kinds of boundary conditions. The thickness of the slab L is taken to be either $0.1~ R_{\odot}$ or $0.2~ R_{\odot}$. The temperatures (based on the second moment of the velocity distribution function) of the boundaries are adjusted to make the mean temperature and the temperature gradient of the system compatible with the previously prescribed plasma conditions. No energy sources or sinks are present in the system. Energy is injected at the boundaries in the form of kinetic energy of the particles. The only way of transporting energy in the system is through a collisional (or collisionless) heat flux which means that all other means (e.g. radiation, waves, internal energy of the particles) are excluded.

The first subsection is devoted to the simulation of the "classical'' case with thermalized (Maxwellian) boundary conditions. We shall see that even in this case the Spitzer and Härm heat flux (Spitzer & Härm 1953) is not able to sustain the prescribed temperature gradient over a distance larger than $0.1~ R_{\odot}$ or so. In the second subsection we shall discuss the case of non thermal velocity distribution functions at the lower boundary. These simulations show that the prescribed temperature gradient can be sustained without local heating provided the number of suprathermal particles is high enough. This number turns out to be much higher than suggested in previous works (e.g. Dorelli & Scudder 1999). In the last subsection we briefly discuss the case of both boundary conditions being non thermal. We consider this case as rather unphysical as it supposes a source of suprathermal particles somewhere above the base depending on the position of the upper boundary. We discuss this case mainly because in the collisionless studies (e.g. Scudder 1992b) and in the reference paper by Dorelli & Scudder (1999) the nonthermal distributions "survive", by construction, across the entire slab of plasma.

3.1 Maxwell-Maxwell boundaries

Figure 3: Maxwell-Maxwell boundaries: in the top panel the electron density (solid line) and the temperature (dashed line) profiles are plotted. The dotted profile in the top panel represents the Spitzer-Härm temperature profile assuming a spatially constant electron heat flux. The dark square on the right indicates the temperature of the upper thermostat, i.e. the temperature of the boundary at $z=0.1~ R_{\odot}$. The bottom panel shows the proton and electron heat flux. As expected, their relative strength is of the order $\sqrt {m_{\rm p}/m_{\rm e}}$. The Spitzer labeled profile in the bottom panel has been computed via Eq. (14) using the measured electron temperature and density profiles.

Open with DEXTER

For the simulations in this section we impose Maxwellian distribution functions at the boundaries, i.e.

$$f_\alpha^{\rm0,L}\left(v\right) \propto \left(\displaystyle{... ...e}^{-\scriptstyle{\frac{m_\alpha v^2}{2k_{\rm B}T_{\rm0,L}}}}.$$ (13)

The normalized temperature and density profiles as well as both the electron and proton heat flux profiles are shown in Fig. 3. As expected, the electron to proton heat flux ratio is roughly equal to $\sqrt {m_{\rm p}/m_{\rm e}}$, which means that the energy is predominantly transported by the electrons. The dotted temperature profile in the upper panel has been obtained by assuming that the electron heat flux $q_{\rm e}$ between the two boundaries is constant and that the heat flux is given by the Spitzer & Härm formula (Spitzer & Härm 1953), which under the above conditions can be approximated by $q_{\rm e} \propto T_{\rm e}^{5/2} \partial T_{\rm e}/\partial z$ (cf. Eq. (14) and Appendix B). From the figure, it appears that if the upper thermostat is located at a distance of $0.1~ R_{\odot}$ from the surface, the measured profiles are in good agreement with the Spitzer-Härm predicted profiles, which seems to indicate that a classical Spitzer-Härm heat flux can sustain the given temperature profile up to a height of $0.1~ R_{\odot}$ without the need for some sort of local heating mechanism (e.g. dissipation of MHD waves). However, before accepting this statement as definitive, we have to ensure that the simulated plasma has the characteristics of the coronal plasma. In order to do so, we have to compute the "real" density of the simulated plasma and compare it to the density of the lower solar corona.

In a collisional plasma, the collision frequencies depend on the density (the higher the density the higher is the rate at which a particle undergoes collisions). On the other hand, the gravitational timescale $\sqrt{L/g}$ does not depend on density, which means that density is not just a free parameter as in the collisionless case (cf. Appendix C). We may then estimate the "real" density of the simulated system from the measured electron-proton collision frequency. The measured electron-proton collision frequency in the simulation is approximately $\nu_{\rm ep} = 717 /t_0^{\rm M}$ near the bottom at z=0. This means that in the average an electron collides 717 times with a proton during the time interval $t_0^{\rm M}$. Based on the measured collision frequency we may now determine the (unknown) number density n of the simulated plasma. In order to do so we make a slight detour in the field of Fokker-Planck models by observing that in Fokker-Planck models of a close to equilibrium plasma, with temperature T, collision frequency $\overline{\nu}_{\rm ep}^{\rm FP}$ and density n are intimately connected through Eq. (B.11) (cf. Appendix B). Thus, by using Eq. (B.10) with $\nu_{\rm ep} = 717 /t_0^{\rm M}$ and setting $\overline{\nu}_{\rm ep}^{\rm FP} = \overline{\nu}_{\rm ep}$ in Eq. (B.11) we can estimate the density at the bottom of our simulation region to be $n_0=1.01 n_0^{\rm M} \approx 10^8~{\rm cm^{-3}}$ which is an acceptable value for the low solar corona.

Similarly, we may compare the electron heat flux observed in our simulation with the Fokker-Planck heat flux for a fully ionized electron-proton plasma (e.g. Spitzer & Härm 1953)

$$q_{\rm SH} = -3.19\times10^{-3}\,\frac{n_{\rm e}k_{\rm B}^2T}... ...line{\nu}_{\rm ep}^{\rm FP}} \nabla T {\rm\ \ \ (SI\ units).}$$ (14)

Assuming, as above, that $\overline{\nu}_{\rm ep}^{\rm FP} = \overline{\nu}_{\rm ep}$ with $\overline{\nu}_{\rm ep}$ given by Eq. (B.10) and taking the temperature gradient measured in Fig. 3, one finds $q_{\rm SH} = -1.0\times10^{-3}q_0^{\rm M}$, which is in good agreement with the heat flux $q_{\rm obs} = -0.9\times10^{-3}q_0^{\rm M}$ measured in the simulation.

We emphasize that even in the collisionless case there is a heat flux flowing from the hot (upper) to the cold (lower) boundary. The collisionless heat flux $q_{\rm NC}$ can be computed analytically by applying Liouville's theorem (e.g. Landau & Lifshitz 1960) to the electron and proton distributions in constant gravitational and electric fields and subject to the above-specified boundary conditions, i.e. $T_{\rm L} = 1.28 T_0$ with $T_0=T_0^{\rm M}$. Given that the net particle flux is zero, the neutralizing electric field is nothing but the familiar gravitoelectric field (Rosseland 1924)

$$E_{\rm g}= {{g\,m_{\rm p}}\over{2e}} \left(1-{{m_{\rm e}}\over{m_{\rm p}}}\right).$$ (15)

In general, for Maxwell-Maxwell boundaries the collisionless electron heat flux is given by (see Appendix C)

$$q_{\rm NC} = \sqrt{{8 \over\pi}}\, \frac{n_{\rm L}k_{\rm B}^{... ...e}^{1/2}} \left(T_0T_{\rm L}^{1/2} - T_{\rm L}T_0^{1/2}\right)$$ (16)

where the electron density at the upper boundary $n_{\rm L}$ is a complicated function of all other parameters, which turns out to be $n_{\rm L} = 0.085 n_0$ for the present case. In the particular case shown in Fig. 3 we have $n_{\rm L}\approx 0.1 n_0$ from which we obtain a theoretical collisionless heat flux $q_{\rm NC} = -8.9\times10^{-3}q_0^{\rm M}$ which is much stronger a flux than either $q_{\rm SH}$ or $q_{\rm obs}$. The collisionless heat flux is an upper limit for the electron transported energy flux.

As already stated, in the collisionless regime $E_{\rm g}$ is the charge neutralizing field. However, in the collisional regime the total electric field is generally made of the sum of $E_{\rm g}$ and the thermoelectric field (see e.g. Golant et al. 1980; Hinton 1983)

$$E_{\rm T}= -{\alpha\over e}~{\partial\over{\partial z}}(k_{\rm B}T)$$ (17)

where $\alpha$ is a constant of order unity. In our simulations, $\alpha$ lies in the range 0.7 to 0.9 when $u_{\alpha\beta}$ in Eq. (9) is smaller than the typical relative velocity between particles of species $\alpha$ and particles of species $\beta$ (Pantellini & Landi 2000). Fokker-Planck models including all kind of collisions (electron-electron, electron-proton and proton-proton) predict $\alpha = 0.71$ (Spitzer & Härm 1953). The thermoelectric field required to ensure quasi-neutrality in the simulation of Fig. 3 is $E_{\rm T} = -0.11$, which is small (but not negligible) compared to the gravitoelectric field $E_{\rm g}=1.17$.

	Figure 4: Maxwell-Maxwell boundaries: Same format as Fig. 3.
Open with DEXTER

The question one may now ask is whether the above results depend on the position of the upper boundary or not. Even if the exact position of the temperature maximum is unknown it is likewise located somewhere between 0.2 and $1 ~R_{\odot}$ above the solar surface (e.g. Ko et al. 1997; David et al. 1998). Without going up to such heights where the zero mass flux hypothesis may no longer be valid, we may just ask the question of what happens if the upper boundary is located at twice the distance, i.e. $L=0.2~ R_{\odot}$. For coherence, the temperature of the upper thermostat has been chosen to be such that the temperature difference between the two thermostats is twice the value in the first simulation. The results are shown in Fig. 4. Interestingly enough, they are not quite the same as for the smaller system of Fig. 3. One observes that the electron temperature gradient is too weak to smoothly connect to the temperature of the upper thermostat (dark square). The heat that flows from the upper to the lower boundary is simply too weak to sustain the prescribed temperature profile: a local heating mechanism (waves?) is required in this case. The lower panel also shows that the observed electron heat flux is rather badly approximated by the Spitzer-Härm heat flux formula Eq. (14), suggesting that even in the most favorable case of Maxwellian boundary conditions, Eq. (14) does not provide a sufficiently reliable estimate of the thermal heat flux in the lower solar corona.

  
3.2 Kappa-Maxwell boundaries

Figure 5: Kappa-Maxwell boundaries with $\kappa =4$. The top panel shows the density (solid line) and the temperature (dashed line). Note that the temperature of the left boundary is only $0.49 T_0^{\rm M}$ (instead of $T_0^{\rm M}$ used in the Maxwellian case of Fig. 3) to compensate for the strong collisional heating of the plasma with height. The dashed line reproduces the temperature profile of the thermal boundaries case shown in Fig. 3. The lower panel shows the heat flux profiles for electrons and protons as well as the total heat flux (protons + electrons). The collisionless electron heat flux has been computed using Eq. (20). Note that while the proton heat flows down the temperature gradient (the "classical" behavior), the opposite is true for the electrons.

Open with DEXTER

In this section we discuss a run with Maxwellian boundary conditions at z=L (as in the previous section) and kappa velocity distribution functions

$$f_{\kappa}(v) = A_\kappa\left[1+\frac{v^2}{\left(\kappa - 3/2\right)v_0^2}\right]^ {\displaystyle -\kappa - 1}$$ (18)

with

$$A_\kappa = \frac{n_0}{2\pi\left(\kappa - 3/2\right)^{3/2}v_0^... ...ight)}{\Gamma\left(3/2\right)\Gamma \left(\kappa - 1/2\right)}$$

as boundary condition at z=0. In order to ensure energy equipartition among species, we use the same $\kappa$ index for both protons and electrons and a velocity v₀ (reminiscent of the thermal velocity in a Maxwellian distribution) $\sqrt {m_{\rm p}/m_{\rm e}}$ smaller for protons than for electrons. The density and temperature profiles as well the v_z distribution profiles for electron and protons are shown in Fig. 5 for the case $\kappa =4$. The upper panel shows that the temperature profile increases very rapidly away from the lower boundary. The rapid rise in temperature is not the manifestation of the collisionless gravitational velocity filtration mechanism first described by Scudder (1992a) but is essentially due to collisional effects. This can be demonstrated easily by estimating the temperature gradient due to the velocity filtration mechanism. According to Scudder (1992a) the latter is given by

 

$$T^{\displaystyle\ast} (z) T_0\left[1+\frac{\psi(z)} {(\kappa - 3/2)k_{\rm B}T_0}\right]$$ =

$$T_0\left[ 1+{{gz/v_0^2}\over{\kappa-3/2}} {{m_{\rm p}+m_{\rm e}}\over{m_{\rm e}}} \right]$$ = (19)

where $\psi(z) = m_{\rm p} g z + e\phi_{\rm g}(z) = m_{\rm e} g z -e\phi_{\rm g}(z)$ is the total potential energy of a particle with $e\phi_{\rm g}(z) = (m_{\rm e}-m_{\rm p})g z/2$ being the gravitoelectric potential energy (cf. Eq. (15)). According to Eq. (19), the temperature rises linearly with height provided $3/2 < \kappa < \infty$. For $\kappa =4$ and $T_0=0.49T_0^{\rm M}$ one has $T^{\displaystyle\ast} (z) = T_0(1+1.95 z)$ and at z=0.3 the temperature should be $T^{\displaystyle\ast} (0.3)= 1.59 T_0 = 0.78 T_0^{\rm M}$ only. Such a temperature is well below the temperature $T(0.3)\approx 1.1 T_0^{\rm M}$ observed in the simulation (see Fig. 5) from where we conclude that velocity filtration is not the principal reason for the rapid rise of the temperature profile inward from the z=0 boundary. Indeed, the strong inward heating of the plasma is neither due to the gravitational field nor to the temperature gradient imposed by the boundary conditions. The heating is essentially due to the effect of collisions on the $\kappa$ velocity distribution injected from the z=0 boundary. Of course the heating effect decreases with increasing $\kappa$ and in the limit $\kappa \rightarrow \infty$ (the Maxwellian case) it disappears completely, as shown in Fig. 5. Let us give a qualitative physical interpretation of the strong temperature gradient near the kappa boundary. If there were no collisions in the system, a collisionless electron heat flux would flow from the lower to the upper boundary (see Eq. (20)). However, the system is collisional and collisions do always act in the sense of a reduction of the heat flux. This is clearly visible in Fig. 5, where the electron heat flux is seen to decreases inwards from the right hand boundary over a distance of the order of 0.3L, inducing a strong heating of the plasma over this very same distance. Since protons become thermalized much more rapidly than electrons, the former carry energy in a "classical'' way, i.e. protons transport heat down the temperature gradient against the electron heat flux. We note that the total (proton + electron) heat flux is constant despite the fact that this is not so for individual species. Thus energy is transferred from electrons to protons and vice versa.

This peculiar behavior of the plasma near a kappa boundary does not show on the temperature profiles in Fig. 2 of Dorelli & Scudder (1999). The difference is due to the restrictions imposed upon the general shape of the electron velocity distribution functions by Dorelli and Scudder. In their paper, the general form of the electron velocity distribution function is a first order truncated Legendre polynomial expansion. Thus, at any given height z, the electron velocity distribution is a superposition of an isotropic kappa distribution function and an odd function of v_z (the first order correction) which does not affect the even moments of the velocity distribution function so that the temperature is by construction determined by the zero order distribution only, i.e. by the isotropic kappa distribution. Our simulation shows that restricting the Legendre expansion to the first order term only does not allow for a correct description of the temperature profiles especially for boundary conditions with low $\kappa$ indices. However, as we shall see in the remainder of this section, and in the next section, the unusual behavior of the heat flux (unusual from a fluid point of view) described by Dorelli & Scudder (1999), remains qualitatively valid for small values of $\kappa$.

The average total heat flux observed in our simulation is $q_{\rm obs} \approx 1.6\times10^{-3}q_0^{\rm M}$ and is mainly carried by the electrons (cf. Fig. 5). The very fact that $q_{\rm obs} >0$ means that energy flows upwards, i.e. from the cold to the hot thermostat. This may appear to be a surprising result as it seems to contradict the second law of thermodynamics (cf. the Introduction in Scudder 1992b). However, the behavior of our system can be described by the Boltzmann equation with a particular scattering operator, defined by the rules outlined in Sect. 2.3, and must therefore obey Boltzmann's H-theorem (Boltzmann 1872) and all fundamental laws of thermodynamics.

As a guiding reference for future discussion, we compute the collisionless electron heat flux $q_{\rm NC}$ by applying Liouville's theorem to the proton and electron velocity distribution functions imposed by the boundary conditions at z=0,L and the constraint of zero bulk velocity. Under these conditions, the charge-neutralizing electric field is precisely the gravitoelectric field $E_{\rm g}$ given by Eq. (15). Straightforward application of Liouville's theorem then leads to

$$q_{\rm NC} = \sqrt{{8}\over{\pi}}\, \frac{n_{\rm L}k_{\rm B}^... ...a - 2}\right)T^ {\displaystyle\ast}_{\rm L} - T_{\rm L}\right]$$ (20)

where the density $n_{\rm L}=n(L)$ is a complicated function of the other parameters of the problem (cf. Appendix C), $T^{\displaystyle\ast}_{\rm L} \equiv T^{\displaystyle\ast}(L)$ and $B_\kappa$ is a $\kappa$ dependent constant

$$B_\kappa \equiv \frac{\Gamma\left(\kappa - 1/2\right)}{\left... ...2\right)^{1/2} \Gamma\left(\kappa - 1\right)}\cdot \nonumber$$ (21)

With the parameters of the present simulation (i.e. $\kappa =4$, $T_{\rm L} = 1.28 T_0^{\rm M}$) one finds $n_{\rm L}=0.1 n_0^{\rm M}$, and $q_{\rm NC} = 5.4\times10^{-3} q_0^{\rm M}$ which, as expected, is stronger than the observed collisional electron heat flux observed in the simulation (Fig. 5). However, the general behavior of the system is neither that of a strongly collisional plasma with a Spitzer & Härm heat flux (1953) nor that of a collisionless plasma since the electron heat flux intensity is both spatially variable and substantially smaller than the collisionless value $q_{\rm NC}$.

The v_z velocity distributions for both electrons and protons at z=0.5L are shown in Fig. 6. From the figure it appears that while the proton distribution is essentially Maxwellian, the electron distribution has still substantial suprathermal tails, particularly for large positive velocities $v_z \rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$>$ }\;2.5$. The obvious reason is that the collisional cross section for a suprathermal proton vs. thermal electron collision is only weakly velocity dependent given that the relative velocity is always approximately $v_{\rm e}$, independent of the proton's velocity. Thus suprathermal protons are efficiently thermalized by collisions with thermal electrons. This is not so for a suprathermal electron since its velocity with respect to either a thermal proton or a thermal electron is, by definition, larger than $v_{\rm e}$. Given that the collisional cross section decreases as the forth power of the relative velocity it then follows that the thermalization of the suprathermal electrons is much less efficient than the thermalization of the suprathermal protons.

Figure 6: Velocity distribution functions for kappa-Maxwell boundary conditions ($\kappa =4$) at height z=0.5L for the case $L=0.1~ R_{\odot}$ shown in Fig. 5. Velocities have been normalized to the local thermal velocity $v_T\equiv \sqrt{2k_{\rm B}T/m}$ of the relevant species. The dotted curves are Maxwellians normalized to the local temperature and density. The dashed curves are kappa distributions with $\kappa =4$ where n0 and v0 are determined by the local density and temperature (cf. Eq. (18)).

Open with DEXTER

Let us conclude this section with a short discussion of the simulation in the light of the Dorelli & Scudder (1999) model (we shall call it the DS model). As already stated, there are some qualitative similarities between the behavior of the plasma observed in our simulations and the behavior of the plasma in their model. However, the very particular form of the distribution function in the DS model implies that their temperature profiles do differ significantly from ours. As already stated, the difference stems from the fact that in the DS model the just described collisional heating near a kappa boundary is missing, essentially because by construction the temperature in the DS model is that of a $\kappa$ distribution function with the same $\kappa$ index throughout the whole plasma slab. The reason for the DS temperature profile not being the collisionless temperature profile is due to the fact that their electric field (which has not been computed explicitly by the authors) is not Rosseland's gravitoelectric field (cf. Eq. (15)), which happens to be charge neutralizing in the collisionless case only. Despite these substantial differences, we do observe an upward-directed heat flux for the $\kappa =4$ in accordance with the DS model which predicts an upward directed heat flux for $\kappa \rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$<$ }\;10$.

Figure 7: Kappa-Maxwell boundaries with $\kappa =4$. The format of the figure is the same as in Fig. 5. Note that the upper boundary is located at $z=0.2 ~R_{\odot}$ and the temperature of the upper thermostat is the same as for the Maxwell-Maxwell case of Fig. 4. Again, the temperature of the lower thermostat is $0.49 T_0^{\rm M}$, as in Fig. 5. The dashed line in the top panel reproduces the temperature profile for the thermal boundary case of Fig. 4. As in the shorter system of Fig. 5, heat flows upward in the temperature gradient.

Open with DEXTER

How sensitive are the above results on the rather arbitrary position of the upper boundary? Figure 7 shows that doubling the size of the system and changing the temperature of the upper thermostat accordingly does not change the conclusions in a very substantial way. The main difference is that the average temperature gradient is slightly reduced with respect to the shorter system. As a consequence, the temperature of the plasma near the upper boundary is clearly below the prescribed temperature of the boundary (dark square on the figure). Thus, even with a kappa index as small as $\kappa =4$ one has to invoke a local heating mechanism to sustain the prescribed temperature gradient up to the z=0.2 level. For comparison, if the system was entirely collisionless, the temperature near the upper boundary would be as high as $1.7 T_0^{\rm M}$, i.e. above the temperature of the boundary.

	Figure 8: Temperature and density profiles for kappa-Maxwell boundary conditions. Each profile corresponds to a different kappa index ranging from $\kappa =4$ (dotted line) to $\kappa=\infty$ (solid line). The heat flux associated with each profile are shown in Fig. 10.
Open with DEXTER

	Figure 9: Same as Fig. 8 for the case extending up to $0.2~ R_{\odot}$.
Open with DEXTER

The temperature profiles for different kappa indices of the z=0 boundary distribution functions are plotted in Figs. 8 and 9. For each run the temperature of the z=0 boundary has been adjusted to obtain equal mid-box temperatures and temperature gradients. The collisional heating near the z=0 boundary is clearly visible on all plotted profiles except the $\kappa=\infty$ case. Figure 8 shows that if the upper boundary (the source of energy) is located at $z=0.1~ R_{\odot}$ the system is able to sustain the prescribed temperature profile independently of the $\kappa$ index. On the other hand, Fig. 9 shows that if the upper boundary is located at a height $z=0.2 ~R_{\odot}$ the system is no longer able to sustain the $1.4\times10^6~{\rm K}/R_{\odot}$ temperature gradient unless some local heating is at work. Indeed, all temperature profiles reach the $z=0.2 ~R_{\odot}$ level with a temperature which is clearly below the value imposed by the boundary. We note in passing that the steepening of the temperature profiles above $z\approx 0.17~ R_{\odot}$ is not due solely to the vicinity of the hot boundary but also to the collisionless gravitational velocity filtration. The reason for the collisionless filtration to become more efficient above $z\approx 0.17~ R_{\odot}$ is that at such heights the density has become extremely low (of the order 10^-2 times the density at z=0) and collisions much less effective in thermalizing the electron distribution function which still have suprathermal tails at a non negligible level (cf. Fig. 6) going into the heating via the collisionless gravitational velocity filtration mechanism. Of course gravitational filtration does not work in the Maxwellian case, which is the reason for the temperature to grow more slowly for $z~\rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$>$ }\;1.7~ R_{\odot}$ in the Maxwellian case than in the $\kappa =4$ case (cf. Fig. 7).

Even though the temperature and density profiles appear to be quite similar over the major part of the simulation domain for all cases shown in Figs. 8 and 9, the transport properties are different. This is particularly evident for the heat flux.

	Figure 10: Heat flux measured in the simulations for kappa-Maxwell boundaries shown in Fig. 8. Solid line represents the Spitzer-Härm value for the Maxwell-Maxwell case (cf. Fig. 3).
Open with DEXTER

In Fig. 10 are plotted the values of the total heat flux observed in the simulations for different values of the kappa index. The horizontal solid line represents the total heat flux for the Maxwell-Maxwell boundaries case, i.e. the classical Spitzer & Härm (1953) heat flux for the given temperature gradient and plasma parameters. For the kappa-Maxwell simulations, the energy flows depend sensitively on $\kappa$; for $\kappa\rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$<$ }\;4$ the heat flux is positive and flows from the cold to the hot boundary. For $\kappa=5$ the heat flux is negative but its intensity is still significantly smaller than the classical Spitzer-Härm value. For $\kappa\rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$>$ }\;6$ the heat flux is essentially Spitzer-Härm. The index $\kappa_{\rm o}$ below which the energy flows in the upward direction can be determined in the collisionless limit using Eq. (20) and the condition $q_{\rm NC} \ge 0$ from where one obtains

$$\kappa \le \kappa_{\rm o} = {1\over 2} {{T_{\rm L}}\over{T_{\... ...rm B}T_{\rm L}}} +3 {{T_{\rm L}-T_0}\over{T_{\rm L}}} \right).$$ (22)

Plugging the parameters of the above simulations into Eq. (22) one finds $\kappa_{\rm o} \approx 12$ which is significantly larger than the value we find in Fig. 10. We conclude by noting that in the DS model, reversal occurs for $\kappa\approx 10$ which is also the value for which the collisionless heat flux reverses in the kappa-kappa case (see Eq. (24) below). Our simulations indicate that the effect of collisions on both the heat flux and the temperature profile is much stronger than suggested by the DS model.

  
3.3 Kappa-kappa boundaries

In this section we briefly discuss the case where the velocity distribution functions are kappa at both boundaries z=0 and z=L. The collisionless electron heat flux $q_{\rm NC}$ can be computed analytically using the general expressions given in Appendix C

$$q_{\rm NC} = \sqrt{{8}\over{\pi}}\, \frac{n_{\rm L}k_{\rm B}^... ...sqrt{T^{\displaystyle\ast}_{\rm L}} - \sqrt{T_{\rm L}}\right).$$ (23)

As for the kappa-Maxwell case, one can determine the limiting value $\kappa_{\rm r}$ below which heat flows in the direction of the temperature gradient by setting $q_{\rm NC} \ge 0$ in Eq. (23). The result is

$$\kappa \le \kappa_{\rm r} = {3\over 2} + {{\psi_{\rm L}}\over{k_{\rm B}(T_{\rm L}-T_0)}}\cdot$$ (24)

For the parameters we use (which are the same as in the Dorelli & Scudder 1999, paper) one finds $\kappa_{\rm r}\approx 10$ which, of course, is precisely the value predicted by the DS model. Indeed, the main effect of collisions in the DS model is to reduce the intensity of the heat flux, not its sign. We note that $\kappa_{\rm r}$ is always smaller than the heat flux reversal value $\kappa_{\rm o}$ found in the kappa-Maxwell case by an amount $T_{\rm L}/2(T_{\rm L}-T_0)$. This difference is due to the fact that heat flows more easily up the temperature gradient if there is no downward-directed suprathermal tail due to the kappa boundary condition at z=L.

We have opted not to present simulations with kappa-kappa boundary conditions for two reasons. The main reason is that from a conceptual point of view it does not make much sense to suppose that there is a generator of kappa distributions located at an arbitrary height L above the coronal base. How shall one choose this point? What is the most appropriate value of the kappa index there? Our approach consists of assuming that there is a mechanism (e.g. shocks) capable of generating suprathermal tails at the base of the corona, i.e. at a natural boundary of the solar atmosphere. This is not so for the fictitious upper boundary we suppose to be located at $0.1~ R_{\odot}$ (following Dorelli & Scudder 1999) or $0.2~ R_{\odot}$, given that the solar corona extends out to distances of the order of many tens of AU. Also, given the strong collisionality of the system, we found the choice of the Maxwellian distribution to be the most natural one (or the less artificial one).

Simulations with a kappa boundary located at much larger distances, where the wind is supersonic and essentially collisionless, may be realistic as non thermal distributions are systematically observed there. Such simulations may become possible in the near future.

We conclude this section by noting that in the DS model the lower and the upper boundary conditions are not independent of each other, as the zero order distribution function is supposed to be a kappa distribution, with the same index $\kappa$, in all points of the system. The constant kappa index assumption in the DS model finds its justification in the fact that the kappa index is known to remain unchanged in the collisionless case (Scudder 1992a). As a consequence, the above discussion on the choice of the upper boundary condition for our simulation is irrelevant in the DS model where the two boundaries cannot be treated separately due to the constant constant kappa index assumption.

4 Conclusion

There we summarize the most important aspects arising from our simulations. First of all, as already suggested by other authors, it appears that due to their rapid thermalization, on a scale much shorter than the density or temperature scale of height, proton velocity distributions are close to Maxwellian in the solar corona. Unless some nonthermal local ion acceleration mechanism (e.g. some kind of magnetic turbulence) is at work, the Spitzer & Härm (1953) theory provides a good description of the proton transport properties in the corona. This seems not to be the case for the electron velocity distributions which, if not Maxwellian at a given height (at z=0 in the simulation), remain non Maxwellian over distances greater than the scale of height of macroscopic quantities. More specifically, if the electron velocity distributions have suprathermal tails, there is no hope for the Spitzer-Härm theory to provide the correct value of the electron heat flux. On the other hand, heat flux density and temperature profiles cannot be described by Scudder's collisionless model either (Scudder 1992a). This has been demonstrated by Anderson (1994), who also showed that the effect of collisions on the velocity filtration model is not a minor effect, which means that linearized collisional operators are inadequate for the description of the transport properties in the corona. Dorelli & Scudder (1999) tried to overcome this difficulty by expanding kappa distributions in terms of Legendre polynomials. This approach has the advantage of not requiring the first order Legendre term f₁ to be small compared to the zero order term f₀. However, the limitation of the development to the first order term only turns out to be too restrictive because of the strong up-down anisotropy of the problem. In particular, the temperature profiles cannot be conveniently described using a Legendre expansion truncated after the first order term, given that the latter does not directly contribute to the temperature. This is not very surprising. Anderson (1994) already suggested that the collisionless temperature profile is strongly modified by collisions (see Fig. 4 in his paper). Here we show that the collisionless temperature profiles of kappa distributions are strongly modified by collisions and that the effect is strongest close to the boundaries where the kappa distributions are artificially maintained (the z=0 boundary in our simulations). The simulations show that the collisional heating of the plasma near a kappa boundary increases with decreasing kappa index. The scale height of the collisional heating is determined by the relaxation length of the electron velocity distribution function, which for coronal plasma conditions is much shorter than the assumed temperature scale height $[(\partial T/\partial z)/T]^{-1} \approx 0.4 ~R_{\odot}$. The heating of a plasma near a kappa boundary could not be observed by Dorelli & Scudder (1999) due to the limited impact on the temperature that collisions are allowed to have in their model. One of the key points of the DS model was to show that the predictions of the collisionless gravitational velocity filtration model (Scudder 1992a) for the corona are only weakly modified by collisions. For example, for a given coronal temperature gradient, heat flux reversal occurs at the same kappa value whether collisions are included or not, the only effect of collisions being a reduction of the heat flux intensity. In reality, the effect of collisions on the collisionless results happens to be much more destructive than suggested by the DS model. We find that in the coronal plasma, heat flux reversal already occurs at $\kappa\approx 4$. For $\kappa=5$, the heat flux is already of the Spitzer-Härm type whereas the DS model predicts strong departures from the classical heat conduction even for $\kappa~\rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$>$ }\;10$. Our simulations also suggest that velocity filtration alone is not capable of sustaining the assumed coronal temperature gradient of $1.4\times10^6~{\rm K}/R_{\odot}$ unless $\kappa < 4$. This means that unless some extremely intense source of suprathermal electrons exists near the base of the corona some local heating mechanism (e.g. waves) has to be at work between the base of the corona and the coronal temperature maximum, which we assume to be at a height $z \rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$>$ }\;0.2 ~R_{\odot}$. Thus, gravitational velocity filtration may be capable of sustaining the observed temperature gradient without substantial heating, provided electron distributions near the coronal base have strong suprathermal tails. Even if present observations of the corona do not allow us to exclude the presence of strongly non thermal electron distributions at low altitudes, it seems hard to imagine a mechanism (e.g. Fermi acceleration, Fermi 1954) capable of sustaining such distributions given the high collisionality of the plasma. The conclusion would eventually be different if kappa distributions were injected into the system from the top, i.e. from the solar wind where non thermal electron velocity distributions are commonplace. However, this is much more general problem which cannot be treated in zero mass flux and plane parallel approximation used in this paper.

  
Appendix A: The algorithm

The structure of the algorithm for advancing the particles of the system of Fig. 1 during a given time interval is very similar to the algorithm described in Pantellini (2000) for the case of hard sphere particles of equal mass and no charge. The main difference is that the presence of an electric field makes the integration of the equations of motion slightly more complicated (not all particles feel the same acceleration). In addition, a non negligible fraction of the simulation time must be spent in computing the charge-neutralizing electric field.

The algorithm can be summarized as follows:

1.

Initialize the velocities $\{\vec{v}_i(t=0)\}$ and the positions $\{z_i(t=0)\}$ of N protons and N electrons. Order the particles based on their height such that 0<z₁<z₂<...<z_2N<L. Make an initial guess for the external electric field E₀ and eventually choose a non zero value for the variable $\epsilon$ (cf. Eqs. (4) and (5)) if polarization effects are expected to be important;

2.

  Determine for each pair of neighboring particles, with indices i and i-1 and for $i=2,3,...\,2N$, the time interval $\delta t_i$ until their next collision. If particle i and i-1 do not collide in the future, we shall set $\delta t_i=\infty$. Also compute the time $\delta t_1$ and $\delta t_{2N+1}$until the next collision of the first and last particle with one of the boundaries. The equations of motion to be solved are Eqs. (1) and (2) with the electric field felt by the particle i being $E_i+q_i\,\epsilon/2$, where E_i is recursively given by Eqs. (4) and (5);

3.

Determine the time interval $\delta t_{\rm min}=\min\{\delta t_i\}$ until the next collision in the system. Let $I\in \{1,2,...,2N+1\}$ be the index of the particle making this collision;

4.

  Advance all particles through the time interval $\delta t_{\rm min}$;

5.

 Make the collision between particle I and particle I-1 if $I\in \{2,3,...,2N\}$ according to the prescriptions given in Sect. 2.3. If I=1 or I=2N+1 (particle-wall collision), draw a new velocity vector for the particle using, in particular, Eqs. (10) and (11);

6.

If the system's charge neutrality is not satisfactory, increase or decrease E₀ in order to reduce departures from neutrality. E₀ should not be corrected too often. Ideally one should not update E₀ before the system has become stationary. This can be a very long time since it is of the order of the longest macroscopic relaxation time. For a Maxwellian plasma, the characteristic timescales could be $\sqrt{L/g}$ or L/c, where c is the sound speed. This step is rarely performed;

7.

Repeat steps 2-5 until a given time level has been reached.

The computational time needed to go through one cycle is proportional to N. A more sophisticated version of the algorithm allows us to make the computational time be proportional to $\sqrt{N}$ instead.

  
Appendix B: Comparison with Fokker-Planck model

Let's consider a relaxed system where N particles of species $\alpha$ and N particles of species $\beta$ are uniformly distributed in a box of dimension L. Let's consider particles which have relative velocities in the spherically geometric velocity-space element $2\pi u^2{\rm d}u{\rm d}\mu$, where $\mu \equiv \cos\theta$ ($\theta$ being the angle between $\vec{u}$ and the z direction). The number of collisions per time unit experienced by a particle of species $\alpha$ with particles of species $\beta$ with relative velocities near $\vec{u}$ is then given by

$${\rm d}\nu_{\alpha\beta} = u\left\vert \mu\right\vert R_{\alp... ...eta}(u) f_{\alpha\beta}( u,\mu )\,2\pi u^2 {\rm d}u{\rm d}\mu$$ (B.1)

where $R_{\alpha\beta}$ is the collision probability function defined in Eq. (9) and $f_{\alpha\beta}$ is the distribution function for the relative velocities between $\alpha$ and $\beta$ particles. If the distribution functions for both species are Maxwellians with thermal velocities $v_\alpha = \sqrt{2k_{\rm B}T_\alpha/m_\alpha}$ and $v_\beta= \sqrt{2k_{\rm B}T_\beta/m_\beta}$, respectively, one has

$$f_{\alpha\beta}(u) = \frac{1}{\pi^{3\over2}v_{\alpha\beta}^3} \frac{N}{L}\,{\rm e}^{\displaystyle{-u^2/v_{\alpha\beta}^2}}$$ (B.2)

with $v_{\alpha\beta}^2 \equiv v_\alpha^2 + v_\beta^2$. Integration of Eq. (B.1) over all velocities and directions then gives the total collision frequency

$$\nu_{\alpha\beta} = \frac{4}{\sqrt{\pi}}\frac{N}{L}\frac{1}{v... ...u)u^3{\rm e}^ {\displaystyle{-u^2/v^2_{\alpha\beta}}}{\rm d}u$$ (B.3)

i.e.

$$\nu_{\alpha\beta} = \frac{1}{\sqrt{\pi}}\frac{N}{L}v_{\alpha\beta} \Psi (\tilde{u}^2_{\alpha\beta})$$ (B.4)

where $\tilde{u}^2_{\alpha\beta}\equiv u^2_{\alpha\beta}/v^2_{\alpha\beta}$ and $\Psi$ is the function

$$\Psi (x) \equiv [1-(1+x){\rm e}^{-x} + x^2\Gamma\,(0,x)]$$ (B.5)

with

$$\begin{eqnarray*}\Gamma (0,x) = \int_{x}^{\infty}{\rm e}^{-t}t^{-1}{\rm d}t \end{eqnarray*}$$

being the incomplete gamma function of order 0. Sample values for the function $\Psi(x)$ are given in Table B.1.

   

Table B.1: Sample values for the function defined in Eq. (B.5); $x=\infty$ correspond to hard spheres collision.

x	0.125	0.25	0.5	1	$\infty$
$\Psi(x)$	0.03256	0.0918	0.230	0.484	1

Let us now consider the case of a thermalized and fully ionized electron-proton plasma and let's choose $u_{\alpha\beta}$ in Eq. (9) by setting

$$u_{\alpha\beta} = v_{\alpha\beta}/\sqrt{2}$$ (B.6)

so that the dependency on $\Psi$ is the same for all kind of collisions (i.e. proton-proton, electron-proton, electron-electron). For $m_{\rm p} \gg m_{\rm e}$, we then have $v_{{\rm ee}} = \sqrt{2}v_{{\rm e}}$, $v_{{\rm ep}} \simeq v_{{\rm e}}$ and $v_{{\rm pp}} = \sqrt{m_{{\rm e}}/m_{{\rm p}}}v_{{\rm ee}}$ which leads to the familiar relationship between the collision frequencies

$$\nu_{{\rm ee}} = \sqrt{2}\nu_{{\rm ep}} = \sqrt{\frac{m_{\rm p}}{m_{\rm e}}}\nu_{{\rm pp}}.$$ (B.7)

It is common practice to define the collision frequency $\overline{\nu}_{\alpha\beta}$ as the rate of momentum exchange between particles of species $\alpha$ and particles of species $\beta$. For isotropic distribution functions and elastic collisions one has (Golant et al. 1980)

$$\overline{\nu}_{\alpha\beta} = \left\langle \frac{\partial}{\... ...k}} \left[u_{\rm k}\nu_{\alpha\beta}^*(u)\right]\right\rangle$$ (B.8)

where $\langle...\rangle$ means averaging over the relative velocities distribution Eq. (B.2) and $\nu_{\alpha\beta}^* \equiv (N/L)u\vert\mu\vert R_{\alpha\beta}(u)$ is the collision frequency for the particles in the velocity-space element $2\pi u^2{\rm d}u{\rm d}\mu$. Carring out the integration, using the collision frequency Eq. (B.4), leads to

 

$$\overline{\nu}_{\alpha\beta} \frac{2}{\sqrt{\pi}}\frac{N}{L}v_{\alpha\beta} [1-(1+\tilde{u}^2_{\alpha\beta}){\rm e}^{-\tilde{u}^2_{\alpha\beta}}]$$ = (B.9)

$$\nu_{\alpha\beta}{2\over {\Psi(\tilde{u}^2_{\alpha\beta})}} \left[1-(1+\tilde{u}^2_{\alpha\beta}){\rm e}^{-\tilde{u}^2_{\alpha\beta}}\right].$$ = (B.10)

On the other hand, the Fokker-Planck electron-proton transport collision frequency for a plasma with an electron density n and temperature T is known to be (e.g. Golant et al. 1980)

$$\overline{\nu}_{{\rm ep}}^{\rm FP} = \frac{ne^4}{3\varepsilon_0^2m^{1/2}_{\rm e} (2\pi k_{\rm B}T)^{3/2}}\ln \Lambda$$ (B.11)

where

$$\ln \Lambda=\ln \left[\frac{12\pi(\varepsilon_0k_{\rm B}T)^{3/2}} {n^{1/2}e^3}\right]$$

is the Coulomb logarithm. We emphasize that the collision rate given in Eq. (B.11) is based on the time interval it takes for a thermal electron's trajectory to become strongly deflected by the protons in the system. After setting $\overline{\nu}_{{\rm ep}}^{\rm FP} = \overline{\nu}_{\alpha\beta}$ we find the number of electrons N needed to simulate a plasma with a given Fokker-Planck collision frequency $\overline{\nu}_{\rm ep}^{\rm FP}$ to be

$$N = \frac{1}{24\pi}\frac{ne^4\ln \Lambda}{\varepsilon_0^2k^2_... ...2_{\rm ep}\right) {\rm e}^{-\tilde{u}^2_{\rm ep}}\right]}\cdot$$ (B.12)

For typical coronal conditions, e.g. $n = 10^8 ~{\rm cm^{-3}}$ and $T = 5\times10^5 ~{\rm K}$ one has $\overline{\nu}_{\rm ep}^{\rm FP} = 20.4 ~{\rm s^{-1}}$. Thus, in order to simulate a plasma slab of thickness $L=0.1~ R_{\odot}$ some $N \simeq 4000$ particles are required. We note that the number of particles needed to simulate a given plasma strongly depends on the choice of $u_{\alpha\beta}$ in Eq. (B.6): the smaller the ratio $u_{\alpha\beta}/v_{\alpha\beta}$ the larger N and the heavier the simulation. One should therefore choose $u_{\alpha\beta}/v_{\alpha\beta}$ as large as possible remembering that when $u_{\alpha\beta}/v_{\alpha\beta}\rlap{\lower 2.5pt \hbox{$\sim$ }}\raise 1.5pt\hbox{$>$ }\;1$ the macroscopic characteristics of the plasma (e.g. heat conduction, thermoelectric coefficient, etc.) are no longer those of a Coulomb collision dominated plasma because of the large number of hard-sphere type collisions involving all particles moving at relative velocities $u < u_{\alpha\beta}$. The properties of such a plasma may differ substantially from the properties of a real electron-proton plasma!

  
Appendix C: The collisionless limit

Let us consider a collisionless electron-proton plasma plunged in a gravitational field $g(z)= -\partial \phi_{\rm G}/\partial z$, which is not necessarily constant. At the boundaries, z=0 and z=L, particles are injected with kappa type velocity distributions with temperatures T₀ and $T_{\rm L}$ and kappa indices $\kappa_0$ and $\kappa_{\rm L}$, respectively. Due to the unequal mass of electrons and protons, an electric field $E(z)=-\partial \phi_{\rm E}/\partial z$ is needed to ensure local charge neutrality. The stationary distribution function f(v,z) for particles of mass m and charge q with a monotonical potential energy

$$\psi(z) = m\phi_{\rm G}(z) + q\phi_{\rm E}(z)$$ (C.1)

can be computed using Liouville's theorem, viz.

$$f\left(v,z\right) = \left\{ \begin{array}{lr} f_0\left(v,z\ri... ...aystyle\ast}\left(v,z\right) & ~~~v_z < -w \end{array}\right.$$ (C.2)

where $w^2 \equiv 2\left[\psi_{\rm L}-\psi\left(z\right)\right]/m$ and

  

$$f_0 = \left(\frac{m}{2\pi k_{\rm B}T_0}\right)^{3/2} \frac{A_0\Ga... ... \left[1+\frac{mv^2+2\psi}{(\kappa_0 - 3/2)2k_{\rm B}T_0} \right]^{-\kappa_0-1}$$ (C.3)

$$f_0^{\displaystyle\ast} =\left(\frac{m}{2\pi k_{\rm B}T_0^\ast}\r... ...^2+2\psi}{(\kappa_{\rm L} - 3/2)2 k_{\rm B}T_0^\ast}\right]^{-\kappa_{\rm L}-1}$$ (C.4)

and where we have set $\psi(0)=0$ and $\psi(L)=\psi_{\rm L}$. In Eq. (C.4) $T_0^\ast$ is the temperature of the downward-traveling kappa distribution at z=0 which is linked to the temperature $T_{\rm L}=T(z)$ via (Scudder 1992a)

$$T_0^\ast = T_{\rm L}\left[1 - \frac{\psi_{\rm L}} {\left(\kappa_{\rm L} - 3/2\right)k_{\rm B}T_{\rm L}}\right].$$ (C.5)

The problem has two unknown parameters, A₀ and $A_{\rm L}$, which are determined by the condition of zero bulk velocity and $n_{\rm L}=n(L)$. The electric field profile is obtained by imposing local charge neutrality. We shall show that this field corresponds to the field required to make the potential energy of protons and electrons to be equal (cf. Eq. (C.14)).

Let $\chi\left(v,\mu\right)$ be an arbitrary function of the absolute velocity v and $\mu = \cos\theta$, where $\theta$ is the angle between z direction and the velocity. The mean value of this function is then given by

$$\langle\chi\rangle = 2\pi\int_0^w\int_{-1}^1{\rm d}\mu v^2{\r... ...\infty\int_{-1}^{-w/v}{\rm d}\mu v^2{\rm d}v \chi f_0^{\ast}.$$ (C.6)

For example, the density is obtained by integrating of Eq. (C.6) with $\chi = 1$:

 

$$n\left(z\right) \frac{A_0}{2}\left[1+\frac{\psi} {k_{\rm B}T_0\left(\kappa_0 - 3/2\right)}\right]^{-\kappa_0+1/2}$$ =

$$\left[1 + {\cal K}\left(\kappa_0 - 1,\xi_0\right)\right]$$

$$\frac{A_{\rm L}}{2}\left[1+\frac{\psi} {k_{\rm B}T_0^\ast\left(\kappa_{\rm L} - 3/2\right)}\right]^ {-\kappa_{\rm L}+1/2}$$ +

$$\left[1 - {\cal K}\left(\kappa_{\rm L} - 1,\xi_{\rm L}\right)\right]$$ (C.7)

where

$${\cal K}\left(\kappa,\xi\right) \equiv \frac{2}{\sqrt{\pi}} ... ...ght)}\int_0^\xi\frac{{\rm d}x} {\left(1+x^2\right)^{\kappa+1}}$$ (C.8)

and

$$\xi_0 \equiv \sqrt{\frac{\psi_{\rm L}-\psi} {\left(\kappa_0 - 3/2\right)\left(k_{\rm B}T_0+\psi\right)}}$$ (C.9)

$$\xi_{\rm L} \equiv \sqrt{\frac{\psi_{\rm L}-\psi} {\left(\kappa_{\rm L} - 3/2\right)\left(k_{\rm B}T_0^\ast+\psi\right)}}\cdot$$ (C.10)

In the limit $k\rightarrow\infty$ we have $\xi\rightarrow\sqrt{(\psi_{\rm L}-\psi)/k_{\rm B}T}$ and ${\cal K}\rightarrow{\cal G}$, the standard error function

$${\cal G} \equiv \frac{2}{\sqrt{\pi}}\int_0^\xi {\rm e}^{-x^2}{\rm d}x.$$ (C.11)

On the other hand the zero bulk velocity condition $\langle v_z\rangle =0$ leads to

  

$$A_0 = A \left[\frac{\left(\kappa_0 - 1\right) \kappa_0\Gamma \lef... ...c{\psi_{\rm L}}{\left(\kappa_0 - 3/2\right) k_{\rm B}T_0}\right]^{\kappa_0 - 1}$$ (C.12)

$$A_{\rm L} = A \left[\frac{\left(\kappa_{\rm L} - 1 \right)\kappa_... ...0^\ast}\right]^ {\kappa_{\rm L} - 1} \left({\frac{T_0}{T_0^\ast}}\right)^{1/2}.$$ (C.13)

All the above calculations indicate that the density n(z) only depends on the charge and mass of the particles via the potential energy $\psi(z)$ . Thus, if we use the same $\kappa$ and same boundary temperatures T₀ and $T_{\rm L}$ for both electrons and protons the charge neutrality condition reduces to

$$m_{\rm p} \phi_{\rm G} + e \phi_{\rm E} = m_{\rm e} \phi_{\rm G} - e \phi_{\rm E}$$ (C.14)

which arises from where one easily computes the classical gravitoelectric field (Rosseland 1924)

$$\phi_{\rm E} = \frac{1}{2}\left(m_{\rm p} - m_{\rm e}\right) \phi_{\rm G}.$$ (C.15)

The constant parameter A in Eqs. (C.12) and (C.13) is determined by the condition $n(L)=n_{\rm L}$, i.e.

 

$$\frac{2n_{\rm L}}{\sqrt{T_0}}\frac{1}{(B_{\kappa_0}/\sqrt{T_{\rm L}^\ast}) + (B_{\kappa_{\rm L}}/\sqrt{T_{\rm L}})}$$ (C.16)

with

$$B_{\kappa_0} \equiv \frac{\Gamma\left(\kappa_0 - 1/2\right)} {\left(\kappa_0 - 3/2\right)^{\frac{1}{2}} \Gamma \left(\kappa_0 - 1\right)}$$ (C.17)

$$B_{\kappa_{\rm L}} \equiv \frac{\Gamma\left(\kappa_{\rm L} - 1/2\right)} {\left(\kappa_{\rm L} - 3/2\right)^{\frac{1}{2}} \Gamma\left(\kappa_{\rm L} - 1\right)}$$ (C.18)

and where the quantity

$$T_{\rm L}^\ast = T_0\left[1 + {{\psi_{\rm L}}\over{(\kappa_0-3/2)k_{\rm B}T_0}} \right]$$ (C.19)

has been introduced.

We can now compute explicitly the higher moments for the velocity distribution

function Eq. (C.2). For example, the parallel pressure (with respect to z) turns out to be

 

$$P \left(z\right) = \frac{A_0k_{\rm B}T_0}{2} \left[1+\frac{\psi\l... ...0}\right]^{-\kappa_0+3/2} \left[1+{\cal K} \left(\kappa_0-2,\xi_0\right)\right]$$

$$+ \frac{A_{\rm L}k_{\rm B}T_0^\ast}{2} \left[1+\frac{\psi\left(z\... ...{\rm L}+3/2} \left[1-{\cal K} \left(\kappa_{\rm L}-2,\xi_{\rm L}\right)\right].$$ (C.20)

If we impose $\kappa_{\rm L} = \kappa_0$ and $T_0^\ast = T_0$, using Eq. (C.7) and Eq. (C.20), we obtain the well known result (Scudder 1992a)

$$T^{\displaystyle\ast} \left(z\right) = T_0\left[1+\frac{\psi(z)}{\left(\kappa -3/2\right)k_{\rm B}T_0}\right].$$ (C.21)

Similarly we find the collisionless heat flux by setting $\chi = 0.5mv^3\mu$ in Eq. (C.6), whence

$$q = \sqrt{\frac{8}{\pi}}\frac{n_{\rm L}k_{\rm B}^{3/2}}{m_{\rm e}... ... -\left(\frac{\kappa_{\rm L} - 3/2}{\kappa_{\rm L} - 2}\right)T_{\rm L}\right].$$ (C.22)

The heat flux for the Maxwellian case given in Eq. (16) can be obtained in the limit $\kappa_0\rightarrow\infty$ and $\kappa_{\rm L}\rightarrow\infty$. Equation (20) is obtained in the limit $\kappa_{\rm L}\rightarrow\infty$ and Eq. (23) by imposing $\kappa_0 = \kappa_{\rm L}$.

References

 
• Anderson, S. W. 1994, ApJ, 437, 860 In the text NASA ADS  
• Boltzmann, L. 1872, Wien. Ber., 275 In the text  
• Chandrasekhar, S. 1943, Rev. Mod. Phys., 15, 1 In the text  
• Chapman, S., & Cowling, T. 1970, The mathematical theory of non-uniform gases, 3 edition (Cambridge University Press, Cambridge, England) In the text  
• David, C., Gabriel, A. H., Bely-Dubau, F., et al. 1998, A&A, 336, L90 In the text NASA ADS  
• Dorelli, J. C., & Scudder, J. D. 1999, Geophys. Res. Lett., 26, 3537 In the text  
• Dupree, A. K. 1972, ApJ, 178, 527 In the text NASA ADS  
• Fermi, E. 1954, ApJ, 119, 1 In the text  
• Fludra, A., del Zanna, G., Alexander, D., & Bromage, B. J. I. 1999, J. Geophys. Res., 104, 9709 In the text NASA ADS  
• Golant, V., Zhilinsky, A., & Sakharov, I. 1980, Fundamantals of plasma physics (John Wiley & Sons, New York) In the text  
• Hinton, F. 1983, in Basic Plasma Physics, ed. A. Galeev, & R. Sudan, vol. 1 (North-Holland, Amsterdam), 147-197 In the text  
• Ko, Y., Fisk, L. A., Geiss, J., et al. 1997, Sol. Phys., 171, 345 In the text NASA ADS  
• Landau, L., & Lifshitz, E. 1960, Course of Theoretical Physics, vol. 1, 1st edition (Pergamon Press, Oxford, England) In the text  
• Lie-Svendsen, Ø., Holzer, T. E., & Leer, E. 1999, ApJ, 525, 1056 In the text NASA ADS  
• Pantellini, F., & Landi, S. 2000, Ap&SS, in press In the text  
• Pantellini, F. G. E. 2000, Am. J. Phys., 67, 60 In the text  
• Rosseland, S. 1924, MNRAS, 84, 720 In the text NASA ADS  
• Scudder, J. D. 1992a, ApJ, 398, 299 In the text NASA ADS  
• Scudder, J. D. 1992b, ApJ, 398, 319 In the text NASA ADS  
• Shoub, E. C. 1983, ApJ, 266, 339 In the text NASA ADS  
• Shoub, E. C. 1992, ApJ, 389, 558 In the text NASA ADS  
• Spitzer, L., & Härm, R. 1953, Phys. Rev. A, 89, 372 In the text