A Review on Solar Wind Modeling: Kinetic and Fluid Aspects

Abstract

The paper reviews the main advantages and limitations of the kinetic exospheric and fluid models of the solar wind (SW). The general theoretical background is outlined: the Boltzmann and Fokker–Planck equations, the Liouville and Vlasov equations, the plasma transport equations derived from an “equation of change”. The paper provides a brief history of the solar wind modeling. It discusses the hydrostatic model imagined by Chapman, the first supersonic hydrodynamic models published by Parker and the first generation subsonic kinetic model proposed by Chamberlain. It is shown that a correct estimation of the electric field, as in the second generation kinetic exospheric models developed by Lemaire and Scherer, provides a supersonic expansion of the corona, reconciling the hydrodynamic and the kinetic approach. The modern developments are also reviewed emphasizing the characteristics of several generations of kinetic exospheric and multi-fluid models. The third generation kinetic exospheric models consider kappa velocity distribution function (VDF) instead of a Maxwellian at the exobase and in addition they treat a non-monotonic variation of the electric potential with the radial distance; the fourth generation exospheric models include Coulomb collisions based on the Fokker–Planck collision term. Multi-fluid models of the solar wind provide a coarse grained description of the system and reproduce with success the spatio-temporal variation of SW macroscopic properties (density, bulk velocity). The main categories of multi-fluid SW models are reviewed: the 5-moment, or Euler, models, originally proposed by Parker to describe the supersonic SW expansion; the 8-moment and 16-moment fluid models, the gyrotropic approach with improved collision terms as well as the gyrotropic models based on observed VDFs. The outstanding problem of collisions, including the long range Coulomb encounters, is also discussed, both in the kinetic and multi-fluid context. Although for decades the two approaches have been seen as opposed, in this paper we emphasize their complementarity. The review of the kinetic and fluid models of the solar wind contributes also to a better evaluation of the open questions still existent in SW modeling and suggests possible future developments.

1 Introduction

The existence of a corpuscular solar outflow with velocities of the order of 10² km/s was deduced by Biermann (1953) based on observations of cometary tails. The term “solar wind” was coined by Eugene Parker (1958a) to describe the supersonic expansion of the solar corona. The concept received a theoretical description the same year (Parker 1958b). The solar wind (SW) models describe the physical mechanisms that accelerate solar corona plasma to supersonic velocities, and account for the SW properties at one Astronomical Unit (AU). Indeed, in the Earth’s neighborhood the solar wind has been sampled quasi-continuously by numerous spacecraft from almost the beginning of the space era. One of the reasons for monitoring the solar wind is that its sudden variations have a large impact on the Earth’s middle and upper atmosphere, ionosphere and magnetosphere. In this paper we review the kinetic and fluid models of the solar wind, with emphasis on their historical evolution and the main theoretical concepts stemming from the classical transport theory. We are aware that an exhaustive review of solar wind modeling is an outstanding task, beyond the scope of our study. Therefore we adopt the historical approach, aiming to provide the reader with the evolution of ideas and theoretical concepts and approaches from which emerge kinetic and fluid models of the solar wind.

The main difference between the kinetic and fluid description results from the fundamental physical variables on which is based each representation. The velocity distribution function (VDF) of each species, f _s (r, v, t), provides rich information on the plasma dynamics; it is the key-variable of the kinetic theory both in plasmas and neutral gases. It gives at time t the number of particles of species s within a 6-dimensional volume of the phase space, d r d v, centered on position r and velocity v at the time t. The computation of the VDF of each species from a kinetic equation is sometimes very difficult. The fluid description, based on the lower order moments of f _s , i.e. the density n _s (r, t), the bulk velocity u _s (r, t), the temperature T _s (r, t), the pressure tensor, $\underline{\underline{{p_s}}}$, is simpler and often sufficient for interpretation of experimental results.

Nevertheless, a thorough theoretical description of the physics governing the dynamics of an ensemble of ionized molecules moving almost without collisions in the gravitational and electromagnetic field must necessarily take into account the smallest scales and the corresponding velocity distribution functions. Therefore, the kinetic formulation should not be overlooked but comprehensively developed. In this paper we discuss the merits and limitations of the fluid and the kinetic approach, applied to the problem of the supersonic expansion of the solar corona. The complementarity of the two has been also discussed recently by Parker (2010). The models discussed in the following chapters are adapted to plasma and electromagnetic/gravitational field parameters typical for our Sun; their application can be extended to the case of stellar winds (Meyer-Vernet 2007).

The validity of any physical description or model of a dynamical system is determined by the length and time scales of the physical processes controlling the state of the system at the smallest scales. In a dense neutral gas or dense non-magnetized plasma the latter are determined by the binary collisions between molecules (e ⁻, H ⁺, etc.). In plasma the smallest spatio-temporal scales depend on the density, temperature, as well as on the properties of the embedded electromagnetic field. In neutral gases and dense plasmas the basic length scale is the mean-free-path, λ _c —the distance traveled by molecules between two collisions; the basic time scale is the collision time, τ _c —the time interval between two encounters. In plasmas an additional fundamental scale is the Debye length, λ _D , the screening distance for the Coulomb force. When the plasma is magnetized the (Larmor) radius of gyration, R _L and the gyration period T _L are also fundamental scales. For a broad range of densities, thermal energies and ambient electromagnetic field, the spatial and temporal kinetic scales, typical for the dynamics of individual particles, can be very small compared to the conventional, fluid macroscopic scales, but sometimes their ranges overlap, giving rise to complex, multi-scale processes. In this case the validity of the fluid approach becomes questionable.

Several reviews on solar wind modeling have been published in the past. Two of the most recent are the works of Marsch (2006) and the book by Meyer-Vernet (2007), where the basic concepts of solar wind modeling are illustrated and discussed from various perspectives. Marsch (2006) emphasizes the role of wave-particle interactions and their kinetics for solar wind acceleration and dynamics; Meyer-Vernet (2007) gives a physical perspective on the kinetics of solar wind models. In the following sections we provide the reader with a historical perspective on the solar wind modeling for over half a century, emphasizing the classical transport theory, a domain where the authors of this review contributed the most. A briefer account of this history can be found in Lemaire and Pierrard (2001).

The paper is organized as follows: in Sect. 2 we outline the theoretical background by discussing the fundamental equations to be considered in solar wind modeling. It shows that the kinetic and fluid approaches have the same theoretical root: the Boltzmann equation for gases and plasmas. Section 3 gives an overview of the early solar wind models, and put them into an historical perspective. In Sect. 4 we discuss the exospheric solar wind models, emphasizing their recent evolution during the last decades. The merits and limitations of modern kinetic models are pointed out. In Sect. 5 we provide a technical review of the multi-fluid models of the solar wind. We discuss their fundamental properties and results. The list of mathematical definitions and symbols used throughout the paper can be found in Appendix 1. Appendix 2 includes a discussion on the Chapman–Enskog and Grad solutions of the Boltzmann equation. It outlines the limitations of these two approaches when the velocity distributions departs significantly from displaced Maxwellians as is the case for the solar wind plasma at larger heliocentric distances.

2 Fundamental Kinetic and Hydrodynamic Equations

The equations relevant for the solar wind modeling are the fundamental equations of plasma physics. In Sects. 4 and 5, we discuss particular solutions of these equations, applied to the supersonic coronal expansion.

2.1 The Boltzmann Equation

A formal derivation of the Boltzmann equation from Liouville’s theorem of statistical mechanics can be found, for instance, in the paper by Grad (1958) or in the monographs by Jancel and Kahan (1966, pp. 262–264) or Uhlenbeck and Ford (1963, pp. 118–138). A more intuitive derivation, close to the original heuristic argument of Boltzmann, can be found in the monographs by Chapman and Cowling (1970, pp. 46–68), Montgomery and Tidman (1964, pp. 4–8) and in several recent textbooks like those of Gombosi (1994) and Schunk and Nagy (2004).

Let f _s (v, r, t) d r d v denote the number of particles of species s inside a phase space volume d r d v located at (r, v), where r denotes position and v velocity. The time evolution of the phase space density f _s is then given by the continuity equation

$$ \frac{\partial{f_s}}{\partial{t}} = - \nabla_{\bf r} \cdot( {\bf v} f_s) - \nabla_{\bf v} \cdot ({\bf a}_s f_s) +\sum_{t=1}^{N}\left(\frac{\delta f_{st}}{\delta t}\right)_{\rm coll} $$

(1)

where the operators ∇ _r and ∇ _v denote derivatives with respect to spatial and velocity coordinates and a _s  = F _s /m _s is the acceleration, where F _s is the total external force acting on the particle and m _s is the particle mass.

The first term on the right-hand side expresses the change in the number of particles in the phase space volume because particles enter and leave (at the velocity v) the volume d r. The second term expresses the analogous change because particles enter and leave the volume d v in velocity space. The term denoted δf _st /δt _coll corresponds to the rate of change of the VDF f _s (r, v, t) due to collisions between species s and t. It includes self-collisions. In the most general case the VDF can also change due to other processes, e.g. ionisation and recombination. These processes may be important particularly in the solar corona but they are not discussed here. Note that (1) is nothing but the expression of the conservation of particles, and as such bears a close resemblance, both in form and content, to the standard continuity equation (20) that we shall encounter later on in the fluid description.

Assuming that the force F _s is either independent of v, or is the Lorentz force, the Boltzmann equation (1) may be cast in a concise form:

$$ \frac{{\fancyscript{D}}f_s}{{\fancyscript{D}}t}= \sum_{t=1}^{N}{J_{st}(f_s,f_t)} $$

(2)

where we use the notation ${\fancyscript{D}}/{\fancyscript{D}}t=\partial/\partial t + {\bf v} \cdot \varvec{\nabla}_{{\bf r}} + {\bf F}_s/m_s \cdot \varvec{\nabla}_{{\bf v}}$ and J _st (f _s , f _t ) denotes the right-hand side of the Boltzmann equation, or the Boltzmann collision integral, that is discussed in some detail in the following paragraphs.

In practical situations it is useful to define the “peculiar” velocity, ${\bf c}_s\,{\equiv}\,{\bf v}-{\bf u}_s,$ in addition to the actual velocity v (Grad 1958; Schunk 1977); u _s is the mean velocity of species s. After this change of variable and considering that the force term, F _s  = F _sG  + F _sEM , is the sum of gravitational, F _sG  = m _s g and electromagnetic, F _sEM  = e _s E + e _s v × B forces, the Boltzmann equation (2) takes the form (Schunk 1977):

$$ \begin{aligned} &\frac{\partial {f_s}}{\partial {t}} + ({\bf c}_s +{\bf u}_s)\cdot \varvec{\nabla} f_s -\frac{D_s {\bf u}_s}{Dt}\cdot \varvec{\nabla}_{{\bf c}_s}f_s - {\bf c}_s\cdot \varvec{\nabla} {\bf u}_s \cdot \varvec{\nabla}_{{\bf c}_s}f_s \\ & \quad +\left[ {\bf g} + \frac{e_s}{m_s}\left( {\bf E} + {\bf u}_s \times {\bf B}\right) \right]\cdot \varvec{\nabla}_{{\bf c}_s}f_s +\frac{e_s}{m_s}\left( {\bf c}_s \times {\bf B} \right) \cdot \varvec{\nabla}_{{\bf c}_s}f_s = \frac{\delta f_s}{\delta t}_{\rm coll}, \end{aligned} $$

(3)

where the convective derivative is defined by:

$$ \frac{D_s}{Dt}\equiv \frac{\partial }{\partial {t}}+ {\bf u}_s \cdot \varvec{\nabla}_{{\bf r}}. $$

In this formulation the m-order moment of the velocity distribution function is defined by:

$$ n_s({\bf r}, t){\fancyscript{{\bf M}}}_s^{(m)}({\bf r}, t) = \int \int \int {{{\bf c}}_s^m f_s({\bf c}_s, {\bf r}, t)}{\rm d} {\bf c}_s $$

(4)

where n _s(r, t) is the zero-order moment (m = 0) giving the number density.

The first order moment (m = 1) corresponds to the average velocity, the second-order moment determines the pressure tensor, etc. For convenience, a list of definitions of the moments of order m ≤ 3 and of their corresponding plasma macroscopic quantities is given in Appendix 1.

A key problem for solving the Boltzmann equation is the treatment of the collision term, δf _s /δt _coll in (2) and (3). Classical review papers and textbooks (e.g. the paper by Grad 1958, or the monographs by Chapman and Cowling 1970, pp. 56–66, or Jancel and Kahan 1966, pp. 264–267), give a formal derivation of the collision term, based on the hypothesis of molecular chaos and elastic binary collisions. The following expression is found:

$$ J_{st}(f_s)= \int\limits_{0}^{2\pi} \int\limits_{0}^{\pi} \int\limits_{-\infty}^{+\infty} \int\limits_{-\infty}^{+\infty} \int\limits_{-\infty}^{+\infty} \left( f'_s f'_t - f_s f_t\right) \left| {\bf g}_{st}\right| \sigma_{st}(g_{st}, \theta) \sin \theta d\theta d\phi d{\bf v}_{t} $$

(5)

where g _st  = v _s  − v _t is the relative velocity between particle s and $t; \theta$ is the zenithal/deviation angle of g _st , and ϕ is the azimuthal angle; σ _st is the differential collision cross section; primed and non-primed quantities are considered after and before, respectively, collisions between particles s and t ($s \equiv t$ in case of self-collisions).

The functional relation between σ _st , |g _st | and $\theta$ is determined by the potential energy of the system of colliding particles and by the collision kinematics. When the colliding molecules/particles can be approximated by rigid spheres, the differential collision cross-section, $\sigma_{st}(|{\bf g}_{st}|, \theta)$ is a constant.

Note that in the lower corona and the inner solar wind additional processes take place, like ionisation, recombination, charge exchange and radiation; they are not included in the collision term (5) and are therefore neglected in models of the more distant solar wind. Wave–particle interactions cannot be neglected when the energy density of the resonant electromagnetic waves (MHD, ion cyclotron, whistlers, etc.) is comparable to the kinetic energy of the plasma particles. The effects of wave–particle interactions is addressed in the review by Marsch (2006) and in monographs of the solar corona (e.g., Aschwanden 2009).

2.2 The Fokker–Planck Equation

In a plasma, in addition to binary encounters, the charged particles interact/collide through the long range Coulomb potential. The differential cross section of Coulomb collisions between particles with charge e _s and e _t is given by:

$$ \sigma_{st}(g,\theta)=\frac{e_s^2 e_t^2(m_s+m_t)}{4g^4 m_sm_t \sin^4(\theta/2)} $$

(6)

Therefore the Boltzmann collision term described by (5) diverges in the case of Coulomb interactions for small impact parameters, ${\theta}\,{\to}\,0,$ and the hypothesis of an infinitely small time of the binary interaction is no more justified. One can approximate, however, that each long range Coulomb collision leads to a tiny deflection of the particle trajectory, described by the Boltzmann integral (5), and then simply add the effect of all collisions (Rosenbluth et al. 1957). In order to avoid the divergence for small deflection angles the integral (5) is cut-off at some angle $\theta_{min}$ that depends on the Debye length. In a proton-electron plasma this relationship is given by (Montgomery and Tidman 1964):

$$ \sin\frac{\theta_{min}}{2} \approx\frac{e^2}{mg^2L_D} $$

(7)

where the Debye length:

$$ L_D=\sqrt{\frac{kT_e}{8\pi n_e e^2}} $$

is defined as the distance over which the Coulomb field of a test-charge is screened-off by collective effects of the plasma charges.

Following this line of thought, the distribution function f _s in the Boltzmann integral can be expanded in terms of ${\Updelta}{\bf v}_s=\frac{m_t}{m_s+m_t}\Updelta {\bf g},$ up to the second order. The collision term on the right hand side of the Boltzmann equation can be rewritten as:

$$ J_{st} = \frac{\partial }{\partial {v_i}}\left[ -A_i^{s t} f_s + \frac{1}{2}\frac{\partial }{\partial {v_j}}\left( B_{ij}^{s t} f_s \right)\right] $$

(8)

where tensor notation is used with summation convention over identical indices; i, j correspond to Cartesian components, s, t indicate species. The term (8) is known as the Fokker–Planck collision term; $A_i^{st}$ is the coefficient of dynamical friction and corresponds to slowing-down effects; the term $B_{ij}^{st}$ is the coefficient of diffusion in velocity space. By assuming an inverse square law for the inter-particle force and a small Debye length, $A_i^{st}$ and $B_{ij}^{st}$, from (8) can be written (Rosenbluth et al. 1957; Clemmow and Dougherty 1969, pp. 423–431):

$$ A_{i}^{s t} = \frac{1}{4 \pi}\left( \frac{e_s e_t}{\varepsilon_0 m_s} \right)^2 \frac{m_s + m_t}{m_t} \ln \Uplambda \frac{\partial }{\partial {v_i}} \int{\frac{f_t ({\bf v}')}{\left|{\bf v} - {\bf v}' \right|} d^3 v'} $$

(9)

$$ B_{ij}^{s t} = \frac{1}{4 \pi}\left( \frac{e_s e_t}{\varepsilon_0 m_s} \right)^2 \ln \Uplambda \frac{\partial^2}{\partial v_i \partial v_j} \int{{\left|{\bf v} - {\bf v}' \right|}f_t({\bf v}') d^3 v'} $$

(10)

where primes denote post-collision quantities and

$$ \Uplambda = \frac{12 \pi \left( \varepsilon_0 kT \right)^{3/2}}{n^{1/2} e^3} $$

is the plasma parameter, with ε₀ the electric permittivity of vacuum, k the Boltzmann constant, T the plasma temperature, n the plasma density, and e is the elementary charge. The form of the collisional integral given in (8)–(10) is precisely the same as the one derived by Landau (1936) for the kinetic equation of plasma state, from statistical mechanics arguments; therefore it is sometimes also called the Landau collision term.

2.3 The Liouville and Vlasov Equations

When the mean-free path of plasma particles is large compared to the characteristic spatial dimension of the system itself, and the time between collisions is larger than the characteristic time scale, and due to the long-range Coulomb interactions between charged particles, the many-body interactions cannot be neglected anymore in the description of plasma dynamics. The effect of collective encounters is not appropriately evaluated by the Boltzmann integral (5) describing binary, elastic collisions. The Liouville theorem then gives the appropriate theoretical framework for the description of plasma dynamics.

The fundamental concept of the statistical approach is the phase space probability density, D _N (X ₁, X ₂, ..., X _N , t), where ${\bf X}_i\equiv({\bf r}_i, {\bf v}_i)$ and N is the total number of particles. D _N denotes the probability that particle 1 is located in [r ₁, r ₁ + d r] and its velocity pertains to the interval [v ₁, v ₁ + d v], while particle 2 is located in [r ₂, r ₂ + d r ₂] and its velocity pertains to [v ₂, v ₂ + d v], etc. Thus D _N is defined in the 6N dimensional space of the positions and velocities of the N plasma constituents (electrons, protons, etc.). The Liouville theorem states that the statistical ensemble is represented in the 6N dimensional space by a “cloud” whose volume does not change with time, as in incompressible flows:

$$ \frac{\partial {D_N}}{\partial {t}}+\sum_{j=1}^{N}\left\{ \frac{\partial {D_N}}{\partial {{\bf r}_j}}\cdot \dot{{\bf r}}_j + \frac{\partial {D_N}}{\partial {{\bf v}_j}}\cdot \dot{{\bf v}}_j \right\} = 0 $$

(11)

where $\dot{{\bf r}}_j$ and $\dot{{\bf v}}_j$ denote time derivatives. Note that (11) gives only a formal expression of the Liouville theorem, in its most general form, also known as a master equation. General solutions of (11) have not been found yet. A discussion of the physical content of D _N and particular forms and solutions of the master equation (11) for different interaction potentials between particles can be found in Balescu (1963).

The key-elements of the statistical theory of plasma physics are described in classical monographs (e.g., Montgomery and Tidman 1964, Chap. 4–6; Balescu 1963, pp. 26–55). In the BBGKY (Born-Bogoliubov-Green-Kirkwood-Yvon) approach the strategy to solve (11) is based on a hierarchy of equations derived by integration of (11) in subdomains of the phase space. The BBGKY hierarchy of equations couples the lower order, or reduced, probability distribution functions to the higher order ones (see Bogoliubov 1962; also Montgomery and Tidman 1964, pp. 41–50). The reduced one-particle distribution function, f ₁, is computed by integration in the 6(N − 1) dimensional sub-domain defined by the coordinates and velocities of the other N − 1 particles:

$$ f_1 ({\bf X}_1)= \int\int \ldots \int{D_N({\bf X}_1, {\bf X}_2, \ldots, {\bf X}_N)d{\bf X}_2, d{\bf X}_3, \ldots, d{\bf X}_N} $$

(12)

f ₁ defined by (12) must be identified with the velocity distribution function described by the Boltzmann equation (1), (2) or (3). The BBGKY chain of equations relates f ₁(X ₁) with f ₂(X ₁, X ₂), f ₃(X ₁, X ₂, X ₃) and all the higher order reduced probability density. The first equation of this chain can be written in a simplified form, when there is no magnetic field and the factor $G=1/(n_0L_D^3) {\to}0$ (see, e.g., Montgomery and Tidman 1964, pp. 48–49):

$$ \frac{\partial {f_1}}{\partial {t}} + {\bf v}_1 \cdot\frac{\partial {f_1}}{\partial {{\bf r}_1}} -\frac{n_s}{m_s}\left[ \int \frac{\partial {\phi_{12}}}{\partial{{\bf r}_1}} f_1({\bf r}_2, {\bf v}_2) d{\bf r}_2 d{\bf v}_2\right]\cdot \frac{\partial {f_1}} {\partial {{\bf v}_1}} = 0 $$

(13)

with ϕ₁₂ the two-particle interaction potential energy and n _s the mean particle density. When (13) is valid, all the higher order corrections of the BBGKY hierarchy are equal to zero. Equation (13) is also known as the Vlasov equation.

Since the one-particle distribution function f ₁ from (12) and (13) has the same physical content as the velocity distribution function from the Boltzmann equation (2), the Vlasov equation is also called the Boltzmann equation without collisions. Note also that Balescu (1963) developed an alternative way to solve the Liouville equation (11) based on Prigogine’s (1963) method of diagrams for non-equilibrium statistical mechanics. Balescu’s approach is based on a direct solution for D _N (X ₁, X ₂, ..., X _N , t) from (11). When G → 0 the Prigogine-Balescu solution leads to the same Vlasov equation (13).

In the following we adopt the kinetic notation introduced in Sect. 2.1 and write the Vlasov equation in the standard form, taking into account a gravitational and magnetic field:

$$ \frac{\partial f_s} {\partial t} + \left[ \frac{q_s}{m_s}\left( {\bf E} + {\bf v} \times {\bf B} \right) + m_s {\bf a}_g\right] \cdot \nabla_{\bf v} f_s + {\bf v} \cdot \nabla_{\bf r}f_s = 0 $$

(14)

where a _g is the gravitational acceleration. In (14) the electric (E) and magnetic (B) fields are derived from the electric charge (ρ) and electric current density (j); thus (14) must be coupled to Maxwell’s equations:

$$ \varvec{\nabla} \cdot {\bf B} = 0 $$

(15)

$$ \varvec{\nabla} \times {\bf E} + \frac{\partial {{\bf B}}} {\partial {t}}= 0 $$

(16)

$$ \varepsilon_0 \varvec{\nabla} \cdot {\bf E} = \rho_c^{ext} + \sum_s{e_s n_s} $$

(17)

$$ \frac{1}{\mu_0} \varvec{\nabla} \times {\bf B} - \varepsilon_0 \frac{\partial {{\bf E}}}{\partial {t}} = {{\bf j}}_{ext} + \sum_s{{\bf j}_s} $$

(18)

where e _s is the electric charge of species s. The external electric charge and electric current densities¹, $\rho_c^{ext}$ and ${{\bf j}}_{ext}$, have been separated from n _s and j _s , the internal contribution of the plasma itself, given by the zero and first order moments of f _s (see (4) and Appendix 1). A Vlasov-Maxwell equilibrium is defined by solutions to the coupled system of equations (14)–(18) with the definitions (4).

In summary, the chain of equations described in this section illustrates the microscopic description of plasmas. The kinetic approach considers the micro-physics of particle dynamics, for each plasma component species. Therefore, the kinetic theory has the following advantages:

• Kinetic models describe the velocity distribution for each component species separately and treat self-consistently the coupling between the particle dynamics, external forces and the electromagnetic field;

• Collisions and wave–particle interaction can be included; relevant solutions have thus been developed based on the Fokker–Planck equation;

• The entire set of moments equations is satisfied by the solutions of the kinetic solution;

• The moments of steady state velocity distributions are analytical functions determined by the electromagnetic field and the VDF parameters at the boundaries of the integration domain;

• The kinetic solutions are self-consistently coupled to Maxwell’s equations;

The advantages and the limitations of solar wind models based on the kinetic theory are discussed in Sect. 4.

2.4 Plasma Transport Equations

There are physical situations, including the case of the solar wind, where the details of the velocity distribution functions are not measured directly. Only the spatial and temporal variation of the lower order moments can then be evaluated, the density, n _s (r, t), the temperature, T _s (r, t), the bulk velocity, u _s (r, t), the pressure, P(r, t), etc.

These observables can be also determined from the solutions of the Boltzmann, Vlasov or Fokker–Planck equations. Their temporal and spatial evolutions are derived by integrating the kinetic equations over the velocity space. One then solves the so-called “equations of change”:

$$ \int\int\int\Upphi_s({\bf v})\frac{{\fancyscript{D}}f_s}{{\fancyscript{D}}t}d{\bf v}= \sum_{t=1}^{N}\int \int \int {\Upphi_s({\bf v})J_{st}(f_s,f_t)} d{\bf v} $$

(19)

where $\Upphi_s({\bf v})$ is a generic notation standing for various powers of the velocity. From (19) one retrieves the partial derivative equations describing the spatio-temporal evolution of the moments of the VDF, also known as the plasma transport equations. If one replaces in (19) $\Upphi_s=m_s, {\Upphi_s=m_s}c_{si},$ and $\Upphi_s=\frac{1}{2}m_s c_s^2,$ where c _si are the components of the peculiar velocity c _s  = v − u _s (with u _s the average velocity of species s defined in Appendix 1) and after integration over velocity one obtains the continuity, momentum and energy transport equations for species s (Schunk 1977):

$$ \frac{\partial {n_s}}{\partial {t}}+ \nabla \cdot (n_s {\bf u}_s) = 0 $$

(20)

$$ m_s n_s \frac{D_s {\bf u}_s}{Dt} + \varvec{\nabla} \cdot \underline{\underline{{p_s}}} - n_sm_s{\bf g} -n_s e_s( {\bf E} + {\bf u}_s \times {\bf B}) = \frac{\delta {\bf M}_s}{\delta t}_{coll} $$

(21)

$$ \frac{D_s}{Dt} \left( \frac{3}{2} n_skT_s \right) + \frac{5} {2} n_skT_s\left(\nabla \cdot {\bf u}_s\right) +\varvec{\nabla} \cdot {\bf q}_s + \underline{\underline{{p_s}}}:\varvec{\nabla} {\bf u}_s = \frac{\delta E_s}{\delta t}_{coll} $$

(22)

where the right hand sides denote the collision terms, and n _s is the number density, u _s is the average velocity, $\underline{\underline{{p_s}}}$ is the pressure tensor, k is the Boltzmann constant, T _s is the temperature (see Appendix 1 for definitions of these macroscopic variables). In (20)–(22) the production and loss of particles are disregarded, which is important in the transition region and below. If these terms are included, the right-hand side of (20) becomes nonzero, and additional terms are added to (21)–(22). In (20)–(22) the force term is of electromagnetic and gravitational origin. Since m _s , m _s v _i , and $\frac{1}{2}m_s v^2$ are conserved in binary, elastic collisions, the corresponding moments of the Boltzmann collision operator (5) vanish when collisions between particles of the same species are considered. Therefore the right hand-side term in (21) and (22) quantify only effects of collisions between different species, s and t.

When equations (20)–(22) are summed over all species s one obtains, after some algebra, the equations of conservation of mass, momentum and energy for the plasma (Montgomery and Tidman 1964, p. 198):

$$ \frac{\partial {\rho_m}}{\partial {t}} + \nabla \cdot (\rho_m {\bf U}) = 0 $$

(23)

$$ \rho_m \frac{D{\bf U}}{Dt} + \varvec{\nabla} \cdot \underline{\underline{{P}}} -\rho_m{\bf g} - \rho_c\left({\bf E} + {\bf J} \times {\bf B} \right) = 0 $$

(24)

$$ \frac{3}{2} nK\frac{DT}{Dt} - \frac{3}{2}kT\nabla \cdot \sum_s{n_s{\bf U}_s} +\varvec{\nabla} \cdot {\bf q} + \underline{\underline{{P}}}:\varvec{\nabla} {\bf U} -({\bf J} - \rho_c {\bf U})\cdot\left[ {\bf E} + {\bf U} \times {\bf B} \right] =0 $$

(25)

where ρ _m is the plasma mass density, U is the plasma bulk velocity, J is the net current density, $\underline{\underline{P}}$ is the plasma pressure tensor, q is the plasma heat flux (see Appendix 1 for the definitions of these plasma parameters).

Both systems of equations, (20)–(22) and (23)–(25) are not closed; they contain more unknowns than equations. Indeed the heat flux q _s , q and pressure tensor $\underline{\underline{{p_s}}}$, $\underline{\underline{{P}}}$ are still functions of the velocity distribution function. They cannot be determined from the first three transport equations given above, without making additional assumptions and approximations. In Appendix 2 we discuss two classical methods for closing this chain of equations (Chapman–Enskog, Grad). One possible simplifying assumption is to consider that the VDF of all species is an isotropic Maxwellian (in the frame of reference comoving with the velocity u _s  = U); in this case q _s  = 0 and the pressure tensor $\underline{\underline{{p_s}}}$ is isotropic. Under these assumptions the (20)–(22) become the one-fluid Euler or 5-moment equations:

$$ \frac{\partial {n}}{\partial {t}} + \varvec{\nabla}_{{\bf r}}\left(n{\bf U}\right) = 0 $$

(26)

$$ \frac{\partial {{\bf U}}}{\partial {t}} + {\bf U} \cdot \varvec{\nabla}_{{\bf r}} {\bf U} = \frac{{\bf F}} {m} - \frac{1}{mn}\varvec{\nabla}_{{\bf r}}\left( nkT \right) $$

(27)

$$ \left(\frac{\partial }{\partial {t}}+ {\bf U} \cdot \varvec{\nabla}_{{\bf r}} \right)\left( Tn^{-2/3} \right) = 0 $$

(28)

The first hydrodynamic models of the solar wind were based on the one-fluid Euler approximation that will be discussed in more details in Sects. 3.2 and 5.

The next approximation considers the electrons and protons as separate fluids with different temperatures and isotropic pressure tensors for each species, but the same bulk velocity, u _s  = U. (e.g. Hartle and Sturrock 1968). The corresponding two-fluid equations and their relationship to the Euler set are discussed in Sects. 5.3.1 and 5.3.2.

The general theoretical framework outlined in this section reveals the complementarity between kinetic and fluid approaches. The kinetic exospheric models are applications of the Vlasov and/or Fokker–Planck equations and are discussed in Sect. 4. Multi-fluid models may be also viewed as applications of the Boltzmann and/or Fokker–Planck equations, based on various expansions of the velocity distribution functions and corresponding to different approximations of the transport equations. The latter are reviewed in Sect. 5.

3 A Historical Account of Modeling the Solar Corona and the Solar Wind

During all solar eclipses a bright halo (the solar corona) can be seen around the Sun when it is masked by the Moon. The coronal luminosity extends then sometimes up to 10 solar radii (10 R _S ). This brightness in the visible spectrum is mainly due to scattering of the photospheric light by free electrons of the fully ionized coronal gas (Thomson scattering).

It has been established since 1942 (Edlén 1942) that this extended solar coronal plasma is heated to more than one million degrees. The high temperature was established from the discovery of the high degree of ionisation of the atoms emitting the yellow, green and red coronal emission lines. It was later confirmed by the broad Doppler width of these spectral lines. The heating mechanism is still not fully understood (see, e.g., Gomez 1990; Heyvaerts 1990; Hollweg 1991; Zirker 1993; Narain and Ulmschneider 1996).

The coronal plasma is constituted of fully ionized hydrogen, helium (5–10%), as well as smaller traces of highly ionized heavier atoms. Up to 1958, it was considered that the density and kinetic pressure of the coronal plasma is maintained in hydrostatic equilibrium by the Sun’s gravitational field. This scenario would correspond then to u _s  = 0 in equations (20)–(22).

Due to this high temperature ($T_e {\approx}10^6$ K), the plasma thermal conductivity is extremely large. Therefore the coronal temperature was assumed to be almost uniform, i.e., independent of heliographic altitude. The radial distribution of the coronal brightness was generally fitted to theoretical electron density distributions considered to be in isothermal hydrostatic equilibrium (van de Hulst 1953).

3.1 Chapman’s Hydrostatic Models of the Solar Corona

The first non-isothermal model of the corona was proposed in the late 1950s by Chapman (1957). Arguing that thermal conductivity of the coronal plasma is large, but not infinitely large and proportional to $T_e^{5/2}$, where T _e is the free electron temperature (equal to T _p , the proton temperature), Chapman (1957) obtained a distribution for T _e (r) as a function of r, the radial distance; he postulated that heat conduction is the only mechanism to transport heat away from the base of the corona into interplanetary space. By requiring that T _e  = 0 at infinity (i.e. tending to that of the interstellar medium), he determined that T _e (r) should decrease as r ^−2/7 from over T _e  = 10⁶ K at the base of the corona, to zero for r → ∞. Assuming T _e (r ₀), the coronal temperature at r ₀ = 1.06 R _S , is equal to 10⁶ K, the electron temperature at 1 AU would still have a high value (2 × 10⁵ K) in Chapman’s conductive model. This surprising result led him to argue that the terrestrial thermosphere could be heated from above by thermal conduction, i.e., by heat being conducted down from interplanetary space into the upper atmosphere of the Earth (Chapman 1957, 1959).²

Applying the theoretical r ^−2/7 temperature distribution Chapman (1957) calculated the distribution of the electron pressure (p _e  = n _e k T _e ) and density (n _e ) in the corona assuming it is in hydrostatic equilibrium. He found that his theoretical electron (and ion) density had a minimum value at 0.81 AU. Using a density of 2 × 10⁸electrons/cm³ at the base of the corona—obtained from eclipse observations—he found that the minimum value of the electron density should be equal to 340 cm⁻³ at r = 174 R _S ; beyond this heliocentric distance in Chapman’s hydrostatic model, n _e (r) increases indefinitely with altitude; at 1 AU, n _e  = n _p  = 342 cm⁻³ with the same density of protons and electrons due to the quasi-neutrality of the interplanetary plasma. The properties at the base of the solar corona and at 1 AU corresponding to Chapman’s conductive and hydrostatic model are listed in Table 1 (second column; model {0}).

Table 1 Slow (or quiet) solar wind observations and model at 1 AU and at a coronal reference level, r ₀

Such a theoretical density profile (i.e. a density increasing versus altitude) implies that the assumed hydrostatic equilibrium is convectively unstable. In one of his subsequent papers, Chapman (1961) argued that the interplanetary medium might become turbulent beyond 0.81 AU (174 R _S ) but he did not envisage the more realistic hypothesis that the corona might not be in global hydrostatic equilibrium, but expanding continuously as observed later on by in situ measurements in the interplanetary medium.

3.2 The First Hydrodynamic Model of the Solar Wind

Parker’s argument against hydrostatic equilibrium and in favor of a continuous expansion of the solar corona is based on a hydrodynamic description. Parker (1965, p. 669) argued that “if an atmosphere is sufficiently dense that the mean free path is small compared to the scale height and/or the radial distance r from the parent body, then the ordinary hydrodynamic equations are appropriate”. And according to him “this is the situation in the corona”. On the other hand “if the atmosphere is so tenuous that the mean free path is long, the universal presence of weak magnetic fields and attendant instabilities maintain approximate statistical isotropy in the thermal motions, and it turns out that again the large-scale, low frequency variations of the gas are describable, in at least an approximate way, by conventional hydrodynamic equations with isotropic pressure. This appears to be the situation in the expanding corona at the orbit of Earth” (Parker 1965, p. 669).

In his pioneering paper, Parker (1958b) demonstrated that hydrostatic models of the corona predict kinetic pressures and densities exceeding those existing in the interstellar medium. In his well known monograph, Parker (1963) estimated the interstellar density to be “1 particle/cm³ and the kinetic pressure to 1.4 × 10⁻¹⁴ dyne/cm²” (Parker 1963, p. 42). These values are much smaller than those predicted at infinity by hydrostatic equilibrium for any coronal temperature decreasing with r slower than 1/r, as it is the case for instance in Chapman’s conductive model ($T\approx r^{-2/7}$), or in an isothermal corona. When the high coronal temperature extends deep into the interplanetary medium, declining asymptotically slower than 1/r, “the atmosphere will extend to infinity with a finite pressure” which is much larger then that of the interstellar medium. “Since there is nothing at infinity to contain such a pressure, the atmosphere would expand, rather than being static” (Parker 1963, p. 43, par. 1).

This reasoning based on a mechanical pressure balance equilibrium argument led Parker (1958b) to propose that the plasma is driven out of the corona, gaining momentum as a result of a large kinetic pressure gradient: i.e. due to the large pressure imbalance at the base of hot corona and in the cold interstellar medium. In other words, the excessively large kinetic pressure due to the high coronal temperature generates a pressure gradient in the momentum equation, accelerating the coronal plasma from slow (subsonic) expansion to a supersonic flow regime beyond some critical altitude where the radial bulk speed becomes equal to the sound speed.

It must be emphasized that the essential part of the supersonic outflowing is the high temperature of the corona extending far out in the solar gravitation field due to the high thermal conductivity of the coronal plasma and ad hoc in situ heating. Maintaining the high temperature in the escaping gas is responsible for the ultimate supersonic velocity (Parker 2009, private communication).

Parker’s first solar wind model was based on the assumption that the coronal plasma is a fully ionized ideal gas, that the uniform temperature is the same for the electrons and protons, and that the radial expansion is not intermittent but stationary. The following one-dimensional equation was derived by Parker from the general Euler equations [see (26)–(28)]:

$$ \frac{r}{U}\frac{dU}{dr} = \frac{r^3}{\left(2KT\right)/m-U^2} \left[ \frac{d}{dr} \left(\frac{2KT}{mr^2} \right) + \left(\frac{GM_S}{r_0}\right) \frac{1}{r^4}\right] $$

(29)

where the bulk speed, U, and the temperature, T, depend on the radial distance, r; $\sqrt{\frac{GM_S}{r_0}}$ is the escape velocity at radial distance r ₀. The implicit assumptions for solving equation (29) are that (1) the corona (and solar wind) are isothermal; (2) the pressure tensors are diagonal and identical for the electrons and protons and (3) the average velocities of the electrons and protons are equal to U, the bulk (mass averaged) velocity of the plasma as a whole, implying no electric current since the plasma is quasi-neutral. The modeling of plasma flow was restricted to the equatorial region.

Furthermore, due to the quasi-neutrality of plasmas, the electron number (n _e ) density is necessarily equal to the total ionic charge densities ($\sum_{i \ne e} Z_i n_i = n_e$) and the net current must also be zero. Under these additional restrictive conditions, the different particle species are not expected to diffuse with respect to each other; the plasma can be modeled as a whole, almost like a neutral gas, with a unique bulk (mass averaged) speed U, a unique characteristic temperature T, and a mass density ρ. In Parker’s SW model, the plasma distribution is spherically symmetric, and U has only one component in the radial direction (U _r ).

Equation (29) has a family of solutions with one critical point r _c , defined by the solution of the equation (Parker 1965):

$$ r^3\frac{d}{dr}\left( \frac{2KT}{mr^2} \right) + \left(\frac{GM_S}{r_0r}\right) = 0 $$

(30)

The critical point is located somewhere between 4 R _S and 8 R _S . It is at this heliocentric distance that the bulk velocity (U _r ) changes from subsonic values to supersonic ones (see Fig. 1). The critical point is unique if T decreases with r less rapidly than 1/r (Parker 1965). Only the critical solution, passing through the critical point, achieves the supersonic acceleration of the solar wind; it is the only one that satisfies obvious boundary conditions at r ₀, the base of the corona (where the density and kinetic pressure are large, and bulk velocity small), and at infinity (where the mass density and kinetic pressure are small). The SW bulk velocity (v in Parker’s original paper or U in this work) varies from very small subsonic velocities at the base of the corona (branch “A” below the critical point in Fig. 1) to supersonic speeds of more than 300 km/s along the branch “B”. Typical subsonic solutions are also illustrated by the lines displayed in Fig. 1 below the critical point; in this family of hydrodynamic solutions, the expansion velocity increases from small values at the base of the corona (r = r ₀ = a), but u _r never reaches the velocity of sound; eventually u _r decreases to zero as r → ∞. In all subsonic solutions, the plasma density and kinetic pressure decrease asymptotically to constant values, which are much too large compared to the corresponding value in the interstellar medium, just as in the hydrostatic solutions (U = 0). The lines asymptotic to the branch “C” correspond to mathematical solutions of the stationary Euler differential equations but with no physical relevance. They correspond to shock like acceleration with density decreasing to zero at altitudes below that of the critical point. When U(r ₀), the bulk velocity at the reference level r ₀, exceeds the value corresponding to Parker’s critical solution, the hydrodynamic equations have no stationary solution extending all the way through to r → ∞; in these instances only time-dependent hydrodynamic solutions can exist over the whole range of altitudes.

Fig. 1

The figure shows the variation of the solar wind bulk velocity with the radial distance as inferred from several classes of solutions for the Euler equation (27) in the form (29) derived by Parker. The position of the critical point is derived from (30); the critical solution, AB, is the only one that gives the supersonic acceleration of the solar wind (Parker 1965)

The existence of the supersonic plasma flow described by Parker’s hydrodynamic solar wind model is supported by Biermann’s (1953) observations of comets and also by in situ observations of Mariner 2 (Neugebauer and Snyder 1962; Snyder and Neugebauer 1964). Missions like IMP 1 & 2, Vela 2 & 3 confirmed the permanence of the supersonic solar wind. Continuous in situ observations reveal that the solar wind plasma is not stationary but that its average bulk velocity changes irregularly between 300 km/s and 450 km/s in what was called by Hundhausen (1968, 1970) the “slow” (or “quiet”) solar wind streams. By increasing the coronal temperature from 1.5 × 10⁶ K to over 3 × 10⁶ K, Parker’s hydrodynamic solutions are able to account for such a range of proton bulk velocities. Later on, however, it has been recognized that a coronal temperature exceeding 3 × 10⁶ K would be needed for Parker’s hydrodynamic model to account for bulk velocities of 600–900 km/s (see Fig. 6.1 in Parker 1963, p. 75), which are sometimes observed at 1 AU in fast speed streams of the solar wind.

Parker’s pioneering theoretical model of the interplanetary plasma inspired other subsequent hydrodynamic models of the solar wind. It influenced the theoretical modeling of other plasma flows such as plasma escaping out of the topside terrestrial ionosphere above the polar caps, along open magnetotail field lines (Banks and Holzer 1968, 1969), that is the polar wind. The history of the polar wind models and their development parallel to that of the solar wind models has recently been reviewed by Lemaire et al. (2007).

3.3 Convective Instability in Chapman’s Hydrostatic Corona

Lemaire (1968) showed that the temperature gradient in Chapman’s conductive models, dT _e /dr, becomes super-adiabatic at 0.22 AU (44 R _S ), i.e. steeper than the adiabatic temperature lapse rate (see also Lemaire 2010). Hence, Chapman’s hydrostatic model is necessarily convectively unstable at an altitude below the minimum of the density in the hydrostatic model of the solar corona. This means that the mechanism of heat conduction alone is not efficient to carry away all the energy deposited at the base of the corona, and transport it out into interplanetary space.

Furthermore, Lemaire (1968), showed that, beyond the radial distance of 0.22 AU (44 R _S ), turbulent convection (within a plasma in average hydrostatic equilibrium, as is the case in the Sun’s convection zone) is also insufficient to evacuate the coronal energy flux toward outer space. Lemaire concluded that a continuous expansion of the coronal plasma (i.e. a laminar hydrodynamic advection/expansion) is required to keep the actual temperature gradient smaller than the adiabatic temperature lapse rate at all altitudes in the corona and in the interplanetary medium.³

This new inference led Lemaire (1968) to conclude that Parker’s hydrodynamic expansion of the solar corona is needed from a thermodynamic point of view. He argued that only a steady state “explosion” is able to evacuate the energy out of the corona. Neither heat conduction, nor the Böhm-Vitense (1954) turbulent convection (i.e. hot elements raising and cooled ones falling within an atmosphere in average hydrostatic equilibrium as in the solar convection region) proved to be efficient enough.

It can be seen that Lemaire’s (1968) argument is different from but complementary to the mechanical one proposed by Parker (1958b) to support the existence of a radial supersonic expansion of the solar corona. Parker’s theoretical argument in favor of a continuous radial expansion of the corona is based on pressure imbalance conditions and the boundary conditions for the momentum transport equation, as was pointed out in the previous section.

3.4 A Brief Historical Perspective on the Hydrodynamic Models of the Solar Wind

In a monograph (Parker 1963) and in subsequent publications, Parker (1965, 1967, 1969) generalized his first isothermal hydrodynamic model of the solar wind. He replaced the hypothesis of a uniform coronal temperature by the assumption that the temperature decreases with radial distances according to a polytropic relationship between the plasma temperature and its density.

It would be too long to review all the alternative hydrodynamic SW models that flourished during four decades. Following the early theoretical generalization of the first SW model, more sophisticated hydrodynamic models have been formulated. Various formulations and approximations for the energy transport equation were coupled to the hydrodynamic continuity and momentum equations (23)–(24), to compute steady state radial distributions of the solar wind temperature as function of the radial distance (Noble and Scarf 1963; Parker 1964; Whang and Chang 1965; Weber 1970; Cuperman and Harten 1970a; Durney 1971, 1972). In these latter one-fluid models, the effect of thermal conductivity was taken into account to determine the radial distribution of the plasma temperature, as Chapman did for the conductive model discussed above. In hydrodynamic models heat transport by conduction is the dominant mechanism close to the base of the corona, while at radial distances close to and beyond the critical sonic point the energy flux is predominantly carried outwards by advection, i.e., by the solar wind bulk motion. The energy is then primarily carried away by the supersonic expansion (or “stationary explosion”) of the coronal plasma.

The rather small effects of viscosity (due to Coulomb collisions) have been considered in Navier-Stokes approximations of the hydrodynamic momentum equations (Scarf and Noble 1965; Whang et al. 1966; Konyukov 1969; Eisler 1969; Dahlberg 1970). Time-dependent hydrodynamic models of the solar wind expansion have also been developed. The first non-stationary model of the solar corona was a self-similar, polytropic radial expansion by Lemaire (1966a).

Instead of single-fluid models of the solar wind, two-fluid and multi-fluid models have been developed. The effects of non-radial magnetic fields, of adiabatic cooling for both electrons and protons, the acceleration by Alfvén waves, electron and proton heating have been investigated by a number of authors (Sturrock and Hartle 1966; Weber and Davis Jr., 1967; Hartle and Sturrock 1968; Urch 1969; Cuperman and Harten 1970b; Cuperman and Harten 1971; Hartle and Barnes 1970; Barnes and Hartle 1971; Whang 1971b; Wolff et al. 1971; Hansteen and Leer 1995; Esser and Habbal 1995; Habbal et al. 1995a; Tu and Marsch 1997; Tu and Marsch 2001; Tam and Chang 1999; Lyngdal Olsen and Leer 1999; Li 1999; Lie-Svendsen et al. 2001; Kim et al. 2004; Janse et al. 2006).

These SW models are listed in the successive columns of Tables 1 and 2 in a more or less chronological order. The solar wind parameters at 1 AU as well as the boundary conditions at the base of the corona for all these different types of stationary solar wind models are compiled in Table 1 for the slow/quiet solar wind, and in Table 2 for the fast solar wind. The average solar wind parameters observed at 1 AU are given in italics. The synthetic characteristics of each model are summarized: i.e. the one-fluid hydrodynamic (1F-H) models, two-fluid hydrodynamic (2F-H) models, three-fluid hydrodynamic (3F-H), three component hybrid/semi-kinetic (3Hb), two fluid sixteen moment (2F-16M) models, exospheric (E) models.

Table 2 Fast solar wind observations and model at 1 AU and at a coronal reference level, r ₀

Coupled multi-fluid models are more difficult to integrate numerically. Furthermore, they rely on an increased number of free parameters and boundary conditions to be imposed to the set of differential equations; this leads to a wider variety of solutions and therefore of a greater chance to fit the SW measurements at 1 AU.

As a consequence of the nonlinearity of the set of transport equations, relatively small changes of the boundary conditions in the corona can produce large amplitude changes in the values of the higher order moments (parallel and perpendicular temperatures, stress tensors, energy and heat fluxes) at 1 AU. Note that in situ heating of the SW plasma influences the terminal flow speed of the wind at 1 AU. The temperature anisotropy at 1 AU is also influenced heavily by in situ heating or momentum transfer mechanisms (i.e. heat deposition, pitch-angle scaterring, wave-particle interactions, momentum transfer by MHD waves).

Multi-fluid models consider more than one critical point (where the bulk speed becomes transonic and where dU/dr may have two different values as illustrated in Fig. 1). The existence of these additional mathematical singularities makes the search of relevant solutions much more difficult and uncertain, if not questionable. Indeed, when the mathematical singularities are located beyond the exobase, where the Knudsen number becomes larger than unity, the physical significance of these singularities becomes uncertain.

The effect of non-radial (non-spherically symmetric) expansions has been extensively modeled. Several radial distributions for the extended coronal heating source and accelerating mechanism have been proposed with the hope to reach hydrodynamic models with very high SW bulk velocity, as observed in fast speed streams (Whang 1971a; Hartle and Barnes 1970; Hansteen and Leer 1995). Since some of these hydrodynamic multi-fluid solar wind models assume that the kinetic pressure tensors are isotropic, these models cannot predict the observed anisotropies of the electrons and protons temperatures at 1 AU.

One-fluid models predicting the same temperatures for the protons and electrons are in principle less adequate than two-fluid models where the boundary conditions and source terms have been adjusted so that the electron temperature at 1 AU is higher than the proton temperature, as consistently observed in the slow solar wind. Two-fluid models have been developed by Hartle and Sturrock (1968) accounting for different temperature profiles of electrons and protons; the model by Leer and Axford (1972) allows an anisotropic proton temperature producing more reasonable temperatures at the Earth’s orbit. Other two-fluid models have been developed by Cuperman and Harten (1970a), Hartle and Barnes (1970), Cuperman and Harten (1971), Habbal et al. (1995a), Tu and Marsch (1997), 2001), Esser and Habbal (1995), Kim et al. (2004). The gyrotropic two-fluid model of Li (1999) and Janse et al. (2006) based on the 16-moment transport equations with an “improved treatment” of heat conduction and/or proton heating (see Sect. 5.4) achieves a rather satisfactory fit to fast SW observations. Some data about these two models, including SW properties at 1 AU, are presented in Table 2 and Fig. 4b, see models {26} and {28}.

The results predicted at 1 AU by the most representative of these models are synthetically reported in Tables 1 and 2 respectively for the slow/quiet solar wind and fast speed stream observations. Various aspects of these solar wind models have already been discussed in previous reviews by Marsch (1994, 2006).

In order to reproduce the high speeds of the fast solar wind, the fluid models must either assume a very high proton temperature in the corona, of the order of 10⁷ K, or extended heating beyond the critical point. Simple energy conservation arguments can be used to show that this is not a shortcoming of the fluid models, but a necessary requirement for any solar wind model, kinetic or fluid, unless an extremely large outward heat flux can somehow be formed in the corona.

Furthermore, in situ observations indicate that fast solar wind speed streams do not originate in the coronal region where the temperature is highest, but on the contrary, from polar coronal holes where the coronal electron temperature is significantly lower than in the equatorial regions, the source of the slow SW. In the “old” hydrodynamic solar wind models, which were driven mainly by electron heating, this was a serious problem as the low electron temperature should produce a slower solar wind from coronal holes. However, more recent observations of polar coronal holes by the UVCS instrument on the SOHO satellite (Kohl et al. 1997) have resolved this problem. They showed that hydrogen, and even more so heavier ions, are much hotter than electrons in coronal holes. This would imply that electrons play a smaller role in the fast solar wind acceleration and that it is driven mostly by proton heating.

As early as 1960, it was claimed that beyond a certain heliospheric distance, the coronal plasma is becoming almost collisionless. As a consequence, it was expected that the VDF can strongly deviate from a Maxwellian. Hence higher-order moments of the VDF can become large since they depend critically on the departure from a drifting Maxwellian. As shown in Sect. 5, such higher order fluid models can become laborious. The diffusion, viscosity and heat conductivity coefficients used in the hydrodynamic equations deviate from the standard expressions derived for a simple Maxwellian VDF, f ⁽⁰⁾(v, r, t). Within the Chapman–Enskog or Grad theories of non-uniform gases, f ⁽⁰⁾ is a zero-order approximation of the actual VDF, f(v, r, t). Furthermore, in the energy transport equation the divergence of the third order, Q _ijl , and fourth order moments, R _ijlk , (see definitions in Appendix 1) have been truncated in order to close the hierarchy of moment equations. Since such a truncation is dictated merely by the convenience to limit the mathematical complexity of the formulation, it lacks a sound physical argument. These particular methods of truncating the moment equations are not identical in the Chapman–Enskog theory (Chapman 1918; Enskog 1917) and in the one developed by Grad (1958). The main purpose of these ad hoc approximations was to obtain a description of the gas as close as possible to the transport equations used in classical hydrodynamics, exclusively designed to model transport in collision-dominated flows. Another more pragmatic reason to use such truncated representation is dictated by the complexity of the hydrodynamic transport equations and of their numerical solutions.

4 Kinetic Exospheric Modeling of the Solar Wind

Two years after the publication of Parker’s first solar wind model, Chamberlain (1960) proposed an alternative kinetic theory for the coronal expansion. Thus began a standing controversy between the proponents of hydrodynamic solar wind models and those arguing in favor of kinetic models.

4.1 On the Existence of a Collisionless Region

Chamberlain argued that beyond a heliospheric distance of 2.5 R _S in the solar corona, the Coulomb collision mean free path of thermal protons becomes larger than H, the atmospheric density scale height of the coronal plasma; H is the characteristic range of altitudes over which the density decreases by a factor e = 2.71. Chamberlain’s argument infers that K _n , the Knudsen number of plasma particles:

$$ K_n=\frac{\lambda_c}{H} $$

(31)

the ratio between λ _c , the mean free path of the particles and H, becomes larger than unity for r > 2.5 R _S . The density scale height in an atmosphere is defined by:

$$ H = {\frac{{kT}}{{\left\langle m \right\rangle g}}} $$

(32)

where g is the gravitational acceleration, and 〈m〉 is the mean atomic mass of the gas. In a fully ionized hydrogen plasma, like the solar corona, 〈m〉  = m _p /2.

When K _n  > 1 the VDF deviates from a Maxwellian VDF, the Chapman–Enskog and Grad expansions of the velocity distribution function (see Appendix 2) can then become inadequate. It is often considered that when K _n  > 1 higher-order fluid models are needed to capture the non-Maxwellian character of the velocity distribution function. However, on top of their complexity reported in Sect. 5.4, it is not a priori obvious how the higher-order models should be closed. In neutral classical gas hydrodynamics, it is considered that a fluid description is inadequate when K _n  > 1; indeed, there is no unique way based on sound physical arguments to close the system of moment equations.⁴

Based on SW models and observations, Hundhausen (1968) estimates that K _n  > 1 above 7 R _S and pointed out that, in order to explain the significant temperature anisotropy observed at 1 AU, the solar wind protons and heavier ions should be collisionless beyond 15 R _S . Brasseur and Lemaire (1977) checked that in Parker’s isothermal SW hydrodynamic model the Knudsen number becomes already larger than unity at r = 4 R _S ; this radial distance is lower than the sonic critical point which is located at r _c  = 6 R _S in this hydrodynamic SW model.

Chamberlain (1960) developed an exospheric model for the coronal ion-exosphere. His model which is known as the “solar breeze model” is based on Jeans (1923) theory for planetary exospheres. Depending on their velocity and of their angular momentum four classes of particles can be identified above the exobase:

• the escaping particles, which have sufficiently large velocities to escape from the Sun’s gravitational potential well,

• the ballistic or captive particles that fall back into the corona,

• the satellite or trapped particles which are continuously bouncing up and down between a magnetic mirror point and a gravitational turning point,

• the incoming particles arriving from the interplanetary regions and penetrating into the collision-dominated region below the exobase.

Taking into account the conservation of the total energy and of the magnetic moment $\mu_s=mv^2_{\bot}/(e_sB)$, and using the Liouville theorem the VDF of collisionless electrons and protons can be determined everywhere in the exosphere provided this function is specified at the exobase. When the VDF at the exobase is an analytical function of v, the VDF at any point in the exosphere is then also an analytical function of v. Furthermore the expressions of all the moments of the VDF are analytical expressions of the constant of motions (and implicitly of the electromagnetic and gravitational potentials) and can be calculated numerically everywhere in the exosphere (Lemaire and Scherer 1971).

Öpik and Singer (1959), as well as Brandt and Chamberlain (1960), applied Jeans’ collisionless (or zero-order) kinetic theory to calculate the density distributions of neutral gases in the terrestrial exosphere. Chamberlain (1960) applied a similar exospheric model for the ion-exosphere of the solar corona where the protons move without collisions along radial magnetic field lines in the gravitational field of the Sun and a polarization electric field. The latter is needed to maintain quasi-neutrality in planetary and stellar ionospheres in hydrostatic equilibrium.

4.2 The Pannekoek–Rosseland Electrostatic Field and Potential Distributions

The Pannekoek–Rosseland (PR) electrostatic field is induced in plasmas by the tendency of gravitational force to polarize the plasma by separating the ions and the less massive electrons. The presence of this charge separation E-field in a plasma was first introduced by Pannekoek (1922) and Rosseland (1924) in their studies of ionized stellar atmospheres. The gravitational force acting on the massive ions is much larger than that acting on electrons. Thus, like in dielectric material, the minute charge separation induced by the gravitational forces polarizes the plasma and generates an E-field which is oppositely directed to the gravitational acceleration, g. This polarization electric field prevents the diffusion in the gravitational field of positive charges with respect to the negative ones, due to their mass difference.

When the electrons and protons are in hydrostatic equilibrium in the gravitational field, and when their temperatures are equal and independent of altitudes (isothermal), the PR electrostatic field is determined by:

$$ {\bf E}_{PR}= - \nabla \Upphi_{PR}= -\frac{\nabla p_{e}}{e n_e} = - \frac{(m_p - m_e) {\bf g}} {2e} = \nabla \left[\frac{(m_p- m_e) \Upphi_g}{2e}\right] $$

(33)

where e is the electron charge; p _e  = n _e k T _e is the electron kinetic pressure; m _p and m _e are the proton and electron masses respectively; $\Upphi_g$ and $\Upphi_{PR}$ are respectively the gravitational and the electrostatic potential. These potentials are functions of the altitude in the corona. Equation (33) determines the PR electric field and the electric potential in a pure hydrogen plasma. In the solar corona $\Updelta \Upphi_{PR}$, the electrostatic potential difference between Chamberlain’s exobase (r = 2.5 R _S ) and r = ∞, is equal to 150 V.

In the case of multi-ionic plasmas, similar analytical relationships can be derived for the electrostatic potential, $\Upphi_{PR}$, even when the ion temperatures (T _i ) and the electron temperature (T _e ) are not equal (cf. Sect. 5.2.2 in Lemaire and Gringauz 1998). The most general expressions for this polarization electric field, including the effect of field-aligned flows of ions and electrons can be found in the paper by Ganguli and Palmadesso (1987). An expression for the parallel electric field sustained by two-dimensional cross-B sheared plasma flows has been derived by Echim and Lemaire (2005).

The density of the polarization charges responsible for the PR electrostatic field is extremely small: (n _p  − n _e )/n _e  = 4 × 10⁻³⁷ (Lemaire and Scherer 1970). This indicates that the quasi-neutrality equation (n _p  = n _e ) is indeed a good approximation in most space plasmas; therefore, instead of solving Poisson equation to calculate the electrostatic field distribution, it is rather satisfactory (and much simpler from a computational point of view) to solve the quasi-neutrality equation by an interative method in order to derive the spatial distribution of the electrostatic potential, $\Upphi_E(r)$.

When (33) is applicable, the sum of the gravitational potential energy and electric potential energy, R(r), is the same for the protons and for the electrons:

$$ R_p(r) = m_p \Upphi_g + e \Upphi_{PR} = m_e \Upphi_g - e \Upphi_{PR} = R_e(r) = (m_p + m_e) \Upphi_g(r) $$

(34)

As a consequence of this equality, H _p , the density scale height of the protons, is precisely equal to H _e , the density scale height of the electrons (see Sect. 5.2.2 in Lemaire and Gringauz 1998). This has to be so, since the quasi-neutrality condition has to be satisfied at all altitudes.

Another important consequence of the presence of the PR electric field is that the minimum energy for an electron to escape out of the potential well (i.e. the critical escape energy, $0.5m_e {v^2_{e\infty}}$ at any given altitude), is precisely equal to that of a proton at the same altitude:

$$ \frac{1}{2} m_p {v^2_{p\infty}} = - R_p(r_0) = \frac{1}{2} (m_p+ m_e) g_0 r_0 = \frac{1}{2} m_e {v^2_{e\infty}} = - R_e(r_0) $$

(35)

Using Chamberlain’s notation, $U_p = (2 k T / m_p)^{1/2}$ corresponds to the proton thermal velocity, and $\lambda_{p0} = m_p {v_{p\infty}}^2 / m_p {U_p}^2$ is the dimensionless ratio between the critical escape energy and the thermal energy of a proton at the exobase altitude. A similar dimensionless parameter can be defined for the electrons: $\lambda_{e0} = m_e {v_{e\infty}}^2 / m_e {U_e}^2$.

4.3 The First Generation Kinetic/Exospheric Models for the Coronal Evaporation

Assuming that the proton VDF at the exobase of the collisionless region can be approximated by a truncated Maxwellian, the Jeans evaporative flux of particles with velocities exceeding the critical escape velocity at the exobase is given by:

$$ F_p(r_0) = n_p(r_0) U_p \left[ 1 + \lambda_{p0}\right] \frac{ e^{-\lambda_{p0}}}{ 2 \sqrt{\pi}} $$

(36)

Chamberlain (1960) also derived analytical formulas for n _p (r), the density of all particles populating the region above the exobase (i.e. the sum of the density of the escaping particles, of the ballistic or captive ones, as well as trapped or satellite particles). From these analytical expressions (which will not be given here), he was able to calculate the radial distribution of n _p (r) and of u _r (r) = F _p /n _p , the average proton velocity:

$$ u_r(r) = (U_p/ 2 \pi^{1/2} ) \left[ 1 + \lambda_{p0} \right] e^{ -\lambda_p} \left(\lambda_p/\lambda_{p0}\right)^2 $$

(37)

where

$$ \lambda_p(r) =\lambda_{p0} \frac{r_0}{r} = 2 \frac{R_p(r)}{m_p {U_p}^2} $$

(38)

is the ratio of the total potential energy and average thermal energy of the protons at the radial distance r. A similar function of r can be defined for the electrons:

$$ \lambda_e(r) =\lambda_{e0} \frac{r_0}{r} = 2\frac{R_e(r)}{m_e {U_e}^2} $$

(39)

As a consequence of (35) and the definition of U _i , the functions λ _e (r) and λ _p (r) have the same value everywhere in the exosphere when the exobase temperature of the electrons is equal to that of the protons (T _e  = T _p  = T).

From (37) it follows that the value of the expansion velocity of the protons is subsonic (u _r  < U _p ) everywhere above the exobase. This is why Chamberlain called his model “solar breeze,” in contrast to the supersonic bulk velocity inferred by Parker in his solar wind model. In the exospheric solar breeze model u _r (r) decreases monotonically with r, from a maximum value, ($U_p/ 2 \pi^{1/2} \left[1 + \lambda_{p0}\right] e^{-\lambda_{0p}}$, at the exobase, to zero at r = ∞.

For an exobase proton temperature of T _p  = 2 × 10⁶ K and an exobase density of n _p (r ₀) = 10⁶ cm⁻³, one has λ_p0 = 2.3, and the average expansion velocity at 1 AU is 20 km/s. Chamberlain’s exospheric calculations predict also a proton density of 370 cm⁻³ at 1 AU, which he considered then to be in fair agreement with values derived by Behr and Siedentopf (1953) from zodiacal light measurements.

It can be seen that the value of plasma bulk velocity (20 km/s) predicted by the solar breeze model is more than an order of magnitude smaller than that predicted by Parker’s hydrodynamic solar wind model at the orbit of the Earth, i.e., 500 km/s for a comparable coronal temperature (see models {7} and {10} in Table 1). The bulk velocity predicted by the exospheric solar breeze solution is also much too small compared to those which were measured later on by the first interplanetary missions (300–700 km/s) at 1 AU. Furthermore, the proton density and temperature predicted at 1 AU by Chamberlain’s evaporative model were quite different from the measurements made in interplanetary space in 1962 and afterwards. This led the space physics community to disregard the solar breeze model, and by the same token, all other kinetic or exospheric models of the solar corona like those discussed later on by Jenssen (1963), Brandt and Cassinelli (1966). In these alternative exospheric models the ballistic/captive and trapped/satellite protons are missing; only the escaping protons contribute to the exospheric density in Jenssen (1963) or Brandt and Cassinelli (1966) evaporative models. This rather questionable assumption lead to smaller exospheric proton densities than in the solar breeze model and therefore to higher values for u _r (respectively 290 and 266 km/s at 1 AU). Nevertheless, these new exospheric models for the evaporation of the solar corona were also based on the postulate that the polarization electric field is given by the Pannekoek-Rosseland distribution. As a matter of consequence they were unable to predict bulk velocities comparable to those observed in the solar wind at 1 AU.

This failure restricted the popularity of kinetic exospheric SW models and fed the persistent and overwhelming belief that only hydrodynamic solutions could offer satisfactory descriptions for the supersonic coronal expansion. It was generally claimed that the exospheric solutions do not satisfy the plasma transport equations. The latter belief is clearly incorrect since all the moments of the VDFs obtained by the kinetic exospheric models necessarily satisfy the entire set of moment equations, including (20)–(22). Despite the undeniable evidence that beyond 4 R _S or so, the Knudsen number of the proton and electron gases becomes larger than unity, kinetic approaches were generally dismissed and merely considered as “academic exercises” of no or little relevance in modeling the solar wind.

At about the same epoch Lemaire (1968) had studied the formal asymptotic behavior (for r → ∞) of the radial distributions of the density, bulk velocity, and temperature in Chamberlain’s evaporative solar breeze model, as well as in Jenssen (1963) alternative ion-exosphere model. He compared these two different types of exospheric distributions with Parker’s hydrodynamic solar wind models. This comparative study did not explain, however, why the first generation kinetic models of the solar wind predict much lower expansion velocities than the hydrodynamic models existent at that time. The origin for this basic disagreement was discovered later on by Lemaire and Scherer (1969, 1971). An exospheric model giving a supersonic terminal velocity was found independently by Jockers (1970). This is recalled in detail in the following section.

4.4 The Second Generation of Kinetic/Exospheric Models

4.4.1 Modification of the Pannekoek-Rosseland Electric Field Distribution by the Zero Flux Condition

One common feature of the first generation exospheric models discussed in previous paragraphs was the Pannekoek-Rosseland electric field given by (33). It should be reminded that this electrostatic potential distribution had originally been introduced to model stellar ionospheres in hydrostatic equilibrium, like Chapman’s conductive model of the solar corona.

But we presently know that the coronal plasma is not in hydrostatic equilibrium, and that protons and electrons continuously evaporate from the solar corona. Since the basic hypothesis used to derive (33) fails to be satisfied, the polarization electric field keeping the coronal plasma quasi-neutral ought to be different from the Pannekoek-Rosseland (PR) electric field.

Lemaire and Scherer (1969) showed that the PR electric field (33) is inadequate because, according to (36), the evaporative flux of the protons would be 43 times smaller than the escape flux of the electrons; this results from the much larger electron thermal speed. Indeed, when T _e  = T _p and n _e (r ₀) = n _p (r ₀), (38) and (39) give λ_e0 = λ_p0 and λ _e (r) = λ _p (r), and according to (36),

$$ F_e(r) / F_p(r) = U_e/U_p = (m_p / m_e)^{1/2}=43. $$

(40)

This situation is not sustainable since it would imply a huge radial electric current flowing out of the solar corona, thus an unphysical electric charging of the Sun. Despite that Pikel’ner (1950) discovered the same inconsistency already two decades ago, this was not generally known to the community until 1969 when rediscovered independently by Lemaire and Scherer. Jockers (1970) developed independently a kinetic exospheric solar wind model quoting the earlier finding by Pikel’ner (1950).

In order to balance the net electric current, a larger potential difference, $\Updelta\Upphi_E$, is required between the exobase and infinity than that corresponding to the PR electrostatic potential distribution, $\Updelta\Upphi_{PR}$. Considering that the exobase is located at a radial distance of 6.6 R _S , that the exobase temperature of electrons is 1.4 × 10⁶ K and 0.9 × 10⁶ K for the protons (these values were inferred from Pottasch’s observations, 1960), Lemaire and Scherer (1971) calculate that $\Updelta\Upphi_{E}$ had to be 410 V, in order to balance the net escape of the electrons and protons, and therefore to keep the radial component of the electric current equal to zero. It can be seen that in this second generation of exospheric models the total field-aligned electric potential drop is more than two times larger than that of the PR field implicitly inferred in Chamberlain’s 1960 first generation exospheric solar breeze model (Lemaire and Scherer 1971).

The much larger value of $\Updelta\Upphi_E$ compared to $\Updelta\Upphi_{PR}$ accelerates the coronal protons to significantly larger velocities at 1 AU than in the first generation exospheric models, which were based on $\Updelta\Upphi_{PR}$ electric potential. For the exobase conditions listed above, Lemaire and Scherer (1971) obtained at r = 215 R _S the following values for the bulk velocity of the protons and electrons: u _r  = 320 km/s and n _p  = n _e  = 7.2 cm⁻³ for their densities (see model {3} in Table 1). Both set of values were in better agreement with the average quiet time solar wind observations reported by Hundhausen (1970) (column {1} in Table 1) than the values predicted by earlier exospheric and hydrodynamic models listed in Table 1 (see also Table 2 in Lemaire and Scherer 1971).

For completeness, it should be added that the actual values of the electric potential distribution, $\Upphi_E(r)$, at all altitudes in the exosphere, are uniquely determined by solving the quasi-neutrality equation:

$$ n_e\left[r, \Upphi_g(r), \Upphi_E(r)\right] = n_p\left[r, \Upphi_g(r), \Upphi_E(r)\right] $$

(41)

The exospheric densities of the protons and of the electrons in (41) are analytical functions of $\Upphi_E$, which are determined by the exospheric theory. Therefore, the radial distribution of $\Upphi_E(r)$ can be calculated for any radial distance, r, in the ion-exosphere. An iterative method had first been used to determined $\Upphi_E(r)$ by Lemaire and Scherer (1969, 1971); a direct integration method was proposed later on by Lemaire et al. (1991).

Figure 2 shows the ratio of the electric force (eE) and gravitational force (m _p g) acting on protons in the exospheric model (LSd) of Lemaire and Scherer (1973). The dots are values for |eE/m _p g| deduced from an empirical electron density distribution derived by Pottasch (1960) from eclipse observations. This comparison shows that the values of |eE/m _p g|, of the exospheric models increase in the inner corona from 0.5 (the PR value shown by the dashed horizontal line), to a maximum value larger than 1 at r = 50 R _S . Note that beyond the exobase level in this kinetic model, the electric force which accelerates the protons to supersonic speed is larger than the gravitational force.

Fig. 2

Variation with radial distance of the ratio between the electric and gravitational force on a proton deduced from an empirical coronal density distribution (dots); the solid line is obtained by the kinetic model of Lemaire and Scherer (1973)

4.4.2 Revisited Second Generation Kinetic/Exospheric Models: Comparison with Previous Models and Data

The radial distribution of the plasma density in the exospheric model published by Lemaire and Scherer (1972) is illustrated by the solid line in the top panel of Fig. 3. The dashed line indicates the r ⁻² density profile corresponding to the asymptotic variation of the solar wind density at large distances where the expansion velocity is almost independent of r. The empirical densities determined by Pottasch (1960) from eclipse observations are displayed by square symbols up to r = 16 R _S , while the range of solar wind measurements at the Earth’s orbit is indicated by a vertical bar at r = 215 R _S .

Fig. 3

Exospheric distribution of the density (n _p(e)), bulk velocity (W _p(e)), electron ($T_{e\bot}$) and proton ($T_{p\bot}$) perpendicular temperatures and omnidirectional temperatures (〈T _e 〉, 〈T _p 〉) of the LSc model (adapted from Lemaire and Scherer 1972)

The second panel in Fig. 3 shows the flow velocity in this exospheric solar wind model. The value of u _r (labeled W _p(e)) increases from 160 km/s at the exobase (r ₀ = 6.5 R _S ), to 307 km/s at 1 AU. This flow velocity is supersonic in almost the whole ion-exosphere. The range of bulk speeds measured in the quiet solar wind regimes are indicated by a vertical bar at r = 215 R _S (Hundhausen 1970).

It can be seen that the theoretical predictions of these new kinetic models fit well the average measurements of n and u _r at 1 AU. These values are reported in columns {22} and {1} of Table 1. The pitch angle averaged electron and proton temperatures, 〈T _e 〉 and 〈T _p 〉, predicted by this kinetic model fall within the range of values measured at 1 AU and reported by Hundhausen (1970) for the quiet solar wind. This is also indicated in Table 1.

The range and mean values of the corresponding parameters for the quiet solar wind reported by Hundhausen (1970) are given in column {1}. The agreement with observations of the revisited exospheric models is definitely much more satisfactory than for the solar breeze model based on the Pannekoek-Rosseland model of the electric field.

Thus, the results of exospheric models LSb and LSI included in Tables 1 and 2 abrogate Chamberlain’s claim “that Parker’s supersonic solution is not supported by some simple evaporative explanation.” Indeed, if a correct polarization electric field distribution would have been used by Chamberlain, he would have obtained results which are in surprisingly good agreement with the observations at 1 AU. The data in Tables 1 and 2 contradict a persistent but false consensus that “exospheric” models would only be “academic exercises” of little relevance for the physics of the solar wind. The predictions of the revisited exospheric model of Lemaire and Scherer (1971, 1972) were even in a better agreement with the observations than the predictions of hydrodynamic models available at that time.

From the second order moments of the exospheric VDFs, the pitch angle averaged temperature, $T =(T_{\parallel}+2 T_\perp)/3$, as well as the temperature anisotropies, $T_\parallel / T_\perp$, of the protons and electrons can be calculated. In Lemaire and Scherer (1971, 1972)—model LSb—the averaged proton temperature is 4.8 × 10⁴ K at the orbit of Earth, and the electron temperature 11.7 × 10⁴ K. As indicated above, these theoretical values fit well the quiet solar wind observations: T _p  = (4.4 ± 1.8) × 10⁴ K, T _e  = (14 ± 5) × 10⁴ K (Hundhausen 1970).

The hybrid or semi-kinetic model of Jockers (1970), denoted {4} in Table 1, gives a less satisfactory agreement at 1 AU: u _r  = 288 km/s, n _p  = n _e  = 12 cm⁻³, T _p  = 6.7 × 10⁴ K, T _e  = 46 × 10⁴ K compared to column {1} of Table 1. Similar lack of agreement exists for the hybrid/semi-kinetic models {5} of Hollweg (1970), {6} Chen et al. (1972), and for the exospheric models {8} and {9} of Jenssen (1963), and Brandt and Cassinelli (1966) not reported in Table 1 for lack of space.

Non-Maxwellian electron distribution functions have been observed in fast and slow solar wind, for a large range of solar conditions, from quiet to active. Three main electron components have been identified: (i) the “core”, i.e. the thermal part of the electron VDF, (ii) the “halo”, i.e. the suprathermal part of the electron VDF and (iii) the “strahl”, i.e. an antisunward magnetic-field-aligned population (Feldman et al. 1975; Pilipp et al. 1987b). Maksimovic et al. (2005) and Štverák et al. (2009) have studied the radial evolution of the relative contribution of each of the three electron populations. The model of Štverák et al. (2009) based on data from Helios, Cluster and Ulysses suggests that the relative density of the “strahl” decreases while the relative density of the “halo” increases with increasing radial distance, for both fast and slow solar wind regimes. It was also shown that in the fast SW the nonthermal tails of the electron VDF form close to the Sun; this is a confirmation of the assumptions made on the boundary conditions of the third generation exospheric models.

The asymptotic total electron temperature at large distances has been computed from second gerneration exospheric models by Meyer-Vernet and Issautier (1998). They established that the average electron temperature radial profile may be approximated by the sum of a term proportional to r ^−4/3 plus a constant term, with both terms being of the same order of magnitude at about 1 AU. They showed also that this remarkable asymptotic behavior is independent of the velocity distribution function adopted at the exobase altitude. This typical behavior seems to be in agreement with the results obtained for the temperature of the “core” electrons determined from a statistical study by Issautier et al. (1998) based on Ulysses observations. Meyer-Vernet and Issautier (1998) pointed out that their exospheric computation do not predict correctly the separate contributions to the temperature of arbitrary parts of the velocity distribution function, like, for instance, the “halo” population. Issautier et al. (2001) have compared the total SW electron temperature from second generation exospheric models and high-latitude Ulysses data, showing a reasonable agreement.

In summary, the density, bulk velocity, and pitch angle averaged temperature of electrons in the second generation of kinetic exospheric models are in agreement with observations. This is a consequence of replacing in the exospheric solution the Pannekoek-Rosseland polarization electric field by the self-consistent Lemaire-Scherer (LS) electric field. The electric potential is derived from the quasi-neutrality equation and the additional constraint that the thermal electrons do not evaporate out of the corona at a higher rate than the protons.

Some renewed interest for kinetic descriptions of the solar wind appeared when Scarf et al. (1967) and Hundhausen (1968) reported on the proton temperature anisotropy, with $T_{\perp}$ larger than $T_{\parallel}$ at 1 AU. They concluded that the interplanetary plasma should be more or less collisionless beyond a radial distance of 10 or 20 R _S . Indeed, none of the hydrodynamic SW models available at that epoch did account for the existence of temperature anisotropies for the protons and electrons, nor did they account for skewed and double-hump VDFs that can hardly be fitted by displaced Maxwellian velocity distribution functions.

Using still another exospheric approach, Griffel and Davis (1969) calculated that the temperature anisotropy resulting from a collision-free flow of protons, along spiral interplanetary magnetic field lines, is much larger than the corresponding values observed at 1 AU. Thus, these authors consolidated the latent argument that, between 0.1 and 1 AU, there should be some relaxation mechanism producing about 2.5 collisions/AU or pitch-angle scattering to reduce the excessive $T_{\parallel}/T_{\perp}$ in collisionless solar wind models. Eviatar and Schulz (1970) and Schulz and Eviatar (1970) reached similar conclusions based on their own kinetic description of solar wind protons and electrons.

Hundhausen’s (1970) observations show that the velocity distribution functions of electrons and protons are not isotropic. Similar results were obtained by Marsch et al. (1982) and Pilipp et al. (1987a) from Helios data. Indeed, a larger dispersion is generally observed for the velocities parallel to the magnetic field lines than for the perpendicular ones, i.e. the observed electron and proton temperatures are anisotropic. In the quiet solar wind, measurements consistently indicate that $T_{e\parallel}/ T_{e\perp} = 1.1 - 1.2$, $T_{p\parallel}/ T_{p\perp} = 2\pm 1$, on the average. Similar value are obtained in fast speed streams. Observations of polar coronal holes from the UVCS SOHO have evidenced temperature anisotropies in protons and heavy ions, but with $T_{p\parallel}/T_{p\perp} < 1$ at the base of the corona (Cranmer et al. 1999).

In all exospheric models, the temperature anisotropies are larger than observed. For instance in Lemaire-Scherer’s kinetic model LSb, at 1 AU: $T_{e\parallel}/ T_{e\perp} = 3.05$, and $T_{p\parallel}/ T_{p\perp} = 164$; in Jockers (1970) hybrid/semi-kinetic model at 1 AU: $T_{e\parallel}/ T_{e\perp} = 1$ and $T_{p\parallel}/ T_{p\perp} = 900.$ (see models {2}, {3} and {4} in Table 1). Chen et al. (1972) and Pierrard et al. (2001c) demonstrated that when a spiral interplanetary B-field is adopted, instead of being radial, the anisotropy of the proton temperature at 1 AU is significantly reduced, but not enough, however, to match the observed values.

The extreme values predicted for the exospheric anisotropies of electron and proton VDFs at 1 AU naturally suggest that some pitch angle scattering mechanisms operate in the solar ion-exosphere, as inferred in the early study of Griffel and Davis (1969). Lemaire and Scherer (1971) suggested for instance that Coulomb collisions, ignored in their exospheric formulation, might indeed change the pitch angles of the escaping particles, and will therefore attenuate the too large exospheric anisotropies of the electron and ion VDFs. According to Lemaire and Scherer (1971), the mean pitch-angle scattering times of solar wind protons and electrons are respectively 6 × 10⁵ s (7 days), and 4 × 10⁴ s (0.5 days) at 1 AU. These time scales are slightly smaller than the time required for the particles to travel a distance of 0.5 AU.

On the other hand, gyrotropic fluid models of the solar wind, which did include the effect of Coulomb collisions, predict that collisions were far from sufficient to prevent extremely large proton temperature anisotropies in the outer solar wind (Lie-Svendsen et al. 2002). The collision time for energy equipartition between electrons and protons is about 100 times larger than the time required for protons to travel a distance of 0.5 AU. As a consequence, at 1 AU no significant energy can be transferred from the hotter electrons to the solar wind protons nor to the other minor heavier ions present in the SW. Of course, wave-particle interactions may play a role in reducing the anisotropy of the VDF. This was already suggested by Lemaire and Scherer (1971), but to quantify their effect the actual frequency spectrum of these waves, their polarization, and radial distributions should be known with some confidence.

The effect of Coulomb collisions can be worked out by solving the Fokker–Planck equations for the different particle species in the coronal and solar wind plasma including for the minor ion species. Some preliminary attempts along these lines will be discussed later in this review. Before that, let us now examine some characteristic correlations observed in the solar wind and compare them to the predictions of exospheric and hydrodynamic models.

In the LSb model developed by Lemaire and Scherer (1971), the exobase altitude was taken to be at r ₀ = 6.5 R _S . This input parameter depends on the coronal electron density and of its scale height taken from eclipse observations up to 16 R _S (Pottasch, 1960). Since the mean free path of the coronal electrons and protons is proportional to the square of the temperatures of these particles, the altitude of the exobase (i.e. where the pitch angle scattering mean free path becomes equal to the atmospheric density scale height) in Lemaire-Scherer’s models is a decreasing function of the assumed exobase temperatures. Therefore, by decreasing the temperature of the electrons and protons at the exobase the corresponding exobase altitude and density can be varied in a self-consistent manner as described by Lemaire and Scherer (1971). As a result, a whole family of exospheric models has been generated by Lemaire and Scherer (1971) showing that an increase of the exobase temperatures in association with a decrease of the exobase altitude results in an increased proton temperature, T _p , and an enhanced bulk velocity, u _r , at 1 AU. However, a change of the exobase temperature did not significantly affect T _e , at the orbit of Earth (cf. Fig. 5 in Lemaire and Scherer 1971).

The positive correlation between the SW expansion velocity and the proton temperature found with Lemaire–Scherer’s exospheric models is illustrated in Fig. 4. It can be seen that over the range of velocities from 250 to 400 km/s, corresponding to quiet solar wind velocities,⁵ the positive correlation predicted by Lemaire-Scherer’s exospheric models (green dashed line) fits rather well the observations of Vela 3 and Explorer 34. Note that this agreement is well within the range of the standard deviations of observed SW parameters.

Fig. 4

Correlation between average proton temperature (T _p ) and solar wind bulk velocity (u _r ) at 1 AU for slow solar wind (panel a), fast solar wind (panel b). The green dashed line shows the positive correlation between the calculated values of u _r and T _p at 1 AU determined in an unpublished parametric study by Lemaire (1971), reported in Fig. 12 of Lemaire and Scherer (1973). The pink dots correspond to Vela 3 measurements reported by Hundhausen (1970). The orange solid line and the two orange dotted lines correspond respectively to the average relationship, the minimum and maximum standard deviations, as reported by Burlaga and Ogilvie (1970) based on their statistical study of Explorer 34 measurements. The star symbol labeled “1” corresponds to the average values of the quiet solar wind reported by Hundhausen (1972); it corresponds also to the mean values given in the column 1 of Table 1. The other numbered star symbols correspond to the ID numbers of various kinetic (E), hydrodynamic (H) and semi-kinetic/hybrid (Hb) models listed and referenced in Tables 1 and 2

Fig. 5

Exospheric distribution of the density (upper left panel), electric potential (upper right panel), bulk velocity (lower left panel), and total potential (lower right) for solar protons from the model by Lamy et al. (2003a) . The dotted curves refer to a second generation exospheric model based on a Maxwellian VDF at the exobase (r ₀ = 6 R _S ); the thin solid curves correspond to a third generation exospheric model for which the exobase is at the same altitude but for which the electron and proton VDF are Lorentzian functions whose index κ = 3. The thick solid line corresponds to an exospheric model whose VDF are also Lorentzian functions (κ = 3) but for which the exobase altitude is at 0.2 R _S , i.e. below the altitude where the total potential energy of the protons has its maximum value. In all the three models the coronal temperature of electrons and protons is the same, T ₀ = 1.5 × 10⁶ K

A positive correlation between u _r and T _p was also obtained by Hartle and Barnes (1970) with their set of two-fluid hydrodynamic models. To obtain a good agreement with SW observations at 1 AU extra heating rates of the protons had however to be added by Hartle and Barnes (1970) in a region extending up to 25 R _S . Indeed, in the hydrodynamic models, the increased proton temperatures and wind velocity at 1 AU are direct consequences of the enhanced heat deposition into the solar corona, while in Lemaire–Scherer’s exospheric models, it is a simple consequence of a parametric increase of the exobase temperatures.

Lopez and Freeman (1986) had interpreted the positive correlation between 〈T _p 〉 and u _r by a variation with velocity of the polytropic index in their hydrodynamic SW models. This tended to indicate that a larger energy input is plausibly the cause of both higher proton temperatures and bulk velocities at 1 AU. A similar theoretical interpretation for the observed (almost quadratic) relationship between 〈T _p 〉 and u _r has been proposed also by Matthaeus et al. (2006). They suggest that MHD turbulence should provide a heating depending only on the age of solar wind plasma elements (since their release from the coronal base). However, after a comprehensive quantitative analysis, Démoulin (2009) showed that MHD turbulence gives a much weaker dependence on u _r than the actual observations at 1 AU.

4.4.3 Limitations of Second Generation Exospheric Models

Like Euler’s hydrodynamic models, exospheric models should be viewed as zero-order approximations of more general kinetic models of the solar wind. However, so far the agreement shown in Figs. 3 and 4 is restricted to quiet/low-speed solar wind conditions; it does not apply to fast/high-speed streams where the bulk velocity exceeds 500–650 km/s. Such a high expansion velocity can possibly be attained for exobase electron temperature higher than 3 × 10⁶ K; unfortunately, this is much higher than usual coronal temperatures inferred from spectroscopic or white light eclipse observations.

Note that a similar restriction is applicable to hydrodynamic models of the solar wind; indeed non-realistic high coronal temperatures (larger than 3 × 10⁶ K) would be required in hydrodynamic models to attain solar wind velocities larger than 500 km/s at 1 AU. Otherwise, conjectured heat deposition or/and ad hoc momentum boosting of the coronal plasma are needed up to 25 R _S or beyond, in order to reach supersonic velocities of 700 km/s or more. Unfortunately, this contradicts the observations collected in fast speed streams. Indeed, it was discovered (see, e.g. Schwenn et al. 1978) that fast speed streams are consistently originating from coronal holes, i.e., from regions of the solar corona where the electron temperature is smaller than in the equatorial regions, i.e. the source of the quiet solar wind.

Therefore, neither the standard collision-dominated hydrodynamic models nor the collisionless exospheric models of the solar wind developed so far can represent adequately the phenomenon of fast speed solar wind streams; their applicability must be limited to the quiet or slow solar wind flow originating from the hotter and denser equatorial regions of the corona.

4.5 The Third Generation of Kinetic/Exospheric Models

In kinetic exospheric models of the SW the altitude of the exobase is energy dependent. This has been already pointed out by Pikel’ner (1950). In the exospheric models discussed in the previous paragraphs, the mean energy of particles was considered to define a sharp cut exobase level: a step like transition region between the collision dominating plasma and the collisionless exosphere. This procedure is a heritage of Jeans’ exospheric theory; to the best of our knowledge there are no successful attempt to develop exospheric SW models with energy dependent exobase levels.

The early kinetic exospheric models were based on the assumptions that: (1) the Knudsen number changes from a small value to a value much larger than unity over a relatively thin layer, (2) the velocity distribution functions at the exobase could be approximated by truncated Maxwellians, and (3) the exobase level is at altitudes high enough in the corona, so that the ratio between the electric force and gravitational force acting on protons |e E/m _p g| is larger than unity in the whole solar ion-exosphere, i.e., that the repulsive electric force accelerating the protons in the exosphere is everywhere larger than the attractive gravitational force slowing them down (cf. Fig. 2).

The kinetic modeling of the solar wind received new credit in the second half of the 1990s. Some of the early, simplifying assumptions used in exospheric models were questioned and less restrictive constraints have been adopted to reach results in better agreement with the challenging fast speed streams observed in the solar wind.

A first improvement was proposed by Pierrard and Lemaire (1996), who developed an exospheric theory based on a Lorentzian VDF (or “kappa” function), instead of the usual Maxwellian VDF. Isotropic Lorentzian VDFs, f _κ(v), are characterized by three independent parameters: a “kappa” index, κ, a density, n ₀, and a most probable thermal speed, v _th , defined by :

$$ v_{th}^2 = \frac{\left(2 \kappa -3\right) k T_0} {\kappa m} $$

(42)

where, T ₀, is the exobase temperature. The Lorentzian velocity distribution function is written:

$$ f_\kappa(v, r, t) = n_0 \frac{\Upgamma(\kappa+1)}{\Upgamma (\kappa - 0.5 ) (\pi \kappa {v_{th}}^2)^{3/2}} \left( \frac{1+ v^2}{\kappa {v_{th}}^2} \right)^{-(\kappa+1)}. $$

(43)

where $\Upgamma$ is the Gamma function. Similar expressions with different values for κ _e and κ _p can be employed for the electron and proton VDFs at the exobase.

Lorentzian (kappa) VDFs decrease as a power law of v, instead of varying exponentially like a Maxwellian; a relatively large fraction of particles populate the tail of the VDF (i.e. with velocities larger than the critical escape velocity). In the limit κ → ∞ the Lorentzian function approaches a Maxwellian distribution function with the same values of n ₀ and T ₀. A review of kappa distribution functions and their applications in space plasma physics has been recently published by Pierrard and Lazar (2010).

The existence of Lorentzian VDF characterized by given values of κ, is treated in third generation exospheric models without discussing the fundamental physical mechanism that would be able to account for the non-Maxwellian power law tails. That fundamental yet unresolved question has been addressed in a few recent theoretical papers. A comprehensive study reporting earlier attempts can be found in Shizgal (2007). It is shown that, in a uniform plasma, the relaxation toward equilibrium does not necessarily tend to kappa distributions, even when a certain type of wave-particle interaction is taken into account. But this does not exclude that such suprathermal tails could not be produced by other processes, including the presence of a gravitational force in association with vertical density gradients. Another recent study by Livadiotis and McComas (2009) is addressing this fundamental question: what could energize the particles inside the solar corona, and produce a population of suprathermal electrons and protons at the exobase, or even closer to the Sun? These authors show how such non-Maxwellian distributions may arise naturally from Tsallis statistical mechanics. Another explanation for the formation of suprathermal electrons, in particular the magnetic-field-aligned component (or “strahl”) was proposed by Vocks et al. (2005) in terms of resonant interaction with anti-sunward propagating whistler waves. A recent review on wave particle interactions and whistler waves can be found in the book by Trakhtengerts and Rycroft (2008).

As a consequence of a postulated enhanced population of suprathermal electrons with energies above the critical energy, the Jeans evaporation flux of electrons will be larger for a Lorentzian VDF than for a Maxwellian VDF given by (36). Since the evaporation flux of coronal electrons emitted by the corona is then larger, a larger electrostatic potential difference, $\Updelta\Upphi_E$, is needed to cancel the net electric current out of the corona. As a consequence of the larger electric potential difference between the exobase and infinity, the protons are accelerated to higher supersonic velocities.

This is illustrated in the four panels of Fig. 5 taken from Lamy et al. (2003b), where the curves show the exospheric distributions of the density, n(r), electric potential, $\Upphi_E(r)$, plasma bulk velocity, u _r (r), and total potential (the sum of the gravitational and electric potential) of protons, R _p (r), in three different exospheric models. It can be seen that the bulk velocity in the exospheric model based on kappa VDF provides larger supersonic velocities (>450 km/s) than in the case of a Maxwellian VDF for which κ→∞. Note that changing the value of κ does not change significantly the plasma density in the ion-exosphere: the latter remains in excellent agreement with the average density observed both in the quiet or slow solar wind and in fast/high-speed solar wind regimes (Hundhausen 1970; Ebert et al. 2009).

An interesting parametric study has been published by Zouganelis et al. (2004) where non-thermal VDFs have been used at the exobase of the solar corona, instead of truncated Maxwellians. Their study illustrates the key role played by the population of supra-thermal electrons in the process of acceleration of the solar wind ions to large terminal speeds. Zouganelis et al. (2004) showed that the high terminal SW speeds obtained in their revisited exospheric models using Lorentzian (kappa) VDF at the exobase is not just an artifact of the type of kappa function, but is a general feature when the tail of the electron VDF is parametrically enhanced. Indeed, similar high terminal speeds are obtained by using a sum of two Maxwellian VDFs, i.e. a first one with normal coronal temperature, representing the core electrons, and a second one with a higher temperature corresponding to the halo electrons observed in the SW. They found that when the kappa index is reduced, or when the relative abundance of the halo electrons is enhanced with respect to the core electrons, the SW bulk speed at 1 AU is systematically enhanced.

By reducing the index κ from large values to less than 3, or by increasing T ₀, the exobase temperature, above 2 × 10⁶ K, it is possible to increase the solar wind velocity at 1 AU to even higher values. However, it seems unrealistic to push these two parameters to such extreme values in order to obtain bulk velocities as high as 700–800 km/s reported in fast speed solar wind streams. Since fast speed streams originate from coronal holes where the electron temperature is less than 10⁶ K it appears that some additional modification is required to boost exospheric bulk velocities up to 1,000 km/s. This is what is tentatively proposed in the studies reported in the next section.

The mean free path of the particles is larger in coronal holes than in the equatorial corona at the same altitude and for the same temperature. As a matter of consequence the exobase is located at a lower altitude in coronal holes than in the equatorial region of the corona. This implies that the gravitational force is larger than the electric force acting on protons. This is evidenced in Fig. 2. As a consequence of the lowering of the exobase altitude the total potential energy of protons, R _p (r), remains dominated by the gravitational potential energy at the base of the solar ion-exosphere. Therefore, just above the exobase, R _p (r) is increasing with r. It reaches maximum value R _p,max at r = r _max where |e E/m _p g| = 1; beyond this heliocentric distance, R _p (r) decreases with r, just as in the earlier exospheric model where the exobase was located higher in the equatorial corona. If the exobase altitude is significantly lowered, R _p (r) is no more a monotonically decreasing function of r, as assumed in the second generation exospheric models of Lemaire–Scherer. The non-monotonicity of the total potential energy of the protons, anticipated by the early kinetic exospheric model of Jockers (1970), is illustrated by the thick solid line in the lower right hand side panel of Fig. 5. For this coronal hole exospheric model, r ₀ = 1.2 R _S and r _max  = 1.9 R _S . Below r _max , protons are attracted toward the Sun by the dominating gravitational force; ballistic or captive protons are then populating this region in addition to the escaping ones. The former class of particles do not have a kinetic energy large enough to reach over the maximum potential barrier, R _p,max.

In the third generation of kinetic SW models, new analytical expressions were found for the densities, escape fluxes, expansion velocity, pressure tensors, temperatures, energy and heat fluxes by Pierrard and Lemaire (1996), Lamy et al. (2003a). These authors performed an exact exospheric treatment of a non-monotonic total (electric and gravitational) potential for the protons, also treated by the kinetic exospheric model of Jockers (1970).

The analytical expressions derived by Pierrard and Lemaire (1996), Lamy et al. (2003a) depend, of course, on the value of r _max and of R _p,max. These parameters must first be calculated by an iterative procedure which has been developed and described by Lamy et al. (2003a). This is how the exospheric distributions illustrated by the thick solid curves in Fig. 5 were obtained.⁶

The thick solid line in the third panel of Fig. 5 shows the radial distribution of u _r obtained in this new exospheric model for which the exobase is here at r ₀ = 1.2 R _S , the exobase temperature is 1.5 × 10⁶ K and κ = 3 for both the VDFs of the electrons and the protons. It can be seen that the supersonic expansion velocity reaches now more than 600 km/s, which is larger than the values obtained in earlier exospheric models with exobase altitudes much higher up in the equatorial region of the corona. Nevertheless, values as large as 1,000 km/s can hardly be obtained even with this ultimate improvement of the exospheric models for the solar wind expansion. The plasma density in this new kind of exospheric model for coronal holes is shown in the upper left panel of Fig. 5. As a consequence of the lower exobase density adopted to match coronal holes densities, the solar wind density at 1 AU is also lower in this fast speed stream exospheric model, which is consistent with the observations. One limitation of the third generation exospheric models is a relatively high maximum temperature reached by SW electrons at some distance from the Sun, for highly nonthermal electron VDFs at the exobase.

This problem could be due to the approximations used by the model to compute the self-consistent electric potential or a possible superradial expansion of the wind (Zouganelis et al. 2004). Nevertheless, the third generation of exospheric models by Lamy et al. (2003a) and Zouganelis et al. (2004) give a description of the transonic region, a major advance that enabled also to obtain large SW terminal velocities at 1 AU. Further improvement should provide a more accurate representation of the solar wind, including of the electron temperature.

As a consequence of the wide variety of adjustable parameters, the modern exospheric models behold a great flexibility and allow to explain not only qualitatively, but quantitatively as well, a broad range of SW features observed at 1 AU; despite the rather artificial truncation of the exospheric VDFs, exospheric models offer the substantial advantage of giving clues to understand the physics of the acceleration of the solar wind flow to supersonic terminal speeds. Such an advantage is limited, however, in single-fluid hydrodynamic SW models, where the electrostatic force accelerating the protons to such speeds is actually hidden behind the gradient of the kinetic electron pressure (Parker 2010).

Note that in any exospheric models of the solar wind, no adjusted/ad hoc heat or momentum deposition rate has been postulated, so far. But, based on the results reviewed above, additional heating mechanisms and acceleration of coronal protons seem indeed to be in order to drive fast speed streams up to bulk velocities of 1,000 km/s. Let this then be a provisional conclusion of the present review after half a century efforts to model the solar wind flow in the framework of exospheric theory which is inherited from J.H. Jeans.

4.6 The Fourth Generation of SW Kinetic Models: Fokker–Planck Solutions

Like any physics-based model, the exospheric approach has its limitations: (1) it is restricted to stationary models; (2) the subtle effects of rare collisions and of wave-particle interactions above the exobase are ignored and (3) it cannot account for processes near and below the exobase.

Fortunately, the effect of infrequent Coulomb collisions can be dealt with by solving the Boltzmann equation with the Fokker–Planck (FP) collision term (see Sect. 2.2), for electrons or the protons, or for both species together. The FP collision term is, indeed, more appropriate than the Boltzmann collision integral (5). In plasmas the cumulative effects of long range Coulomb collisions is more efficient than close binary collisions to produce momentum transfer, pitch angle deflection, and thermalization of the different species of charged particles.

Lie-Svendsen et al. (1997) solved the Boltzmann equation with the Fokker–Planck collision term (8) to develop the first collisional model of the solar wind. They determined the VDF of the test electrons in the exobase transition region. They used a standard finite difference numerical method to calculate the VDF at preset altitudes. A similar Fokker–Planck equation for the solar wind electrons was solved by Pierrard et al. (1999). They developed a spectral method to determine the polynomial expansion of the electron VDF across the exobase transition and adapted a numerical code developed earlier to solve the FP equation for the polar wind (Pierrard 1997). The model by Pierrard et al. (1999) used correct boundary conditions at the exobase, satisfying the regularity conditions determined by Pierrard (1997). In a more recent application Pierrard et al. (2001b) added collisions with protons. Pierrard et al. (1999) integrated the FP equation for the solar wind electrons from 1 AU to 4 R _S in the solar corona. A typical VDF of core-halo electrons, observed at 1 AU in high speed solar wind stream by the WIND spacecraft, was used as boundary conditions at the top of the integration domain in their FP model calculation. Figure 6 illustrates the VDF of electrons given by this model at two different radial distances in the solar wind (r = 215 R _S and and 4 R _S ). The asymmetry of the electron VDF observed at 1 AU, which can be interpreted as forming the strahl population of high speed SW electrons, maps down into the solar corona as a slight asymmetry of the coronal VDF.

Fig. 6

The contour plots show the iso-density contours of $f_e(y_{perpendicular}, y_{parallel})$, at two different radial distances: r = 215 R _S (left panels), r = 4 R _S (right panels); $y_{perpendicular}$, y _parallel denote the perpendicular respectively the parallel component of the velocity. The dashed circles correspond to the thermal speed of the electrons. The line-plots show the cross sections of $f_e(y_{perpendicular}, y_{parallel})$ in the radial direction, which is assumed here to be the direction of the interplanetary magnetic field; the dashed lines in these two panels show the corresponding cross-sections of the electron VDF in the direction perpendicular to the magnetic field; figure from Pierrard et al. (1999)

From their study based on the Fokker–Planck equation Pierrard et al. (1999) infer that, in order to match the halo and strahl electron VDFs observed at 1 AU, a suprathermal tail should be present in the VDF of the electrons deep into the solar corona. In other words, the features characterizing the suprathermal electrons observed at 1 AU cannot be generated within the interplanetary medium solely by Coulomb collisions on the way between the corona and 1 AU; according to their results an origin deep in the corona has to be postulated. Alternatively, kinetic instabilities or wave-particle interactions may be responsible for producing specific features in high speed SW streams.

Models of solar wind VDFs based on coupled Fokker–Planck equations for electrons and protons remain to be developed in the future. Using coupled and time dependent Fokker–Planck equations for solar wind electrons and protons is a widely open field of investigation: an unexploited avenue for scientists interested in advancing basic plasma physics for applications to space weather issues.

5 Fluid Models of the Solar Wind

Instead of considering direct solutions to the fundamental Boltzmann equation (with varying degrees of simplification), to calculate the velocity distribution f(r, v, t) and its moments, one can consider models that “only” provide averages of the VDF, such as particle density, mean flow velocity and temperature (or mean kinetic energy). This is obviously a much more “coarse grained” picture than provided by the full VDF. One may therefore be misled into thinking that fluid models are inferior to kinetic models. This is in general not the case, however. Although the Boltzmann equation (2) looks deceptively simple, obtaining actual numerical solutions for a realistic physical system can be a formidable task. Even in the simplest, spherically symmetric solar wind outflow, at least two velocity components must be included, and obtaining a steady state solution entails solving a partial differential equation in three dimensions (one space, two velocity space). With the full collision term (5) or (8), included, it even becomes a nonlinear integro-differential equation. As we have seen in the previous section, it can only be solved numerically (or analytically in very simple situations), such as steady state, one-dimensional, collisionless flow.

5.1 Advantages and Limitations of Fluid Models

Fluid models, by providing a more coarse grained picture of the plasma, can afford to do much more:

• Collisions can be included (although the terms can become very complicated and approximations have to be made, particularly for higher-order fluid models).

• Processes such as ionization and recombination may easily be included.

• Many particles species, and their interactions, may be included simultaneously.

• Interaction of matter with radiation may be accounted for.

• The fluid equations may be coupled with Maxwell’s equations, creating e.g. magnetohydrodynamic (MHD) equations, and thus accounting for the interaction between the plasma and the electric and magnetic fields.

• A very large range of spatial scales may be included simultaneously, ranging from the solar upper atmosphere and to the heliopause, if need be.

• One is not restricted to stationary one-dimensional flow, and with today’s computational resources obtaining time-dependent solutions for three-dimensional outflow, e.g., to model coronal mass ejections, is feasible.

Because fluid models allow us to include more processes and spatial scales, they are well suited as “global” models of the solar wind, while kinetic models are more suited to specific processes that cannot easily or reliably be described by fluid models. For instance, the actual evolution of the VDF in a collisionless plasma requires solving the Boltzmann (or, more correctly, Vlasov) equation. As we shall see, fluid models have their limitations and uncertainties, which can only be overcome in a kinetic description.

The list above also illustrates that the particle transport equation(s) (either Boltzmann or fluid equations) are but one ingredient in a model describing the solar wind system. Including these other processes requires that the transport equations themselves must be kept as simple as possible. That also applies to the fluid models. Depending on the assumptions made when deriving them, they may quickly attain a complexity rivalling that of the Boltzmann equation, particularly so-called higher-order fluid models. What is potentially gained by their more detailed description of the plasma may quickly be lost in the difficulties obtaining numerical solutions to them, or because we no longer can afford to include some of the other processes listed above. Fluid equations have an additional advantage over kinetic models: We have many decades of experience developing algorithms for solving them numerically. In fact, the vast field of computational fluid dynamics (CFD) is devoted to this.

Fluid equations describe the evolution of the moments of the VDF, such as the density, (mean) flow velocity, and temperature. As discussed in Sect. 2.4, these equations may be formally derived by multiplying (2) by some velocity moment v _i v _j … and integrate over v. However, because of the v · ∇ _r f term in (2), the equation for a moment of order N (in velocities) will contain a moment of order N + 1. Hence this procedure leads to an infinite series of coupled partial differential equations, which has to be closed by expressing a higher-order term (moment) as some analytic function of lower-order terms. In a collision-dominated gas the closure assumption need not be problematic; for instance, one may safely assume that the heat flux (a higher-order moment) is proportional to the temperature gradient (a lower-order moment). However, in collisionless gases, where the departure from a Maxwellian VDF may be large, it is not obvious how the set of equations should be closed.

A second problem with fluid equations is how to evaluate the collision term. Because δf/δt _coll is an integral over the product of the VDFs of two colliding species, it will necessarily involve higher order moments. Moreover, except for the case of Maxwell molecule interactions (in which case the collisional cross section is inversely proportional to the relative velocity of the colliding particles, $\sigma \propto 1/|{\bf v} - {\bf v}'|$), the collision term cannot be evaluated without making explicit assumptions about the shape of the VDFs (Burgers 1969). Since Maxwell molecule interactions are not relevant for a fully ionized plasma (but rather for collisions between ions and neutrals), this means that solar wind fluid equations cannot be developed without assuming an explicit, analytic form for the VDF. This analytic form will also provide the closure for the moment equations. Except for the simplest possible approximation, the so-called 5-moment approximation, which is not adequate for the corona and solar wind as we shall see, the collision terms can become unwieldy and further approximations have to be made in order to obtain analytical expressions. Typically one has to assume that the flow speed difference between species is much smaller than thermal speeds, $|{\bf u}_s .-{\bf u}_t| \ll\sqrt{kT_{s(t)}/m_{s(t)}}$, which limits the applicability to the transition region and corona where collisions are sufficiently strong to prevent large flow speed differences. Finally, particularly for Coulomb collisions the collision terms will be very sensitive to the shape of the VDF. Without a careful choice for the VDF, large errors may result, e.g. in the forces acting between particles and in the heat flux coefficient, even in the lower corona where particles are subsonic.

Although fluid models have their advantages over kinetic models, and for computational reasons we often have no choice but using a fluid description, the above discussion shows that fluid models also have their shortcomings. They depend strongly on the assumptions made on the VDF and the closure relationship. One must carefully select the “right” model for the problem one wants to investigate.

The presentation below must necessarily be short, and will barely do justice to the topic. For the reader searching for a thorough presentation of the topic, the classic textbook of Chapman and Cowling (1970) provides a detailed account of diffusion, thermal conduction, and viscosity in non-uniform gases. The monograph by Burgers (1969) provides a thorough description of collisional effects, and presents the techniques necessary to evaluate the collision integrals analytically. More recently, Gombosi (1994) provides an introduction to the kinetic theory of gases, including so-called generalized transport equations. In addition to the extensive treatment of Grad (1958), discussed in the Appendix 2, Schunk (1977) offers a more “practical” and concise presentation of transport equations for use in aeronomy, at various levels of complexity, while Barakat and Schunk (1982) review transport equations for anisotropic space plasmas.

5.2 Properties of the Solar Wind and Upper Solar Atmosphere Relevant for Fluid Modeling

Many of the basic properties of the solar wind are set in the subsonic region of the flow. Models suggest that the elemental abundances of the solar wind may even be set in the upper chromosphere, where the gas becomes ionized (e.g. Hénoux 1998). Observations indicate that the wind is rapidly accelerated within a few solar radii from the Sun (Habbal et al. 1995b; Grall et al. 1996; Antonucci et al. 2000), implying that most of the energy flux carried by the solar wind has been deposited in the corona.

Moreover, the solar wind mass flux is not only set in the subsonic region of the flow, but will largely be determined by the energy balance in the solar transition region (the interface between the chromosphere and corona). The mass flux depends on the coronal pressure, which is approximately equal to the transition region pressure. In the absence of flow the transition region pressure is in turn directly proportional to the downward heat flux density from the corona (Landini and Monsignori Fossi 1975). Hence the coronal density, and thus the solar wind mass flux, depend sensitively on the energy balance of the transition region (Lie-Svendsen et al. 2001).

This implies that if obtaining, e.g., the elemental composition or the mass flux of the solar wind is an important goal of the model, it must extend all the way down to the upper chromosphere. In such SW models it is of course essential that Coulomb collisions are described accurately, as e.g., heat is transported in a plasma dominated by collisions. Since heat conduction is essential in the transition region and corona, one must also choose a fluid description that allows for heat conduction (but, as we shall see later, the simplest possible fluid description does not).

The Coulomb interaction differs from the other collision types (neutral-neutral and neutral-ion collisions) by being strongly energy dependent. In fact, the cross section for Coulomb collisions is proportional to the velocity difference of the colliding particles to the inverse fourth power (6). As a consequence, the collision frequency, assuming Maxwellian VDFs, is strongly temperature dependent, ${\nu}\,{\propto}\,T^{-3/2}$. The temperature increases, particularly for protons and heavier ions, from the transition region into the corona, while the density decreases with altitude. The combined effect is a rapidly decreasing collision frequency, resulting in a steep transition from collision-dominated to collisionless flow for protons and ions, somewhere in the corona.

Another effect of the velocity dependence of the Coulomb cross section is the strong temperature dependence of heat conduction in a fully ionized, collision-dominated gas, where the electron heat flux is (Spitzer and Härm 1953; Braginskii 1965)

$$ {\bf q}_e = - \frac{\kappa_e'}{\ln \Uplambda} T^{5/2} \nabla T, $$

(44)

with $\kappa_e' = 1.84 \times 10^{-10}\,{{\hbox{W}}\,{\hbox{m}}}^{-1}\,{{\hbox{K}}}^{-7/2}$ and $\Uplambda$ the plasma parameter specified in Sect. 2.2 accounting for the Debye shielding of the electric field. The T ^5/2-dependence of q _e means that heat conduction is virtually absent in a cold plasma, while it becomes “suddenly” very efficient at coronal temperatures of the order of 10⁶ K. Heat conduction thus acts like a thermostat for the electrons: little heating is required to attain temperatures of order 10⁶ K, while it is difficult to attain temperatures much higher than this (except, perhaps, under extreme conditions such as in a solar flare). This is also in agreement with, e.g., observations of coronal holes, the source regions of the fast solar wind, indicating electron temperatures slightly below 1 × 10⁶ K (Wilhelm et al. 1998). The efficient heat conduction at these temperatures also implies that electrons will not become collisionless in the corona, and as the temperature decreases outwards in the solar wind they may even remain collisional far into the solar wind, despite the decreasing density. Note also from (44) that heat conduction is almost independent of the coronal density (with only a very weak density dependence through the $\ln\Uplambda$ term).

Collisions quickly cease to be important for protons and ions at larger distance from the Sun, where the solar wind is instead influenced by external forces (gravity and electric and magnetic forces), pressure forces within the plasma, and wave-particle interactions. Here a fluid model needs to ensure that the effects of collisions are sufficiently small to be unimportant or small first order corrections. However, the lack of collisions means that the VDF may, and will, become highly anisotropic (and possibly non-Maxwellian). In fact, observations of coronal holes from the SOHO satellite show that protons and oxygen ions are highly anisotropic even in the corona, with the temperature along the line of sight ($T_\perp$) higher than the temperature perpendicular to the line of sight ($T_\parallel$) (Cranmer et al. 1999). In situ observations by the Helios spacecraft show that protons are highly non-Maxwellian and anisotropic with $T_\perp > T_\parallel$ also at 0.3 AU in the fast solar wind (Marsch et al. 1982). $T_\perp$ and $T_\parallel$ refer to directions relative to the SW magnetic field. Since adiabatic expansion of the wind should have the opposite effect at large distances, leading to $T_\perp \ll T_\parallel$, because of magnetic moment conservation, this observation shows that some form of wave-particle interaction could have affected the proton VDF. Electrons on the other hand remain close to Maxwellian even at such large distances from the Sun, with the core of the VDF being quite isotropic but with an antisunward high-velocity “strahl” in high-speed solar wind streams (Pilipp et al. 1987a). To capture such proton anisotropies in a fluid description, one must consider so-called gyrotropic equations, in which the VDF is expanded about a bi-Maxwellian with different temperatures parallel and perpendicular to the magnetic field.

Below we shall separate the fluid descriptions into models in which the VDF is expanded about a Maxwellian, and models in which it is not (including expansion about a bi-Maxwellian). Although there is definitely overlap between the two classes of models, one may regard the former as more suited to describe the transition region, corona, and inner solar wind, while the latter group is mostly suited to the collisionless, supersonic solar wind where VDFs deviate strongly from isotropic Maxwellians. We shall only include elastic Coulomb collisions, although models extending into the low corona and transition region should include ionisation, recombination, and charge exchange reactions, as well as neutral-neutral and neutral-ion collisions. This is certainly the case when minor elements are taken into account in the model. Since the VDF should be close to a Maxwellian where and when these processes are important, the corresponding collision terms can be derived fairly straightforwardly if required.

5.3 Models for the Corona and Inner Solar Wind

When collisions, which restore a Maxwellian VDF in a uniform gas, are important, it is generally considered as reasonable to assume a VDF for species s of the form

$$ f_s({\bf r}, {\bf v}, t) = f_{s}^{(0)}({\bf r}, v) \left( 1 + \phi_s({\bf r}, {\bf v}, t) \right) $$

(45)

where $f_{s}^{(0)}$ is a drifting Maxwellian,

$$ f_{s}^{(0)} \equiv n_s \left( \frac{m_s}{2\pi k T_s} \right)^{3/2} \exp\left( - \frac{m_s c_s^2}{2kT_s} \right). $$

(46)

Here m _s is the mass of a particle and k is the Boltzmann’s constant, while c _s ≡v − u _s (r, t) is the “peculiar” velocity as measured in the reference frame moving at the species drift velocity u _s .

Irrespective of the form of ϕ _s , we shall always insist that n _s , u _s and T _s are the density, flow velocity and temperature of species s (moment definitions are summarized in Appendix 1):

$$ n_s = \int \int \int f_s d{\bf v} $$

(47)

$$ {\bf u}_s =\frac{1}{n_s} \int \int \int {\bf v} f_s d{\bf v} $$

(48)

$$ T_s = \frac{m_s}{3kn_s} \int \int \int ({\bf v} - {\bf u}_s)^2 f_s d{\bf v} . $$

(49)

The correction term ϕ _s accounts for the departure from a Maxwellian, and will generally be some polynomial of order N in v [or more correctly, in c ≡v − u _s (r, t)]. However, because it is postulated that u _s should be the mean flow velocity of the gas, higher order terms will have to be included to ensure this. Moreover, to be able to evaluate the collision terms analytically one generally has to assume that |ϕ _s | ≪ 1, so that terms proportional to ϕ _s ϕ _t may be neglected in the collision integrals. This requires that flow velocity differences are small in comparison with thermal velocities.

5.3.1 The 5-Moment Approximation, ϕ _s  = 0

The simplest choice for the VDF that one assumes first is ϕ _s  = 0. Multiplying (2) by 1, v, and v ², and using (47)–(49), one obtains the moment equations (20) and [using the notation of Schunk (1977)],

$$ m_s n_s \frac{D_s {\bf u}_s}{Dt} + \nabla p_s - m_s n_s {\bf g} - e_s n_s ( {\bf E} + {\bf u}_s \times {\bf B}) = \frac{\delta {\bf M}_s}{\delta t} $$

(50)

$$ \frac{D_s}{Dt} \left( \frac{3}{2} p_s \right) + \frac{5}{2} p_s \nabla \cdot {\bf u}_s =\frac{\delta E_s}{\delta t}. $$

(51)

where p _s  = n _s kT _s is the thermal pressure. When summed over the number of component species, (50)–(51) give the Euler equations (27)–(28), see definitions of moments in Appendix 1. When ϕ _s  = 0 the collision terms can be obtained analytically without making further approximations, and the terms on the right hand side of (50)–(51) may be written. These expressions, given below, are valid for arbitrary temperatures and flow velocities:

$$ \frac{\delta{\bf M}_s}{\delta t} = m_s n_s \sum_t \nu_{st} ({\bf u}_t - {\bf u}_s) \Upphi_{st} $$

(52)

$$ \frac{\delta E_s}{\delta t} = \sum_t \frac{m_s}{m_s+m_t} n_s \nu_{st} \left[ 3 k (T_t-T_s) \Uppsi_{st} + m_t ({\bf u}_s - {\bf u}_t)^2 \Upphi_{st} \right], $$

(53)

where the sums extend over the species with which the s-specie collides; ν _st is the collision frequency between species s and t, and $\Upphi_{st}$ and $\Uppsi_{st}$ are correction factors.

For Coulomb collisions the collision frequency reads (in SI units):

$$ \nu_{st} = \frac{1}{3} n_t \frac{m_t}{m_s+m_t} \left( \frac{2\pi k T_{st}}{\mu_{st}} \right)^{-3/2} \frac{e_s^2 e_t^2}{\varepsilon_0^2\mu_{st}^2} \ln\Uplambda $$

(54)

The reduced mass and temperature are defined by:

$$ \mu_{st} \equiv \frac{m_s m_t}{m_s+m_t} $$

(55)

$$ T_{st} \equiv \frac{m_tT_s+m_sT_t}{m_s+m_t} $$

(56)

For Coulomb collisions the velocity dependent correction factors read

$$ \Upphi_{st} =\frac{3}{2\varepsilon_{st}^2} \left( \frac{\sqrt{\pi}} {2} \frac{{erf}(\varepsilon_{st})}{\varepsilon_{st}} - \exp(-\varepsilon_{st}^2) \right) $$

(57)

$$ \Uppsi_{st} = \exp(-\varepsilon_{st}^2), $$

(58)

where erf(x) is the error function and

$$ \varepsilon_{st} = \frac{|{\bf u}_s - {\bf u}_t|}{\sqrt{2kT_{st}/\mu_{st}}}. $$

(59)

To a good approximation the solar wind is a fully ionized electron-proton plasma. Assuming quasi neutrality and no current, and neglecting all terms proportional to the electron mass m _e , the electron and proton momentum equations (50) may be added to read simply

$$ m_p n \frac{D{\bf u}}{Dt} + \nabla(nk(T_e+T_p)) - m_p n {\bf g} = 0, $$

(60)

where now n _e  = n _p  = n, u _e  = u _p  = u and the right-hand side is zero from Newton’s 3rd law (there may still be collisions).

The most serious shortcomings of the 5-moment approximation is the lack of heat conduction. As emphasized in Sect. 5.2, heat conduction, particularly for the electrons, is essential in the corona, as the high temperature there makes it extremely efficient. Mainly for that reason, this description is not acceptable for the corona and solar wind. It has other shortcomings as well; for instance, the frictional force between charged particles is too large (for a given u _s  − u _t ), implying a too large electric resistivity. For an electron-proton solar wind plasma these terms vanish in any case, as assumed in (60), in cases where electric currents or other particles (such as helium) are included, the errors will be large.

The 5-moment approximation is attractive because it is simple, and the collision terms are exact and valid for all flow conditions. When one allows for ϕ _s  ≠ 0, to account for heat conduction the collision terms will only be valid for flow speeds much smaller than the thermal speed, ε _st ≪ 1. When the solar wind becomes supersonic this condition may not be met.

In solar wind models with α-particles (or minor ions) included, spurious coupling between α-particles and protons may occur in the supersonic solar wind. Because the temperature tends to decrease as the wind is accelerated (even with heat conduction), and $\nu_{st} \propto T^{-3/2}$ from (54), the frictional coupling between protons and α-particles may become large if the flow speed difference is large. However, the Coulomb cross section (6) shows that the coupling should be small if the velocity difference is large, which is handled by the $\Upphi_{st}$ term in the 5-moment description. To overcome this problem, the collision terms are “repaired” by introducing ad hoc velocity corrections. An obvious solution is to apply the 5-moment velocity corrections (57) and (58), also to higher-order fluid models where they are not formally correct. Similarly, the second term on the right-hand side of (53) represents frictional heating, which can also be important when, e.g., α-particles stream faster than protons in the outer corona. Again, this effect will be absent in the following higher-order approximations with ϕ _s  ≠ 0 because the computation of the collision terms is difficult when flow speed differences become large. This illustrates that a procedure that is formally (mathematically) correct is not necessarily physically correct, and that the seemingly more accurate higher-order approximations that we shall discuss next may be, in reality, no more accurate than even the “crude” 5-moment approximation.

5.3.2 The 8-Moment Approximation

To account for heat conduction, one must have a skewness in the VDF, so that ϕ _s  ≠ 0. The simplest choice is a function linear in c _s , $\phi_s \propto {\bf a} \cdot {\bf c}_s$. Recall that c _s is the velocity of particles s measured in the frame moving with the average velocity, c _s  = v − u _s . Since ϕ _s must be a scalar function, one takes the scalar product with some vector a. This may be regarded as the first order Taylor expansion, in powers of c _si , of the “true” ϕ _s . However, in order to enforce that the average thermal velocity is equal to zero, i.e. $\int {\bf c}_s f_{s0}(1+\phi_s) d{\bf c}_s = {\bf 0}$, so that u _s is the mean flow velocity, one must add a third order term as well. The 8-moment assumption for ϕ _s may then be written

$$ \phi_s = - \frac{m_s}{kT_sp_s} \left( 1 - \frac{m_s c_s^2} {5kT_s} \right) {\bf q}_s \cdot {\bf c}_s, $$

(61)

where the coefficient in front is chosen such that q _s is the heat flux density,

$$ {\bf q}_s \equiv \frac{1}{2} m_s \int c_s^2 {\bf c}_s f_s d{\bf c}_s. $$

(62)

With this choice, the left-hand side of the momentum equation (50) remains unchanged, while the heat flux divergence is added as a source or sink of thermal energy to (51):

$$ \frac{D_s}{Dt} \left( \frac{3}{2} p_s \right) + \frac{5}{2} p_s \nabla \cdot {\bf u}_s + \nabla \cdot {\bf q}_s = \frac{\delta E_s}{\delta t}. $$

(63)

Taking the appropriate moment of (2), the equation for q _s becomes

$$ \frac{D_s {\bf q}_s}{Dt} + \frac{7}{5} ({\bf q}_s \cdot \nabla) {\bf u}_s + \frac{7}{5} {\bf q}_s (\nabla \cdot {\bf u}_s) + \frac{2}{5} (\nabla {\bf u}_s) \cdot {\bf q}_s + \frac{5}{2} \frac{kp_s}{m_s} \nabla T_s - \frac{e_s}{m_s} {\bf q}_s \times {\bf B} = \frac{\delta{\bf q}_s'}{\delta t}, $$

(64)

where ∇u _s denotes a tensor, and the right-hand side again denotes the contribution from the collision term.

The corresponding collision terms have only been evaluated when flow velocity differences are small compared with species thermal speeds, and are given by Burgers (1969). If, furthermore, it is assumed that temperature differences are small, the collision terms on the right hand-side of the equations of momentum and energy conservation (50) and (63) become

$$ \frac{\delta {\bf M}_s}{\delta t} = m_s n_s \sum_t \nu_{st} ({\bf u}_t - {\bf u}_s) + \frac{3}{5} \sum_t \nu_{st} \frac{\mu_{st}}{kT_{st}} \left( {\bf q}_s - \frac{m_s n_s}{m_t n_t} {\bf q}_t \right) $$

(65)

$$ \frac{\delta E_s}{\delta t} =m_s n_s \sum_t \frac{\nu_{st}} {m_s+m_t} 3 k (T_t-T_s) $$

(66)

With the same approximation, δq _s ′/δt is given in a concise form by Schunk (1977). As mentioned at the end of Sect. 5.3.1 the frictional heating is absent in (66), because of the approximations that had to be made when evaluating the collision terms.

By allowing for heat conduction in the formulation, not only the description of the energy balance in the corona and solar wind is improved, but also the description of forces, by allowing for thermal diffusion. Moreover, by having a separate equation (64) for the heat flux density, this set of equations displays a self-consistent transition to collisionless flow: classical heat conduction (${\bf q}_s \propto -\nabla T_s$) is only retained in the collision-dominated limit of small flow speeds when the term proportional to ∇T _s dominates the left-hand side of (64), and there is thus no need for some ad hoc transition from classical to non-classical heat conduction in this equation set.

Although heat conduction and thermal diffusion are accounted for in this description, they are not accounted for well. In the collision-dominated limit with a small temperature gradient, the electron heat conduction coefficient κ ^′ _e [see (44)] obtained from this approximation is less than half of the correct value obtained, e.g., by Spitzer and Härm (1953), based on a numerical solution of the Boltzmann equation with the Fokker–Planck collision term. Thermal diffusion is off by a similarly large factor. The reason for the discrepancy lies in the choice of a first order (in c _s ) correction term in (61). In a collision-dominated region with a small temperature gradient, the actual electron VDF for small values of c _e is rather proportional to $c_e^3$ or $c_e^4$ (Spitzer and Härm 1953). Since the Coulomb cross section is so sensitive to the relative velocity [(6)], this inaccuracy leads to large errors in the collision integrals. To improve the description of Coulomb collisions it is therefore necessary to include higher-order terms in c _s in ϕ _s . Note also that the term of order $c_s^3$ in (61) is not really a third order term as it does not have an independent coefficient; recall that it was merely introduced to ensure that u _s remained the mean flow velocity.

The second sum in (65) accounts for thermal diffusion, which is important for minor ions (and He ions) in the transition region and inner corona (Nakada 1969), but also in coronal loops. In regions with steep temperature gradients, heavy, ionized particles will collide more often with “cold” electrons and protons streaming upwards than with the (relatively) hotter particles streaming down the temperature gradient, the net effect being that the heavy particles feel a strong upward force which is even stronger than gravity in the transition region. In a pure electron-proton plasma with no currents the thermal diffusion term merely modifies the electric field, but has no impact on the electron-proton flow.

5.3.3 13-Moment Approximation

The next level of approximation is to add a second order (in c _s ) term to ϕ _s (Grad 1958),

$$ \phi_s = - \frac{m_s}{kT_sP_s} \left( 1 - \frac{m_s c_s^2} {5kT_s} \right) {\bf q}_s \cdot {\bf c}_s + \frac{\tau_{sij}} {2kT_sP_s} c_i c_j, $$

(67)

where τ _sij are components of the stress tensor $\underline{\underline{{\tau_s}}}$, given as

$$ \tau_{sij} \equiv m_s \int\left( c_{si} c_{sj} - \frac{1}{3} c_s^2 \delta_{ij}\right) f_s d{\bf c}_s , $$

(68)

and a sum over indices i and j is implied in (67). Because τ _sij  = τ _sji and τ _sii  = 0 it has 5 independent components, hence this is a 13-moment approximation. In addition to the pressure gradient of the momentum equation (50), we now must add the divergence of $\underline{\underline{{\tau_s}}}$,

$$ m_s n_s \frac{D_s {\bf u}_s}{Dt} + \nabla p_s + \nabla \cdot \underline{\underline{{\tau_s}}} - m_s n_s {\bf g} - e_s n_s ( {\bf E} + {\bf u}_s \times {\bf B}) = \frac{\delta {\bf M}_s} {\delta t}, $$

(69)

where $\nabla \cdot \underline{\underline{{\tau_s}}} \equiv {\bf e}_k \cdot \partial(\tau_{sij} {\bf e}_i {\bf e}_j)/\partial x_k$. The equation for the thermal energy reads

$$ \frac{D_s}{Dt} \left( \frac{3}{2} p_s \right) + \frac{5}{2} p_s \nabla \cdot {\bf u}_s + \nabla \cdot {\bf q}_s + \tau_{sij} \frac{\partial {u_{sj}}}{\partial {x_i}} = \frac{\delta E_s} {\delta t}. $$

(70)

Again collision terms are hard to evaluate without making severe approximations (except for the non-relevant case of Maxwell molecule interactions). The momentum and energy collision terms in the limit of small flow speed differences and small relative temperature differences are identical to (65) and (66) of the 8-moment description. Zamlutti (1998) has proposed approximate collision terms for the 13-moment approximation also to be used for large flow speed differences. However, these collision terms have only been obtained from the Boltzmann collision integral, and are not directly applicable to charged particle interactions which are dominated by small-angle collisions (described by the Fokker–Planck collision term).

The equations for $\underline{\underline{{\tau_s}}}$ and q _s are given, e.g., by Schunk (1977). By including the stress tensor, this formalism includes viscous effects in addition to heat conduction, and in a collision-dominated gas the Navier-Stokes equations may be retrieved (see Gombosi 1994).

By including a higher-order component (the stress tensor), this formalism should better describe a weakly collisional solar wind plasma, in which departure from a Maxwellian will be larger. However, since the collision terms have only been evaluated in a limit which requires a strongly collisional plasma, it is far from obvious that these equations will provide a better description of a weakly collisional plasma. And in the fully collisionless solar wind departures from a Maxwellian may become large, so that there is no guarantee that even a second (or higher) order approximation to ϕ _s is any better than the first order approximation. In that case, it may be advantageous to switch to gyrotropic equations in which the VDF is expanded about a bi-Maxwellian, which we shall discuss in the next section.

Note also that the assumption (67) still does not accurately describe the electron VDF in a collision-dominated plasma, where, as we have discussed, ϕ _s should rather be of third or fourth order in c _s . Hence it is expected that the errors in the collision terms are still large; since δM _s /δt remains unchanged from the 8-moment approximation, thermal diffusion will still not be accurately described.

To improve the description of collisions, one should at least include third order terms in ϕ _s . This results in the so-called 20-moment approximation (see Gombosi 1994), in which 20 coupled partial differential equations have to be solved. Needless to say, these are so complex that they seem to be of little practical use.

The higher-order moments, such as the stress tensor, are usually of little physical importance in themselves for the solar wind system. Mostly we are concerned with the basic properties, such as the mass flux, flow speed, and energy flux. The higher-order moments are only of interest as they affect the low-order moments, either through the closure assumption or because they modify the collision terms. One may therefore ask whether solving higher-order equations is worthwhile, given the tremendous effort required to obtain a numerical solution to them, keeping in mind that we still have to use highly simplified collision terms.

5.3.4 8-Moment Approximation with Improved Coulomb Collision Terms

If describing collisions reasonably correctly is important, instead of trying to solve the almost intractable 20-moment set of equations one may instead make a simple assumption for ϕ _S which still captures the essence of Coulomb collisions. Killie et al. (2004) suggested that the 8-moment approximation (61) be replaced by an 8-moment approximation which is third order in c _s ,

$$ \phi_s = - \frac{m_s c_s^2}{5k^2T_s^2P_s} \left( 1 - \frac{m_s c_s^2}{7kT_s} \right) {\bf q}_s \cdot {\bf c}_s, $$

(71)

where the fifth order term and the coefficients have been chosen such that u _s and q _s remain the flow speed and heat flux density, respectively. With this, still 8-moment, approximation, it turns out that the left-hand sides of the transport equations are identical to those of the original 8-moment approximation, given by (20), (50), (63), and (64). However, the collision terms δM _s /δt and δq _s ′/δt will change. Collision terms still have to be evaluated in the “semi-linear” approximation of small flow speed differences. In the linear limit of small relative temperature differences as well, one obtains

$$ \begin{aligned} \frac{\delta{\bf M}_s}{\delta t} &\, =\, n_s m_s \sum_{t} \nu_{st} ({\bf u}_t - {\bf u}_s) + \sum_{t}\nu_{st} \frac{3}{5} \frac{\mu_{st}}{kT_{st}} \left[{\bf q}_s\left( 1 - \frac{5}{7} \frac{m_t}{m_s+m_t} \right) \right. \\ &\quad \left. - {\bf q}_t \frac{m_s n_s}{m_t n_t} \left( 1 - \frac{5}{7} \frac{m_s}{m_s+m_t} \right) \right]. \end{aligned} $$

(72)

The modified terms in the second sum [compared with (65)] means that thermal diffusion is now reduced (when the heat flux remains unchanged). The heat flux collision term, δq _s ′/δt, also changes substantially.

In the collision-dominated limit with small temperature gradients, the new collision terms lead to a heat conduction coefficient κ _e ′ that does not deviate by more than 20% from the correct value (see the expression (44) of the Spitzer-Härm heat flux); note also that the original 8-moment assumption described in Sect. 5.3.2 leads to a more than factor 2 deviation. The coefficients of thermal diffusion agree within approximately 30%. Hence this equation set should be well suited to the part of the solar wind strongly influenced by collisions, while the equations also contain a self-consistent transition to collisionless flow. In the collisionless regime the new equations are identical to the original 8-moment equations.

This shows how it is possible to improve the model without greatly increasing the complexity, by using a carefully chosen ϕ _s . As a formal derivation it seems strange to start the Taylor expansion at third order—what happened to the first and second order terms? A formal expansion to third order implies solving the formidable 20-moment set of equations, as well as deriving the appropriate collision terms. Because the complexity increases drastically when the order is increased, it is imperative to keep the order as low as possible. Hence, in this case it is better to skip the first two terms in the Taylor expansion altogether, knowing that they must be negligible in a collision-dominated plasma.

Note that this approximation is tailored to a fully ionized, collision-dominated plasma. In a weakly ionized or neutral gas, where Coulomb collisions do not dominate, this equation set should offer no improvement over the original 8-moment set.

5.4 Models for the Supersonic Solar Wind

Since collisions quickly cease to be important as the solar wind expands and becomes supersonic (except for the electrons), describing collisions is no longer essential. Rather, the fluid description should describe the (nearly) collisionless evolution reasonably well, and be able to accommodate effects of wave-particle interactions. Since collisions are nearly absent, the departure from a Maxwellian [ϕ _s in (45)] may become large. One may use higher-order fluid models, described in the previous section, but it can be argued that the cost of doing that may be prohibitive. Moreover, if the departure from a Maxwellian is large, it may no longer make sense to expand about a Maxwellian and it may be advantageous to consider expansions about other zero-order assumptions for the VDF.

The magnetic field in the corona and inner solar wind is strong (it is a low-β plasma) and the rapid gyration of charged particles in the magnetic field implies that the VDF should be nearly gyrotropic (i.e., isotropic in the plane perpendicular to B), at least in the rest frame of the plasma (the plasma as a whole may undergo E × B drift). If collisions are few, the motion along and perpendicular to B may become partly decoupled. A more accurate description of such a magnetized plasma will be provided if different “temperatures” are considered along and perpendicular to the magnetic field. Moreover, as also discussed in Sect. 4.4, observations of ions in the corona by the UVCS instrument on the SOHO satellite show that even in this region, where collisions should be frequent, the temperature along the line of sight (assumed to be nearly perpendicular to B in the corona) is much higher than the parallel temperature (Kohl et al. 1997; Cranmer et al. 1999; Antonucci et al. 2000; Zangrilli et al. 2002). To capture such effects, without going to prohibitively high order in ϕ _s , the moment expansion about a zeroth order velocity distribution has to allow for temperature anisotropies and to be sufficiently simple that collision terms can still be developed (since they should not be negligible in the corona).

In the gyrotropic formulation the VDF is expanded about a bi-Maxwellian,

$$ f_s^{\rm bM} \equiv n_s \frac{m_s}{2\pi kT_{s\perp}} \sqrt{\frac{m_s}{2\pi k T_{s\parallel}}} \exp\left[ - m_s \left( \frac{c_{s\perp}^2}{2kT_{s\perp}} + \frac{c_{s\parallel}^2} {2kT_{s\parallel}} \right) \right]. $$

(73)

Here ${\bf c}_{s\parallel}$ and ${\bf c}_{s\perp}$ are the components of the peculiar velocity c _s  ≡ v − u _s parallel and perpendicular to the magnetic field, ${\bf c}_{s\parallel} \equiv {\bf e}_3 {\bf e}_3 \cdot {\bf c}_s$ and ${\bf c}_{s\perp} \equiv (\underline{\underline{{I}}} - {\bf e}_3 {\bf e}_3) \cdot {\bf c}_s$, where $\underline{\underline{{I}}} \equiv {\bf e}_1 {\bf e}_1 + {\bf e}_2 {\bf e}_2 + {\bf e}_3 {\bf e}_3$ is the unit tensor. The unit vectors e ₁, e ₂, and e ₃ are orthogonal, with e ₃ chosen to be parallel to B. The expansion then proceeds as in Sect. 5.3, assuming the full VDF to be

$$ f_s = f_s^{\rm bM} (1 + \chi_s) $$

(74)

where χ _s is assumed to be small (at least when evaluating collision terms).

This choice has the advantage that the base function $f_s^{\rm bM}$ is quite simple, and that it reduces to the isotropic Maxwellian (46) when $T_{s\parallel} = T_{s\perp} = T_s$. Hence it has the potential that a reasonably correct VDF may be retrieved in the collision-dominated limit, and thus that it can reproduce, reasonably well, the transport coefficients of classical transport theory. Hence fluid equations based on (74), with a proper choice for χ _s , should give a reasonable description of the solar wind flow all the way from the chromosphere and into interplanetary space.

5.4.1 Gyrotropic Equations with χ _s  = 0

Again, the simplest possible choice is just to use (73) as the closure assumption for the VDF. In addition to the continuity equation (20), this leads to five equations for the expansion coefficients u _s , $T_{s\parallel}$, and $T_{s\perp}$. In compact notation these may be written

$$ n_s m_s \frac{D{\bf u}_s}{Dt} + \nabla \cdot \underline{\underline{{p_s}}} - n_s m_s {\bf g} - n_s e_s ({\bf E} + {\bf u}_s \times {\bf B}) = \frac{\delta{\bf M}_s}{\delta t} $$

(75)

$$ \frac{D \underline{\underline{{p_s}}}}{Dt}+ \underline{\underline{{p_s}}} (\nabla \cdot {\bf u}_s) + 2( \underline{\underline{{p_s}}} \cdot \nabla) {\bf u}_s = \frac{\delta \underline{\underline{{E_s}}}}{\delta t}, $$

(76)

where $\underline{\underline{{p_s}}} = p_{s\perp} (\underline{\underline{{I}}}-{\bf e}_3 {\bf e}_3) + p_{s\parallel} {\bf e}_3 {\bf e}_3$ is the partial pressure tensor of species s, with $p_{s\parallel(\perp)} = n_s k T_{s\parallel(\perp)}$. Equivalently, expressed in terms of $p_{s\parallel}$ and $p_{s\perp}$ they may be written,

$$ n_s m_s \frac{D{\bf u}_s}{Dt} + \nabla_\perp p_{s\perp} + \nabla_\parallel p_{s\parallel} + (p_{s\parallel} - p_{s\perp}) \nabla \cdot ({\bf e}_3 {\bf e}_3) \\ - n_s m_s {\bf g} - n_s e_s ({\bf E} + {\bf u}_s \times {\bf B}) = \frac{\delta{\bf M}_s}{\delta t} $$

(77)

$$ \frac{D p_{s\parallel}}{Dt} + p_{s\parallel} (\nabla \cdot {\bf u}_s + 2 \nabla_\parallel \cdot {\bf u}_s) =\frac{\delta E_{s\parallel}}{\delta t} $$

(78)

$$ \frac{D p_{s\perp}}{Dt} + p_{s\perp} (\nabla \cdot {\bf u}_s + \nabla_\perp \cdot {\bf u}_s) = \frac{\delta E_{s\perp}}{\delta t} $$

(79)

where $\nabla_\parallel \equiv {\bf e}_3 {\bf e}_3 \cdot \nabla$ and $\nabla_\perp \equiv (\underline{\underline{{I}}} - {\bf e}_3 {\bf e}_3) \cdot \nabla$.

This set of equations was first derived by Chew et al. (1956), although they did not make the assumption (73) explicitly, and they used the ideal MHD approximation E = −u _s  × B and Maxwell’s equations to cast the Lorentz force term in (75) in terms of u _s and B.

With χ _s  = 0 the VDF is sufficiently simple that the collision terms have been calculated without having to make further approximations, even for Coulomb collisions and even valid for (arbitrarily) large flow speed differences and temperature anisotropies (see, Barakat and Schunk 1981). A closed form for the Coulomb collisions terms for a bi-Maxwellian plasma drifting in the direction parallel to the magnetic field has been recently found by Hellinger and Trávnicek (2009). The terms are sufficiently complicated that they will not be reprinted here, though.

Solving the full 3D set of equations, even this fairly simple set, can be a challenging task. Often it is sufficient to consider the flow in only one dimension. For instance, the steady state fast solar wind from polar coronal holes may to a good approximation be regarded as a 1D problem, in which the plasma flows along field lines. In that case the only effect of the magnetic field is to specify the flow geometry. Moreover, far from the Sun the solar wind plasma pressure is substantially larger than the magnetic pressure, in which case the magnetic field plays no role in the wind expansion, which may be regarded as approximately radial. In such cases we may imagine that the flow occurs along a flow tube which has a cross section of, say, $A_0=1\,{{\hbox{m}}}^2$ at the solar “surface” and which is centered on a radial magnetic field line. At a distance r from the center of the Sun the flux tube cross section is A(r). Equations (20) and (77)–(79) then become

$$ \frac{\partial {n_s}}{\partial {t}} + \frac{1}{A} \frac{\partial {(n_su_s A)}}{\partial {r}}= 0 $$

(80)

$$ \frac{\partial {u_s}}{\partial {t}} + u_s \frac{\partial {u_s}}{\partial {r}} + \frac{k}{m_s n_s} \frac{\partial {(n_s T_{s\parallel})}}{\partial {r}} + \frac{1}{A} \frac{d {A}}{d {r}} \frac{k}{m_s} (T_{s\parallel}-T_{s\perp}) - g - \frac{e_s}{m_s} E = \frac{1}{m_s n_s}\frac{\delta M_s}{\delta t} $$

(81)

$$ \frac{\partial {T_{s\parallel}}}{\partial {t}} + u_s \frac{\partial {T_{s\parallel}}}{\partial {r}} + 2 T_{s\parallel} \frac{\partial {u_s}}{\partial {r}}= \frac{1}{n_s k} \frac{\delta E_{s\parallel}}{\delta t} $$

(82)

$$ \frac{\partial {T_{s\perp}}}{\partial {t}} + u_s \frac{\partial {T_{s\perp}}}{\partial {r}}+ \frac{1}{A} \frac{d {A}} {d {r}} u T_{s\perp} = \frac{1}{n_s k} \frac{\delta E_{s\perp}}{\delta t}. $$

(83)

A purely radially expanding outflow (i.e., spherically symmetric outflow) corresponds to $A(r) = A_0 (r/R_S)^2.$

Note from (81) that $T_{s\perp} > T_{s\parallel}$ leads to acceleration of the wind. This is a consequence of the Lorentz force: when the magnetic field changes slowly with radial distance (as we assume), the magnetic moment is conserved, leading to conversion of perpendicular motion into parallel motion and hence acceleration. With an isotropic VDF as in Sect. 5.3.1, we are not able to include this effect of the Lorentz force. Note also that in a steady state with no collisions, (82) and (83) simplify to

$$ \frac{d {(u_s^2 T_{s\parallel})}}{d {r}} =0 $$

(84)

$$ \frac{d {(A T_{s\perp})}}{d {r}} =0, $$

(85)

stating that as the wind is accelerated, $T_{s\parallel}$ will decrease with distance, while $T_{s\perp}$ decreases as the flow tube expands. If the wind is accelerated rapidly near the Sun, and rapidly compared to the flow tube expansion, we expect that $T_{s\parallel} < T_{s\perp}$. On the other hand, in the outer solar wind, where $du_s/dr\,{\approx}\,0$, eventually $T_{s\parallel} \gg T_{s\perp}$. This illustrates that, in the absence of collisions or wave-particle interactions, anisotropies must arise in the solar wind, showing why it is necessary to go beyond the isotropic formulation of Sect. 5.3 describing the supersonic solar wind far from the Sun.

This equation set is simple, collision terms valid for any flow conditions are available, and they capture important effects of the Lorentz force. Their main limitation is that they do not allow for heat fluxes. Although heat conduction is not as critical from the corona and outwards as it is in the transition region, it nevertheless plays a role, particularly for the electrons near the Sun.

Another shortcoming is that the temperature evolution they predict cannot be entirely correct. Imagine that in the corona $T_{s\perp} \gg T_{s\parallel}$, implying that some process (e.g., cyclotron waves) has heated the particles perpendicularly to B, and also supported by observations. Neglecting collisions and interaction with waves, magnetic moment is conserved, and the large perpendicular motion in the corona must be translated into parallel motion in the solar wind. However, all of the perpendicular motion cannot be translated into the flow speed u _s (that is, wind acceleration): a high $T_{s\perp}$ implies that the VDF is very broad in the perpendicular direction, with some particles having a small perpendicular velocity and others a large velocity. Since the magnetic field changes only the direction of the velocity vector, the magnetic field expansion should lead to a VDF that is similarly broad in the parallel direction. If one solves the collisionless Boltzmann equation, mapping the VDF from the corona and outwards, one would thus find that the true VDF in the solar wind is narrow in the perpendicular direction, but very broad in the parallel direction, which means that $T_{s\parallel}$ in the solar wind must become large, as shown by the collisionless kinetic models (Sect. 4.4). This effect is not captured at all with χ _s  = 0. In fact, (84) predicts the opposite effect: $T_{s\parallel}$ should decrease monotonically outwards as the wind is accelerated, and there is no coupling that converts a high $T_{s\perp}$ in the corona into a high $T_{s\parallel}$ in the solar wind. The remedy for this flaw is also to allow for heat flux by having a nonzero χ _s .

5.4.2 The 16-Moment Approximation

Hydrodynamic equations based on (74) that allow for heat fluxes as well as stresses were first developed by Oraevskii et al. (1968) for a collisionless plasma. Their assumption for χ _s is, when written in the form used by Barakat and Schunk (1982),

$$ \begin{aligned} \chi_s &= \frac{\beta_{s\perp}}{2 m_s n_s}\left[ \beta_{s\perp} (c_{s1}^2-c_{s2}^2) \underline{\underline{{\tau_{s}}}} : {\bf e}_{1} {\bf e}_{1} + 2\beta_{s\perp} \underline{\underline{{\tau_{s}}}} : {\bf e}_{1} {\bf e}_{2} c_{s1} c_{s2} + 2\beta_{s\parallel} \underline{\underline{{\tau_{s}}}} : {\bf c}_{s\perp} {\bf c}_{s\parallel} \right] \\ & \quad - \frac{\beta_{s\perp}^2}{m_s n_s} \left( 1 - \frac{\beta_{s\perp} c_{s\perp}^2}{4} \right) {\bf q}_{s}^\perp \cdot {\bf c}_{s\perp} - \frac{\beta_{s\perp} \beta_{s\parallel}}{m_s n_s} \left( 1 - \frac{\beta_{s\perp} c_{s\perp}^2}{2} \right) {\bf q}_s^\perp \cdot {\bf c}_{s\parallel} \\ &\quad -\frac{\beta_{s\parallel}^2}{2m_s n_s} \left( 1 - \frac{\beta_{s\parallel} c_{s\parallel}^2}{3} \right) {\bf q}_s^\parallel \cdot {\bf c}_{s\parallel} - \frac{\beta_{s\perp} \beta_{s\parallel}}{2 m_s n_s}( 1 - \beta_{s\parallel} c_{s\parallel}^2) {\bf q}_s^\parallel \cdot {\bf c}_{s\perp} \end{aligned} $$

(86)

where $\beta_{s\parallel} \equiv m_s/(kT_{s\parallel})$ and $\beta_{s\perp} \equiv m_s/(kT_{s\perp})$ and we have introduced the tensor contraction $\underline{\underline{{A}}} : \underline{\underline{{B}}} \equiv A_{ij} B_{ji}$. In this approximation there are 16 moments of the VDF to be solved for. Note that there are now two heat flux densities: ${\bf q}_{s}^\parallel$ denote the transport of parallel thermal energy while ${\bf q}_{s}^\perp$ denote the transport of perpendicular thermal energy. Since $T_{s\parallel}$ and $T_{s\perp}$ may be very different, the conduction of the two forms of thermal energy may also be different.

With this assumption, the left-hand side of the momentum equation acquires an additional stress tensor term [compared with the χ _s  = 0 case (77)] while the equations for $p_{s\parallel}$ and $p_{s\perp}$ get additional terms from the heat flux vectors and the stress tensor:

$$ \begin{aligned} & n_s m_s \frac{D{\bf u}_s}{Dt} + \nabla_\perp p_{s\perp} + \nabla_\parallel p_{s\parallel} + (p_{s\parallel} - p_{s\perp}) \nabla \cdot ({\bf e}_3 {\bf e}_3) + \nabla \cdot \underline{\underline{{\tau_s}}}\\ &\quad - n_s m_s{\bf g} - n_s e_s ({\bf E} + {\bf u}_s \times {\bf B}) = \frac{\delta{\bf M}_s}{\delta t} \end{aligned} $$

(87)

$$ \begin{aligned} & \frac{Dp_{s\parallel}}{Dt} + p_{s\parallel} (\nabla \cdot {\bf u}_s + 2 \nabla_\parallel \cdot {\bf u}_s) + 2 {\bf e}_3 {\bf e}_3 : ( \underline{\underline{{\tau_s}}} \cdot \nabla {\bf u}_s ) + \nabla \cdot {\bf q}_s^\parallel \\ &\quad - \underline{\underline{{\tau_s}}} : \frac{D({\bf e}_3 {\bf e}_3)}{Dt} -\underline{\underline{\underline{{Q_s}}}}\,\vdots\,\nabla({\bf e}_3 {\bf e}_3) =\frac{\delta E_{s\parallel}}{\delta t} \end{aligned} $$

(88)

$$ \begin{aligned} & \frac{D p_{s\perp}}{Dt} + p_{s\perp} (\nabla \cdot {\bf u}_s + \nabla_\perp \cdot {\bf u}_s) + (\underline{\underline{{I}}} - {\bf e}_3 {\bf e}_3 ) : (\underline{\underline{{\tau_s}}} \cdot \nabla {\bf u}_s) + \nabla \cdot {\bf q}_s^\perp \\ &\quad + \frac{1}{2} \underline{\underline{{\tau_s}}} : \frac{D({\bf e}_3 {\bf e}_3)}{Dt} + \frac{1}{2} \underline{\underline{\underline{{Q_s}}}}\,\vdots\,\nabla({\bf e}_3 {\bf e}_3) = \frac{\delta E_{s\perp}}{\delta t} \end{aligned} $$

(89)

where $\underline{\underline{\underline{{Q_s}}}}$ is a rank 3 tensor given as

$$ \underline{\underline{\underline{{Q_s}}}} = {\bf q}_s^\parallel {\bf e}_3 {\bf e}_3 + e_3 {\bf q}_s^\parallel {\bf e}_3 + {\bf e}_3 {\bf e}_3 {\bf q}_s^\parallel - 2 {\bf q}_{s\parallel}^\parallel {\bf e}_3 {\bf e}_3 $$

(90)

and $\underline{\underline{\underline{{C}}}}\,\vdots\,\underline{\underline{\underline{{D}}}} \equiv C_{ijk} D_{kji}$. The equations for the stress tensor $\underline{\underline{{\tau_s}}}$ and the heat flow vectors are given by Barakat and Schunk (1982) and will not be listed here.

The corresponding Coulomb collision terms have been derived by Chodura and Pohl (1971), again in the limit of small flow speed differences, and are too complex to be listed here. Collision terms for other types of interactions were later derived by Demars and Schunk (1979).

If we limit ourselves to one spatial dimension, with flow along a radial magnetic field with a specified flow tube area A(r), we retrieve the equations (80)–(83), but with additional source terms for $T_{s\parallel}$ and $T_{s\perp}$, and additional equations for $q_s^\parallel$ and $q_s^\perp$ (the radial flow of parallel and perpendicular thermal energy, respectively):

$$ \frac{\partial {T_{s\parallel}}}{\partial {t}} + u_s \frac{\partial {T_{s\parallel}}}{\partial {r}}+ 2 T_{s\parallel} \frac{\partial {u_s}}{\partial {r}} + \frac{1}{n_s k} \frac{\partial {q_s^\parallel}}{\partial {r}} + \frac{1}{A} \frac{d {A}}{d {r}} \frac{q_{s}^\parallel}{n_s k} - \frac{2}{A} \frac{d {A}}{d {r}} \frac{q_s^\perp}{n_s k} =\frac{1}{n_s k} \frac{\delta E_{s\parallel}}{\delta t} $$

(91)

$$ \frac{\partial {T_{s\perp}}}{\partial {t}}+ u_s \frac{\partial {T_{s\perp}}}{\partial {r}}+ \frac{1}{A} \frac{d {A}} {d {r}} u T_{s\perp} + \frac{1}{n_s k} \frac{\partial {q_s^\perp}}{\partial {r}} + \frac{2}{A} \frac{d {A}} {d {r}} \frac{q_s^\perp}{n_s k} = \frac{1}{n_s k} \frac{\delta E_{s\perp}}{\delta t} $$

(92)

$$ \frac{\partial {q_{s}^\parallel}}{\partial {t}} + u_s \frac{\partial {q_{s}^\parallel}}{\partial {r}}+ 4 q_{s}^\parallel \frac{\partial {u_s}}{\partial {r}} + u_s q_{s}^\parallel \frac{1}{A} \frac{d {A}}{d {r}} + 3 \frac{k^2n_sT_{s\parallel}}{m_s} \frac{\partial {T_{s\parallel}}}{\partial {r}} = \frac{\delta q_{s}^\parallel}{\delta t}' $$

(93)

$$ \frac{\partial {q_{s}^\perp}}{\partial {t}}+u_s \frac{\partial {q_{s}^\perp}}{\partial {r}}+ 2 q_{s}^\perp \frac{\partial {u_s}} {\partial {r}} + 2 u_s q_{s}^\perp \frac{1}{A} \frac{d {A}}{d {r}} + \frac{k^2n_sT_{s\parallel}}{m_s} \frac{\partial {T_{s\perp}}}{\partial {r}}\\ + \frac{1}{A} \frac{d {A}}{d {r}} \frac{k^2n_sT_{s\perp}} {m_s} (T_{s\parallel}-T_{s\perp}) = \frac{\delta q_{s}^\perp} {\delta t}'. $$

(94)

Comparing (82) and (91), note that $q_s^\perp$ has now become a “source” for $T_{s\parallel}$. Because of this term, if $T_{s\perp} \gg T_{s\parallel}$ in the corona, magnetic moment conservation leads to an increasing $T_{s\parallel}$ with increasing radial distance. This effect is also demonstrated in numerical solar wind solutions using this equation set (Lyngdal Olsen and Leer 1999), and it is in agreement with solutions to the collisionless Boltzmann equation, and the opposite of the behavior predicted by (84). Hence by including heat fluxes, also the description of the collisionless flow is improved.

In situ observations of the proton VDF between 0.3 and 1 AU by the Helios spacecraft (Marsch et al. 1982) show that in most cases $T_{p\parallel}/T_{p\perp}>1$, and that the ratio tends to increase with radial distance from the Sun, consistent with the effect of magnetic moment conservation described above. However, there are certainly cases where $T_{p\parallel}/T_{p\perp} \le 1$, particularly for the fast solar wind at small heliocentric distances. This shows that the proton magnetic moment is not strictly conserved. Since Coulomb collisions should be quite unimportant for protons so far from the Sun, it indicates that wave-particle interactions must take place in the solar wind. We should also point out that if $T_{p\parallel} \gg T_{p\perp}$ the firehose plasma instability may be triggered, which would tend to lower the temperature anisotropy. Such effects are not included in the models considered here, but can certainly be included by appropriate choices for the collision terms on the right-hand side of the equations (which would then account for both Coulomb collisions and wave-particle interactions). To determine the magnitude of these plasma processes, required to reproduce the observed VDF, it is important that the fluid models describe reasonably well the effect of the expanding magnetic field and the wind acceleration on the temperature anisotropy.

With the actual temperature anisotropies and radial variation of temperatures measured by Helios, Holzer et al. (1986) found that the anisotropy of the pressure tensor had little influence on the solar wind acceleration and heating of the plasma, compared with the pressure gradient force and the pressure work term, respectively.

Although the description of the collisionless solar wind is improved by the 16-moment approach, compared to the χ _s  = 0 case, the description of the heat conduction in a collision-dominated plasma is not accurate. The heat conduction terms in (86) are still first order in c _s , as it is in the isotropic 13-moment approximation (67). In a multi-fluid model of a collision-dominated plasma the electron heat conduction will thus not be described much better than in the simple 8-moment approximation of Sect. 5.3.2.

5.4.3 Gyrotropic Equations with Improved Collision Terms

To improve the description of Coulomb collisions it is possible to repeat the approach of Sect. 5.3.4, namely by letting χ _s be of third order in c _s instead of first order. Janse et al. (2005) developed gyrotropic transport equations assuming

$$ \chi_s = \alpha_{s \perp} {\bf q}_s \cdot {\bf c}_{s \perp} c^2_s \left( 1 + \gamma_{s \perp} c^2_s \right) + \alpha_{s \parallel} {\bf q}_s \cdot {\bf c}_{s \parallel} c^2_s \left( 1 + \gamma_{s \parallel} c^2_s\right), $$

(95)

where the coefficients are chosen such that u _s and q _s are the flow velocity and the heat flux density, leading to

$$ \gamma_{s \perp} = -\frac{m_s}{kT_{s \perp}}\frac{4 T_{s \perp}^2 + T_{s \perp} T_{s \parallel}}{24 T_{s \perp}^2 + 8 T_{s \perp} T_{s \parallel} + 3T_{s \parallel}^2} $$

(96)

$$ \gamma_{s \parallel} = -\frac{m_s}{kT_{s \perp}}\frac{2T_{s \perp}^2+3T_{s \perp}T_{s \parallel}} {8T_{s \perp}^2+12T_{s \perp}T_{s \parallel}+15T_{s \parallel}^2} $$

(97)

$$ \alpha_{s \perp} = - \frac{m_s^2}{n_s k^3} \frac{1}{T_{s \perp}}\frac{24T_{s \perp}^2+8T_{s \perp}T_{s \parallel}+3T_{s \parallel}^2}{96T_{s \perp}^4+48T_{s \perp}^3T_{s \parallel}+4T_{s \perp}^2T_{s \parallel}^2+24T_{s \perp}T_{s \parallel}^3+3T_{s \parallel}^4} $$

(98)

$$ \alpha_{s \parallel} = - \frac{m^2_s}{n_s k^3} \frac{1}{T_{s \parallel}}\frac{8 T_{s \perp}^2+12T_{s \perp}T_{s \parallel} + 15 T_{s \parallel}^2}{16 T_{s \perp}^4+48 T_{s \perp}^3 T_{s \parallel} + 6 T_{s \perp}^2 T_{s \parallel}^2 + 60 T_{s \perp} T_{s \parallel}^3 + 45 T_{s \parallel}^4}. $$

(99)

Note that, compared with (86), the stress tensor has been omitted and there is only one heat flux vector.

For the case of one-dimensional, radial flow, the resulting momentum and energy equations become formally identical to (81), (91), and (92), while the equation for the heat flow becomes

$$ \begin{aligned} \frac{\partial {q_s}}{\partial {t}} + u_s \frac{\partial {q_s}} {\partial {r}} + 2 q_{s}^\parallel \frac{\partial {u_s}} {\partial {r}} + \frac{1}{2} q_{s}^\parallel u_s \frac{1}{A} \frac{d {A}}{d {r}}+ 2 q_{s}^\perp \frac{\partial {u_s}}{\partial {r}} + 2 q_{s}^\perp u_s \frac{1}{A} \frac{d {A}}{d {r}}\\ +\frac{k^2 n_s T_{s\parallel}}{m_s} ( T_{s\parallel} - T_{s\perp} ) = \frac{\delta q_s'}{\delta t}, \end{aligned} $$

(100)

where

$$ q_{s}^\parallel = 30 q_s \frac{T_{s \parallel}^3(4 T_{s \perp}+3 T_{s \parallel})}{16 T_{s \perp}^4+48T_{s \perp}^3T_{s \parallel}+6T_{s \perp}^2T_{s \parallel}^2+60T_{s \perp}T_{s \parallel}^3+45T_{s \parallel}^4} $$

(101)

$$ q_{s}^\perp = 2 q_s \frac{T_{s \perp}^2 (8 T_{s \perp}^2+24T_{s \perp} T_{s \parallel} + 3 T_{s \parallel}^2)}{16T_{s \perp}^4 + 48 T_{s \perp}^3 T_{s \parallel} + 6 T_{s \perp}^2 T_{s \parallel}^2 + 60 T_{s \perp} T_{s \parallel}^3 + 45 T_{s \parallel}^4} $$

(102)

still represent the flow of parallel and perpendicular thermal energy in the radial direction, satisfying $q_s = q_s^\parallel/2 + q_s^\perp$.

In the isotropic limit $T_{s\parallel} = T_{s\perp}$ these equations simplify to the equations of Sect. 5.3.4. Collision terms for the momentum and heat flow equations consistent with the approximation (95) have not been derived. The corresponding isotropic collision terms (Killie et al. 2004) will not be entirely correct in the outer corona where temperature anisotropies become large. However, the collision terms derived by Chodura and Pohl (1971) are also questionable in this region as they assume that flow speed differences are small.

In the collisionless regime these equations give essentially the same behavior as the 16-moment equations of Sect. 5.4.2 (without the stress tensor). In particular, the conversion of $T_{s\perp}$ in the corona into $T_{s\parallel}$ and wind acceleration is well described also with this equation set.

These equations are thus intended to be used for models of a fully ionized solar wind, extending all the way from the collision-dominated transition region and into the magnetized flow of interplanetary space.

5.4.4 Higher Order Gyrotropic Equations Tailored to Spherically Symmetric Plasma Outflow

Cuperman et al. (1980, 1981) developed transport equations that account for temperature anisotropies, heat conduction, and a higher-order term that leads to non-thermal tails in the assumed VDF. The equations are specifically tailored to one-dimensional (i.e., purely radial), spherically symmetric outflow, where, e.g., heat conduction is strictly in the radial direction. Although temperature anisotropies are included, the expansion is actually about a pure Maxwellian, not a bi-Maxwellian, and is hence based on the same starting point as the equations of Sect. 5.3.

In spherically symmetric outflow, the VDF can only depend on three independent variables: the radial distance r, the speed v = |v|, and $\cos\theta$ where $\theta$ is the angle between the velocity and radius vectors. Hence Cuperman et al. assume that the VDF for species s has the form

$$ f_s(r,c_s,\cos\theta) = f_s^{(0)}(r,c_s) + \sum_{n=0}^2 b_n(r,c_s) f_s^{(0)}(r,c_s) P_n(\cos\theta) $$

(103)

where c _s ≡|v − u _s (r, t)| is still the peculiar speed, $f_s^{(0)}$ is the Maxwellian defined by (46) and P _n are Legendre polynomials.

In addition to the familiar moments n _s , u _s , $T_{s\parallel}$, $T_{s\perp}$, $T=(T_\parallel + 2T_\perp)/3$, and q _s , this expansion introduces the fourth order moment ξ _s (r, t) defined as

$$ \xi_s \equiv \frac{1}{n_s} \int\int\int (c_s\cos\theta)^4 f_s d {\bf v} - 3 \left( \frac{kT_{s\parallel}}{m_s} \right)^2. $$

(104)

ξ _s is thus a moment that can account for tails in f _s . In terms of these, the b _n polynomials in (103) are

$$ b_0(c_s) = c_0 \left( 1 - \frac{2m_sc_s^2}{3kT_s} + \frac{4 m_s^2 c_s^4}{15 (kT_s)^2}\right) $$

(105)

$$ b_1(c_s) = c_1 \sqrt{\frac{m_s}{kT_s}} c_s \left( -1 + \frac{m_s c_s^2}{5kT_s} \right) $$

(106)

$$ b_2(c_s) = c_2 \frac{m_s c_s^2}{kT_s}, $$

(107)

and the coefficients contain the moments of the VDF,

$$ c_0 =\frac{5}{8} \left( \frac{m_s}{kT_s} \right)^2 \xi_s + \frac{4}{3} \left( \frac{T_{s\parallel} - T_{s\perp}} {T_s} \right)^2 $$

(108)

$$ c_1 =\left( \frac{m_s}{kT_s}\right)^{3/2} \frac{q_s}{m_s n_s} $$

(109)

$$ c_2 = \frac{T_{s\parallel} - T_{s\perp}}{3T_s}. $$

(110)

Notice that the b ₁ term, which accounts for heat conduction, is identical to the 8-moment assumption (61).

The resulting momentum equation is identical to the gyrotropic momentum equation (81) when $A(r) \propto r^{-2}$, while the other equations are (omitting subscript s to simplify notation):

$$ \frac{\partial {T_\parallel}}{\partial {t}}\,+\,u \frac{\partial {T_\parallel}}{\partial {r}}+ 2 T_\parallel \frac{\partial {u}}{\partial {r}}+ \frac{6}{5} \frac{1}{nk}\frac{\partial {q}}{\partial {r}} + \frac{4}{5} \frac{1}{nk} \frac{q}{r} = \left( \frac{\delta T_\parallel}{\delta t} \right)_c $$

(111)

$$ \frac{\partial {T_\perp}}{\partial {t}}+ u\frac{\partial {T_\perp}}{\partial {r}} + 2 T_\perp \frac{u}{r} + \frac{2}{5} \frac{1}{nk}\frac{\partial {q}}{\partial {r}} + \frac{8}{5} \frac{1}{nk} \frac{q}{r} = \left( \frac{\delta T_\perp}{\delta t} \right)_c $$

(112)

$$ \frac{\partial {q}}{\partial {t}}+ u \frac{\partial {q}}{\partial {r}} + \frac{2}{r} u q + 4q \frac{\partial {u}}{\partial {r}} + \frac{5}{2} \frac{k^2nT_\parallel}{m} \frac{\partial {T_\parallel}} {\partial {r}}\\ - \frac{10}{3} \frac{k^2 n}{m} (T_\parallel - T_\perp)^2 + \frac{5}{6} m \frac{\partial {(n\xi)}}{\partial {r}} = \left( \frac{\delta q}{\delta t} \right)_c $$

(113)

$$ \begin{aligned} &\frac{\partial {\xi}}{\partial {t}}\,+\,u \frac{\partial {\xi}} {\partial {r}} + 4\xi \frac{\partial {u}}{\partial {r}}+ \frac{36}{5} \frac{k}{n m^2} q \frac{\partial {T}}{\partial {r}} - \frac{16}{5} \frac{k}{nm^2} q \frac{\partial {(T_\parallel - T_\perp)}}{\partial {r}}\\ &\quad + \frac{24}{5} \frac{kT}{m} \frac{\partial }{\partial {r}} \left( \frac{q}{nm} \right) - \frac{24}{5} \frac{k}{m} (T_\parallel - T_\perp) \frac{\partial }{\partial {r}} \left( \frac{q}{nm} \right)\\ &\quad - 8 \frac{k}{(nm)^2} q (T_\parallel - T_\perp) \frac{\partial {n}}{\partial {r}} - \frac{64}{5} \frac{k}{nm^2} \frac{q}{r} (T_\parallel - T_\perp) = \left( \frac{\delta \xi}{\delta t} \right)_c, \end{aligned} $$

(114)

where the right-hand side contain collision terms. Comparing with the 16-moment equations (91)–(92) (when $A(r) \propto r^{-2}$), we notice the same temperature dependence in the two sets, while the heat conduction terms differ.

The collision terms for Coulomb collisions, using the Fokker–Planck collision approximation, have been derived in explicit form in the (usual) limit when relative flow differences are much smaller than thermal speeds ($(u_s-u_t)^2 \ll kT_{s(t)}/m_{s(t)}$), and by neglecting the contribution from the higher order term b ₂ (Cuperman et al. 1981). The expressions are sufficiently complex that they will not be reproduced here, but they are still in a form that should be straightforward to include in a numerical model.

This equation set accounts for heat conduction and thermal forces in the collision-dominated regime, and allows for a transition to collisionless flow. It may thus be an alternative to the equation set of Sect. 5.4.3 for a solar wind model that extends all the way from the solar transition region to interplanetary space. A disadvantage is that the first order correction b ₁ is still first order in c _s , and it is therefore not clear how well this set describes heat conduction and thermal forces in the transition region. On the other hand it has the advantage (compared to the set in Sect. 5.4.3) that a higher-order correction (ξ _s ) is included, which can accommodate a high-energy tail in the VDF. A comparison of the two sets would also be of interest; since the closure assumptions are not the same, it would be useful to know, e.g., how the proton heat fluxes differ in the two sets in the transition to collisionless flow in the outer corona.

To our knowledge, the full equation set (81), (111)–(114) has so far not been implemented and solved in any numerical model.

5.4.5 Gyrotropic Models Based on Observed Proton VDFs

So far all models have been expansions about Maxwellians or bi-Maxwellians, with the underlying assumption that collisions are sufficiently frequent that the departure from a Maxwellian will not be very large. However, far from the Sun collisions are so few that the VDF may depart significantly from a Maxwellian; indeed, the in situ solar wind proton observations confirm this (Marsch et al. 1982). Expanding about a base function that is far from the observed VDF requires that we may have to include a large number of terms in the correction term ϕ _s or χ _s to be able to reconstruct the observed VDF. We have already argued that increasing the order of the correction term will drastically increase the complexity of the equations. Moreover, with the forms for ϕ _s and χ _s discussed above, the VDF may even become negative in parts of velocity space when, e.g., the heat flux is large.

Leblanc and Hubert (1997, 1998) and Leblanc et al. (2000) have therefore proposed to use a base function that can reproduce the basic properties of some of the proton VDFs observed by the Helios spacecraft (Marsch et al. 1982), and particularly the high-energy tail of the VDF. The full VDF is assumed to have the shape

$$ f_s = f_s^{G} (1 + \psi_s) $$

(115)

where

$$ f_s^G = n_s \frac{m_s}{4\pi k T_{s\perp} D_s^*} \exp\left( - \frac{m_s}{2kT_{s\perp}} c_{s\perp}^2 - \frac{c_{s\parallel}+D_s^*}{D_s^*} + \frac{1}{E_s^*} \right) {erfc} \left[ \sqrt{E_s^*} \left( \frac{1}{E_s^*} - \frac{c_{s\parallel} + D_s^*}{2D_s^*}\right) \right], $$

(116)

and erfc is the complementary error function, with

$$ D_s^* = \sqrt[3]{\frac{q_s^\parallel}{2 m_s n_s}} $$

(117)

$$ E_s^* = \frac{2 m_s D_s^{*2}}{kT_{s\parallel} - m_s D_s^{*2}}. $$

(118)

In the limit $D_s^* \rightarrow 0$, $f_s^G$ becomes identical to the bi-Maxwellian $f_s^{bM}$ (73).

The expansion ψ _s is then written in a form very similar to (86), except that the coefficients are somewhat different because the base function is different. For one-dimensional flow along the magnetic field—neglecting the stress tensor and assuming that all heat flow is along the magnetic field—ψ _s simplifies to

$$ \psi_s = - \beta_{s\parallel} \beta_{s\perp} \left( 1 - \frac{\beta_{s\perp}}{2} c_{s\perp}^2 \right) \frac{q_{s\parallel}^\perp}{m_s n_s}c_{s\parallel}, $$

(119)

where $\beta_{s\parallel}$ and $\beta_{s\perp}$ are defined in Sect. 5.4.2.

This expansion leads to an equation set very similar to the 16-moment set of Sect. 5.4.2, with the main difference in the collision terms. The main drawback seems to be that collision terms have not been presented in a closed form suitable for inclusion in a numerical model.

To our knowledge, solar wind solutions with this equation set have not been obtained, neither with nor without the collision terms. Without actual numerical solutions, it is difficult to quantify the improvement offered by this approximation.

6 Summary and Conclusions

In this paper we have summarized the main arguments for and limitations of the kinetic-exospheric and multi-fluid models of the solar wind. Our aim was on one hand to outline a brief history of more than 50 years of solar wind modeling, as well as to provide a comparative review of kinetic and multi-fluid approaches emphasizing their advantages and limitations as well as their common theoretical roots—the fundamental equations of plasma state.

On the one hand we discussed the purely mechanical, or hydrodynamic, approach that explains the main properties of the solar wind by its macroscopic evolution in terms of density, bulk velocity and temperature. This description was the first to predict successfully the main bulk properties of the solar wind (Parker 1958b). On the other hand we emphasized the arguments of the holders of a more detailed description, the kinetic theory, that provides a detailed view on how the energy is distributed among plasma particles, pitch angles, from the collisional region out to the collisionless supersonic stream. This approach is built on the fundamental concept of the velocity distribution function, of each of the species streaming in the solar wind. It would be false, however, to state that the kinetic theory disregards the macroscopic plasma transport. Indeed, the entire hierarchy of the transport equations is satisfied by the kinetic solutions. Although the first generation kinetic exospheric models were less successfull to obtain the supersonic solar wind expansion, the corrections found for the Pannekoek-Rossland electric field (Lemaire and Scherer 1969) eventually provided the second generation kinetic exospheric models in agreement with observations and with the results of the hydrodynamic approach (Lemaire and Scherer 1971). At the same epoch Jockers (1970) developed independently a kinetic exospheric model for the supersonic solar wind. The third generation kinetic exospheric models by Lamy et al. (2003a) and Zouganelis et al. (2004) treats the effects of a non-Maxwellian electron velocity distribution function in the corona and gives clues about the kinetic processes in the transonic region. The fourth generation kinetic exospheric models by Pierrard et al. (1999, 2001a) includes Coulomb collisions and their effects on the velocity distribution functions and macroscopic parameters of the solar wind. The kinetic treatment of wave-particle interactions (Marsch 2006; Aschwanden 2009) is not discussed in this review, but gives insight on the anisotropies observed in the solar wind. After several decades of fluid and kinetic solar wind modeling, often leading to controversies, scientists now realize that both approaches are complementary and not opposed (Lemaire and Echim 2008; Parker 2010).

Fluid and kinetic models stem from the same theoretical root: the Boltzmann equation. The Chapman–Enskog and the Grad theories (see Appendix 2) are examples of solutions of the Boltzmann equation bridging the gap between microscopic and macroscopic descriptions. The Chapman–Enskog and the Grad solutions are based on the assumption that the zero-order approximation of the velocity distribution function of the protons and electrons, f(v, r, t), is necessarily a displaced Maxwellian, $f_M^{(0)}(v, r, t)$, characterized by a given number density n(r, t), a given average velocity u(r, t), and a given isotropic temperature T(r, t); the values of these lowest order moments of the actual VDF do not change when higher order approximations are defined and determined for $f(v, r, t) = f_M(v, r, t) \left[1 + g_1(v, r, t) + \ldots\right]$. In other words neither in the Chapman–Enskog approximation nor in Grad’s one the series of functions g _k contribute to changing the value of n, u and T; this means that the higher order terms f _M (v, r, t)g ₁(v, r, t) in the expansion of f, affect only the higher order moments, e.g. the non-diagonal components of the pressure tensor p _ij , the components of the energy flux tensor, etc.

Another key aspect of kinetic and fluid modeling is the treatment of the heat flux. In the kinetic treatment the heat flux is correctly computed from the moments of the velocity distribution function. By having a separate equation for the time evolution of the heat flux vector, thus treating the heat flux on an equal footing to density, flow speed, and temperature, the multi-fluid isotropic equations sets in Sects. 5.3.2–5.3.4, and the gyrotropic multi-fluid equations in Sects. 5.4.2 and 5.4.3, allow for a transition from classical heat conduction proportional to −∇T to a “collisionless” heat flux that will not be proportional to the temperature gradient. This is particularly important in multi-fluid solar wind models driven mainly by proton heating (as the recent UVCS/SOHO observations indicate), since protons will quickly become collisionless due to the strong heating. Whether the proton heat flux is classical or non-classical (collisionless) can have a large impact on the energy budget of the corona/solar wind system. If classical heat conduction dominates in a large region of the corona, much of the energy received by the protons may be conducted downwards and subsequently, through collisions, converted into electron heat conduction which is lost as radiation in the transition region. On the other hand, if the protons quickly become collisionless, more of the heating may be used to accelerate the wind. In solar wind models based on the 16-moment approximation (Lie-Svendsen et al. 2002) it was found that the protons would become almost collisionless quite low in the corona, so that the proton heat flux was positive (pointing outwards from the Sun) even below the proton temperature maximum, leading to a very high speed wind with a very low mass flux.

But the critical question remains. How well do these higher-order fluid models describe the proton heat flux in the transition to collisionless flow? Since the heat flux is a higher-order moment, one might expect that it is sensitive to the shape assumed for the VDF (the closure assumption). At present we have no independent way of evaluating the accuracy of the proton heat flux description of these equations in this transition regime. This is one area in which kinetic models could provide valuable insight. By actually solving the Boltzmann equation for protons in the outer corona, with proton-electron and proton-proton collisions included through the Fokker–Planck collision term, one could evaluate how well the higher-order fluid models do the energy bookkeeping in this region. This is a challenging task for the kinetic models, however: since collisions are still essential in this region, collisionless kinetic models (exospheric models) are inadequate. Moreover, proton self-collisions are needed. Since the proton VDF is expected to deviate strongly from a Maxwellian, in particular if the heating is caused by cyclotron waves heating mostly perpendicular to the magnetic field, the test particle approach developed by Lie-Svendsen et al. (1997) and Pierrard et al. (1999) for electron self-collisions is most likely not directly applicable to protons. Despite the challenges, or because of them, this is a problem that definitively deserves attention in order to improve our understanding of the energy budget of the corona and solar wind.

Our review emphasizes that the kinetic approach is fundamental and provides a rigorous treatment of the physics at all levels, going beyond the limits of an academic exercise. Indeed, kinetic models describe microscopic and macroscopic processes, addressing self-consistently the complex aspects of the interactions between plasma particles and fields. Such a detailed level of description is obtained at the expense of a lower spatial resolution and is currently limited to stationary situations. Future development of new kinetic Vlasov and Fokker–Planck simulations are expected to give new momentum to kinetic modeling of the solar wind.

Appendices

Appendix 1

The following list defines the first moments of the velocity distribution function of each species s, and the main macroscopic properties of the plasma itself. A similar list can be found in Montgomery and Tidman (1964; p. 197), the only difference is that we define the peculiar, or thermal, velocity c _s  = v − u _s , following Grad (1958) and Schunk (1977), i.e. in the reference frame moving with the average velocity, u _s , of each species:

$$ n_s = \int\int\int{f_s\, d {\mathbf{v}}}, \hbox{partial number density of species}\, s $$

(120)

$$ n = \sum_{s}n_s, \hbox{total number density of the plasma} $$

(121)

$$ {\mathbf{u}}_{s} = \frac{1}{n_s}\int\int\int{f_s\, {\mathbf{v}}\, d {\mathbf{v}}}, \hbox{average velocity of species}\, s $$

(122)

$$ \rho_m = \sum_{s}{n_s}m_s, \hbox{total mass density of the plasma} $$

(123)

$$ {\mathbf{U}} = \frac{1}{\rho_c}\sum_{s}{m_s n_s {\mathbf{u}}_s}, \hbox{center of mass velocity or bulk velocity of the plasma} $$

(124)

$$ {\mathbf{w}} = {\mathbf{v}} - {\mathbf{U}}, \hbox{particle velocity with respect to the center of mass velocity} $$

(125)

$$ {\mathbf{c}}_{s} = {\mathbf{v}} - {\mathbf{u}}_{s}, \hbox{random velocity} $$

(126)

$$ {\mathbf{U}}_{s} = \frac{1}{n_s}\int\int\int{f_s {\mathbf{w}}\, d {\mathbf{v}}}, \hbox{average velocity of species } s \hbox{with respect to center of mass velocity} $$

(127)

$$ {\mathbf{J}} = \sum_s{e_s n_s {\mathbf{u}}_s}, \hbox{total electric current density of the plasma} $$

(128)

$$ \rho_c = \sum_s{e_s n_s}, \hbox{total electric charge density of the plasma} $$

(129)

$$ \underline{\underline{p_s}} = m_s \int\int\int{{\mathbf{c}}_s {\mathbf{c}}_s\, f_s\, d{\mathbf{v}}},\;\hbox{kinetic pressure tensor of species } s $$

(130)

$$ \underline{\underline{P}} = \sum_s \underline{\underline{p}}_s, \hbox{total pressure of the plasma} $$

(131)

$$ \underline{\underline{\tau}} = \sum_s m_s \int\int\int {\left({\mathbf{c}}_{s} {\mathbf{c}}_{s} - 1/3c_s^2 \underline{\underline{1}} \right)f_s\, d{\mathbf{v}}}, \hbox{stress tensor, traceless part of}\, \underline{\underline{P}} $$

(132)

$$ \underline{\underline{Q_s}} = m_s \int\int\int{{\mathbf{c}}_{s} {\mathbf{c}}_{s} {\mathbf{c}}_{s}\, f_s\, d{\mathbf{v}}}, \,\hbox{heat tensor of species } s $$

(133)

$$ {\mathbf{q}}_{s} = \frac{m_s}{2} \int\int\int{c_s^2\, {\mathbf{c}}_{s}\, f_s\, d{\mathbf{v}}} \,\hbox{heat flux vector of species } s $$

(134)

$$ {\mathbf{q}} = \sum_{s} \,{\mathbf{q}}_{s}\, \hbox{total heat flux vector of the plasma} $$

(135)

$$ T_{s\parallel} = \frac{m_s}{kn_s} \int\int\int{(c_s)_{\parallel}^2\,f_s\, d{\mathbf{v}}}\, \hbox{the partial parallel temperature of species} s $$

(136)

$$ T_{s\bot} = \frac{m_s}{2kn_s} \int\int\int{(c_s)_{\bot}^2\,f_s\, d{\mathbf{v}}} \,\hbox{the partial perpendicular temperature of species}\;s $$

(137)

$$ T_s = \frac{1}{3}(T_{s\parallel} + 2T_{s\bot}) \,\hbox{the omnidirectional temperature of species} s $$

(138)

$$ T_{\parallel} = \frac{1}{n} \sum_s n_s T_{s\parallel} \,\hbox{the plasma parallel temperature} $$

(139)

$$ T_{\bot} = \frac{1}{n} \sum_s n_s T_{s\bot} \,\hbox{the plasma perpendicular temperature} $$

(140)

$$ T = \frac{1}{3}(T_{\parallel} + 2T_{\bot}) \,\hbox{the plasma temperature} $$

(141)

Appendix 2

In this appendix we briefly review two methods adopted to solve the Boltzmann equation, and which relevant for kinetic and fluid modeling of the solar wind. A historical review of the approaches adopted to solve the Boltzmann equation may be found in the introduction of Chapter 8 of the monograph by Jancel and Kahan (1966).

In order to illustrate the process by which one can reduce the solutions of the Boltzmann equation to fluid solutions, and thus switch from the formalism based on the velocity distribution functions to the one based on their first moments, it is useful to introduce a smallness parameter, ε:

$$ \frac{{\fancyscript{D}}f_s}{{\fancyscript{D}}t}= \frac{1} {\varepsilon}J_s $$

(142)

where ${\fancyscript{D}}f_s/{\fancyscript{D}}t$ is identical to the left-hand side of (2) and J _s is the collision integral defined in (5). ε is associated with the smallest spatial scale (can be the mean-free-path) or the smallest time period (the mean collision time) or both. In this Appendix we shall also denote the left hand side of the equation (142) by $\dot{f}$.

A special class of the Boltzmann equation is the so-called “normal solutions” that are valid only when the mean free path, λ _c , and the collision time, τ _c are small with respect to the spatial and temporal scale, L or ${\fancyscript{T}}$ respectively: i.e. ε≪1. The normal solutions are obtained by expanding f _s as a power series of $\varepsilon = \lambda_c/L$ or $\varepsilon = \tau_c/{\fancyscript{T}}$. The class of normal solutions is built on the rigorous mathematical proof given by Hilbert (1924) for the existence of a solution for a Fredholm type equation like the Boltzmann equation. The expansions proposed first by Chapman (1918) and Enskog (1917) give normal solutions of interest for physicists.

1.1 Chapman–Enskog Expansion

Enskog (1917) and Chapman (1918) developed a method to solve the Boltzmann equation suited for practical situations that require a smaller number of initial conditions, as in fluid dynamics approaches. In addition to the expansion

$$ f_s = f_s^{(0)} + \varepsilon f_s^{(1)} + \varepsilon^2 f_s^{(2)} \ldots $$

(143)

Enskog adds also an expansion of the temporal derivative operator:

$$ \frac{\partial }{\partial {t}} \equiv \sum_{r=0}^\infty\varepsilon^r\frac{\partial {^{(r)}}}{\partial {t}} $$

(144)

and also defines an expansion for the moments $\overline{M}$, based on (143) as given below:

$$ \overline{M} = \overline{M}^{(0)} + \varepsilon \overline{M}^{(1)} + \varepsilon^2 \overline{M}^{(2)} +\ldots = \frac{1}{n}\int \int \int M{\left[ f_s^{(0)} + \varepsilon f_s^{(1)} + \varepsilon^2 f_s^{(2)} \ldots \right]} d{\bf c}_s $$

(145)

When the expansion (145) is inserted into the conservation equations [continuity, momentum and energy, see (20)–(22)] the corresponding results define the formal expansion of the time derivative operator (144) as a function of spatial derivatives only:

$$ \frac{\partial {^{(0)} n^{(0)}}}{\partial {t}} = -\varvec{\nabla}_r \cdot \left(n^{(0)} {{\bf u}^{(0)}} \right), $$

(146)

$$ \frac{\partial {^{(0)} {{\bf u}^{(0)}}}}{\partial {t}} = -\left({{{\bf u}^{(0)}}}\cdot \varvec{\nabla}_r \right) {{\bf u}^{(0)}} + \frac{{\bf F}}{m} -\frac{1} {\rho}\varvec{\nabla}_rp^{(0)}, $$

(147)

$$ \frac{\partial {^{(n)}{\bf u}^{(n)}}}{\partial {t}} = - \frac{1}{\rho} \varvec{\nabla}_r\cdot \underline{\underline{{\overline{p}^{(n)}}}},\quad {n>0} $$

(148)

$$ \frac{\partial {^{(0)}T^{(0)}}}{\partial {t}} =-{\bf u}^{(0)} \cdot \varvec{\nabla}_r T^{(0)} - \frac{2T^{(0)}}{3}\varvec{\nabla}_r \cdot {\bf u}^{(0)}, $$

(149)

$$\frac{\partial {^{(n)}T^{(n)}}}{\partial {t}} = - \frac{2}{3kn^{(0)}}\left[\underline{\underline{{\overline{p}^{(n)}}}}:\varvec{\nabla}_r{\bf u^{(0)}} + \varvec{\nabla}_r\cdot \overline{{\bf q}}^{(n)}\right], {n>0} $$

(150)

The above equations are specific for each species but in order to keep the notation simple the index s corresponding to the species has been dropped. $p=\overline{p^{(0)}}$ is the isotropic pressure corresponding to the isotropic Maxwellian $f^{(0)}$, while $\underline{\underline{{{p}^{(n)}}}}$ and q ⁽ⁿ⁾ are defined respectively by (145) for M = c c and M = c ² c (see also Appendix 1 for the definitions of $\underline{\underline{{p}}}$ and q).

From the so-called integrability conditions (Hilbert 1924) one obtains a set of equations equivalent to:

$$ \overline{n}^{(r)}= 0, {r\ge1} $$

(151)

$$ \overline{{\bf v}}^{(r)}= 0, {r\ge1} $$

(152)

$$ \overline{T}^{(r)}= 0, {r\ge1} $$

(153)

The Eqs. (151)–(153) state that at any time the higher order corrections of f do not modify the density, average velocity and temperature. The latter are fully determined by f ⁽⁰⁾: the zero order approximation of f which is postulated to be an isotropic Maxwellian distribution function. The Eqs. (151)–(153) are restricting arbitrarily the class of normal solutions of the Boltzmann equation. In certain situations, e.g. when the Knudsen number is of the order of unity or larger, the restrictive conditions may not be appropriate and justified. This appears to be the case in the collisionless region of the interplanetary medium where the VDF departs significantly from displaced Maxwellians.

The Chapman–Enskog solutions of the Boltzmann equation are necessarily based on the assumption that ε is small and on the restrictive conditions (151)–(153). As put forward by Grad (1958) an essential consequence of the assumptions of Chapman and Enskog expansion is the existence of two well separated time scales: (1) a (small) collision time, i.e. the time necessary for the non-Maxwellian flow to relax and (2) a macroscopic decaying time where dissipation can be described in terms of traditional thermodynamic processes/variables like diffusion, viscosity, and heat conductivity. When ε is not small, i.e. when the mean-free-path or the collision frequency are no longer small compared to the dimensions of the system, the results obtained based on the expansion (143) and on the constraints (151)–(153) are no longer applicable.

In the Chapman–Enskog kinetic theory outlined above, displaced Maxwellian VDFs are adopted as the zero order approximation of f(v, r, t). Nevertheless, when K _n  > 1, the VDFs of collisionless particles are not well approximated by displaced Maxwellian functions, which correspond to equilibrium VDFs in a uniform collision dominated gas or plasma. For instance in a Knudsen gas where K _n  > 1, the kinetic pressure tensor (i.e. the second order moments of the VDF) departs from its classical Navier-Stokes approximation; furthermore, the expression for the heat flux departs then from the classical Fourier law: q =  −λ∇T, where λ is the thermal coefficient of the gas.

1.2 Expansion About a Local Maxwellian: Grad’s Method

When the smallness factor ε used to expand the distribution function, f(v, r, t) is no longer much smaller than 1 the normal solutions are not valid. In Grad’s method the expansion of the VDF f is again based on a local isotropic Maxwellian in the frame of reference moving with the bulk velocity u(r, t):

$$ f^{(0)}=n\left( {\bf r}, t \right) \left(\frac{m}{2\pi kT({\bf r}, t)} \right)^{\frac{3}{2}} \exp\left[{-\frac{m\left( c({\bf r}, t) \right)^2}{2kT({\bf r},t)}}\right] $$

(154)

The density, n(r, t), temperature, T(r, t), as well as u(r, t) are not constants but functions of time and spatial coordinates. This may be justified when there is a clear separation between spatio-temporal scales specific to collisions and the hydrodynamic evolution of the fluid (Grad 1958). This implies that the mean-free path of particles is much smaller than the scale length of the gas, and that the characteristic collision time is much smaller than the time over which u, n and T change significantly. When this is no longer the case the expansion about a local Maxwellian is not applicable; an expansion about an absolute Maxwellian (i.e. with n, U and T independent of r) might be then in order, as pointed out by Grad (1958).

The method of Grad (1958) consists in expanding the distribution function using the orthogonal base of the Hermite polynomials, with the local Maxwellian, f ⁽⁰⁾, as a weight function:

$$ f=f^{(0)}\sum_{n=0}^{\infty}\frac{1} {n!}a_{i_n}^{(n)}{\fancyscript{H}}_{i_n}^{(n)} $$

(155)

where the term $a_{i_n}^{(n)}{\fancyscript{H}}_{i_n}^{(n)}$ stands for the sum

$$ \sum_{i_1i_2\ldots i_n}a_{i_1i_2\ldots i_n}^{(n)}{\fancyscript{H}}_{i_1i_2\ldots i_n}^{(n)}. $$

${\fancyscript{H}}_{i_n}^{(n)}$ denotes a tensor with n subscripts i _n  = i ₁ i ₂...i _n whose components are polynomials of n-th degree as defined below:

$$ \begin{aligned} {\fancyscript{H}}^{(0)} &= 1,\quad {\fancyscript{H}}_i^{(1)} = c_i, \quad {\fancyscript{H}}_{ij}^{(2)} = c_ic_j - \delta_{ij}, \\ {\fancyscript{H}}_{ijk}^{(3)} &= c_ic_jc_kc_l - \left(c_i\delta_{jk} + c_j\delta_{ik} + c_k\delta_{ij} \right),\,\ldots \end{aligned} $$

The coefficients $a_{i_n}^{(n)}$ are tensors determined by an inversion formula:

$$ a_{i_n}^{(n)}= \frac{1}{n}\int \int \int{f{\fancyscript{H}^{(n)}_{i}}} d{\bf v} $$

(156)

giving a sequence of moments of increasing order:

$$ a^{(0)} = 1,\quad a_i^{(1)} = 0,\quad a_{ij}^{(2)} = \tau_{ij}/p,\quad a_{ijk}^{(3)} = Q_{ijk}/{p\sqrt{kT}} \ldots $$

(157)

where the definitions of moments from Appendix 1 have been used.

The formulas found for a ⁽⁰⁾, $a_i^{(1)}$ and $a_{ij}^{(2)}$ state that the density, n(r, t), the mean velocity, u(r, t), and the temperature, T(r, t), of the distribution function f are uniquely determined by the corresponding moments of the zero order approximation, f ⁽⁰⁾, the Maxwellian distribution function. Consequently $a_{i_n}^{(n)}$ vanish for n > 2 when the VDF is strictly an isotropic and Maxwellian VDF, i.e. when f = f ⁽⁰⁾. In Grad’s class of kinetic solution defined by (155) the higher order coefficients (n > 2) do not contribute to the five first order order moments: n(r, t), u(r, t), and T(r, t). This imposes a serious restriction and limitation to Grad’s method. This restriction is similar to that adopted in the Chapman–Enskog method and described by (143) and (151)–(153).

Grad’s expansion (155) is inserted into the Boltzmann equation (2), and after equating the coefficients of $H_{i_n}^{(n)}$ one obtains an infinite series of differential equations for the linear combinations of the moments of f. This unlimited series of partial differential equations defines the hierarchy of moments equations which leads to an alternative approximation of the hydrodynamic transport equations. We show here the equations corresponding to the second order terms, n = 2:

$$ \begin{aligned} \frac{\partial {a_{ij}^{(2)}}}{\partial {t}}& + u_r \frac{\partial {a_{ij}^{(2)}}}{\partial {x_r}}+ a_{ir}^{(2)} \frac{\partial {u_{j}}}{\partial {x_r}} + a_{rj}^{(2)} \frac{\partial {u_{i}}}{\partial {x_r}}+ \left( a_{ij}^{(2)} +\delta_{ij} \right)\frac{1}{kT}\frac{D(kT)}{Dt}+ \sqrt{kT}\frac{\partial{a_{ijr}^{(3)}}}{\partial {x_r}}\\ &\quad +\frac{\sqrt{kT}}{n} a_{ijr}^{(3)} \frac{\partial{n}}{\partial {x_r}} + \frac{3}{2kT} a_{ijr}^{(3)} \frac{\partial {kT}}{\partial {x_r}} + \frac{\partial {u_i}}{\partial {x_j}}+ \frac{\partial {u_j}}{\partial {x_i}}= J_{ij}^{(2)} \end{aligned} $$

(158)

where tensor notation and summation convention are used. Remember that $D/Dt\equiv(\partial/\partial t)+u_r(\partial/\partial x_r)$ is the convective time derivative. When the force between particles is specified, the contribution of the collision term, denoted $J_{ij}^{(2)}$ in (158), is given by an infinite quadratic form of all the Hermite coefficients (see Grad 1958, Eqs. 28.14–28.17).

An important property of the Grad’s expansion is evidenced by (158): the differential equation for the moment a ⁽ⁿ⁾ involves spatial derivatives of the higher order moment a ⁽ⁿ⁺¹⁾. Nevertheless, the corresponding collision term, $J_{ij}^{(2)}$, is given by sums involving the infinite sequence of Hermite coefficients of all orders, n = 0,…, ∞. In order to obtain tractable solutions for the collision term, the expansion (155) is truncated and replaced by a finite sum of the first n Hermite polynomials. Thus the Grad solution of Boltzmann equation assumes that all the moments, or Hermite coefficients, of order (n + 1) and higher are set equal to zero:

$$ a_{i_{n+1}}^{(n+1)} = 0 $$

(159)

This ad hoc truncation is not based on physics principles other than convenience to limit the mathematical complexity of the system of differential equations to be solved.

The left hand side of only one equation of Grad’s moment hierarchy is affected by the truncation (159), namely the equation corresponding to the highest order considered, n. Nevertheless, the collision integral on the right hand side of all equations of the hierarchy is altered by the truncation (159). A total number of M differential equations results by truncating the Grad’s hierarchy at the (Hermitian) order (n), therefore the corresponding Grad expansion is also called the M-moments method. In Sect. 5 we discussed the application of the 5-, 8-, and 16-moment equations to modeling the solar wind.

When the Grad expansion is limited to a few Hermite coefficients the solution can be interpreted as a generalization of the usual multi-fluid hydrodynamics. When the expansion is more elaborate and more coefficients (or moments) are considered then Grad’s method provides a microscopic estimation of the distribution function described by the subclass of VDF determined by (155).

The convergence of the Chapman–Enskog and Grad expansions of the velocity distribution functions (143) and (155) has never been proven mathematically, not even in the case of gases for which ε is much smaller than unity, i.e. for Knudsen numbers much smaller than unity. A fortiori, the convergence of these expansions is not even valid in collisionless gases and plasmas, as for instance in coronal and stellar ion-exospheres. Therefore, it would be presumptuous to consider that the higher the order of the Chapman–Enskog or Grad expansion is, the better will be the solutions of the system of fluid equations used to model the solar wind at large heliospheric distances as well as any other collisionless astrophysical plasmas.

A relevant example may be found in the field of classical hydrodynamics. In the 1950s and 1960s the propagation and attenuation of sound waves was investigated in laboratory experiments. The wave number and the wave attenuation constant, at a given frequency, were measured experimentally (Greenspan 1956; Toba 1968). These propagation and attenuation constants were also determined theoretically using different hydrodynamical approximations of the transport equations and of the moment equations (Truesdell 1955; Sirovich and Thurber 1965). All these studies confirm that the theoretical predictions of these wave propagation constants (phase velocity or wave number, and attenuation constant) fail to fit the measured values, when the ratio between mean free path of the particles and the wavelength becomes of the order and larger than unity: i.e. when the Knudsen number of the gas is of the order of 1 or larger.