A&A 400, 769-778 (2003)
DOI: 10.1051/0004-6361:20021822
S. Landi1 - F. Pantellini2
1 - Dipartimento di Astronomia e Scienza dello Spazio,
Largo Enrico Fermi 2, 50125 Firenze, Italy
2 -
Observatoire de Paris, LESIA,
5 Place Jules Janssen, 92195 Meudon, France
Received 30 August 2001 / Accepted 9 December 2002
Abstract
We present and discuss a completely self-consistent kinetic simulation of
a steady state transonic solar type wind. The equations of
motion of an equal number of protons and electrons plunged
in a central gravitational field and a self-consistent
electric field are integrated numerically.
Particles are allowed to make binary collisions with
a Coulombian scattering cross-section. The velocity
distributions of the particles injected
at the boundaries of the simulation domain
are taken to be Maxwellian.
As anticipated by previous authors we
find that the transonic solution implies the existence
of a peak in the proton equivalent potential at some
distance above the sonic critical point. Collisions appear to be
the fundamental ingredient in the process of accelerating the
wind to supersonic velocities. For a given temperature at the
base of the simulation domain the acceleration efficiency
decreases with decreasing density. The reason is that
the plasma has to be sufficiently collisional for the
heat flux to be converted efficiently into plasma bulk velocity.
Concerning the heat flux we find that
even when in the vicinity of the sonic point
the collisional mean free path of a thermal
particle is significantly smaller than
the typical scales of variation of
the density or the temperature, the electron heat flux cannot be described
conveniently by the classical Spitzer-Härm conduction law;
not even in most of the subsonic region. Indeed, in the simulations where
a transonic wind forms the heat flux has been found to
strongly exceed the Spitzer-Härm flux, in opposition
to recently published results from multi-moment models.
We emphasize that given the high coronal temperatures we use in our simulations
(3 times the typical solar values) we do not expect the results
presented in this report to be uncritically transposable to the case of
the "real'' solar wind. In particular, the quantitative aspects of our results
must be handled with some care.
Key words: Sun: solar wind - stars: winds, outflows - plasmas - conduction - methods: numerical
At all heights, from the bottom of the corona up into the interplanetary space, the solar atmosphere is a permanently expanding, out of thermodynamic equilibrium and fully ionized plasma. During the 1950's the recognition of the weak collisionality of the solar wind conveyed some doubts concerning the ability of fluid models to describe the solar atmosphere conveniently. As a consequence Chamberlain (1960), largely influenced by the theories on gas evaporation from planetary atmospheres, published the first kinetic model of the solar wind. In Chamberlain's evaporation model the wind is subsonic at the Earth's orbit in clear opposition with the supersonic solution of the fluid equations proposed a few years earlier by Parker (1958). During the early 1960's in situ measurements confirmed the supersonic nature of the solar wind and kinetic models just fell in disuse for some time. In the early 1970's it became clear that Chamberlain's erroneous prediction of a subsonic wind was the consequence of having mistakenly assumed that the charge neutralizing electric field was the Pannekoek-Rosseland field (e.g., Rosseland 1924). The latter is based on the assumption of a static solar atmosphere which has been a privileged working hypothesis since Laplace's Traité de mécanique céleste, published in the early years of the nineteenth century, but has been shown to be completely at odds with observations. After the definitive relaxation of the static approximation for the electric field, kinetic models of the solar wind were back on stage again (Jockers 1970).
The simplest, and most widely used kinetic models
are the so called exospheric models
(e.g., Lemaire & Sherer 1971). In these models
the solar atmosphere is assumed to change
from fully collisional to collisionless
at a sharply defined level called the exobase.
Above the exobase, conventional exospheric
models assume a monotonically decreasing
equivalent proton potential
(cf. Eq. (4))
and no protons coming into the system from infinity, so that,
by construction, all protons have
anti sunward directed velocities.
For non pathological distribution functions this means that
the plasma bulk velocity at the exobase is of the order
of the proton thermal velocity, i.e. of the order of the
sound velocity. For example,
is the mean velocity of a proton population with Maxwellian
velocity distribution
truncated for
(subscripts
and
refer to the radial
direction with respect to the center of the Sun).
Since the typical velocity of a proton at the exobase is
by construction of the order of the radial bulk velocity, it
follows that the exobase is
located at a heliocentric distance comparable
to the distance r* where the subsonic-supersonic transition
is located (the sonic critical point).
However, as demonstrated graphically by Jockers
(1970),
the proton potential cannot be monotonic from deep inside the corona, where
the bulk velocity is supposed to be small compared to sound speed and where
static approximation may apply, out to infinity,
where the wind is supersonic. Jockers anticipated that
on its way from the corona to infinity
a transonic wind must overcome a
maximum in the proton potential
such that
,
where
is the location of the maximum.
In two recent papers Scudder (Scudder 1996a,b)
pushes a step farther by identifying
the critical point of Parker's fluid model
with the location of the maximum of the proton
potential energy
.
Based on that assumption
he derives a number of
constraints on the possible radial variations of both
the proton and the electron temperatures near the sonic point.
However, even though the existence of a transonic
wind seems to be intimately related to the existence
of a peak of the potential
,
there
is no reason for
to coincide with
the sonic point of fluid theories unless very special, and therefore
unlikely, conditions are met there.
For example, in the simulations presented in this
paper we do always find
.
It can be shown analytically that this is indeed
the normal case for radially decreasing temperature profiles
provided T decreases more slowly than r-1 (Meyer-Vernet et al. 2002).
In this paper we present self-consistent kinetic simulations of a stationary solar type wind, where we concentrate on those aspects which cannot be addressed by fluid theories such as the electric field, the heat flux and the collisionality of the plasma. We deliberately treat only the most simple case of Maxwellian boundary conditions for the particles' velocity distribution function. The effects of plasma instabilities and waves are also not included in the model. Such additional "complications'' may hide part of the fundamental physics of the acceleration process and shall be discussed elsewhere. In this respect we do merely mention that the effect of resonant waves on a hybrid (fluid + kinetic) solar wind model has been discussed by Tam & Chang (1999) who conclude that ions may well be accelerated more efficiently by resonant waves rather than by the radial electric field. A similar model has been used by Lie-Svendsen & Leer (2000) to show that the two temperature electron velocity distribution functions often observed in the solar wind can be generated by Coulomb collisions without the need of assuming the presence of non-Maxwellain distribution in the corona. Olsen & Leer (1999) and Li (1999) use a closed system of transport equations based on an anisotropic bi-Maxwellian approximation for the velocity distribution functions to simulate the solar wind from the lower corona outward. Lie-Svendsen et al. (2001) extend the model down to the chromosphere based on the argument that chromosphere, transition region, corona and solar wind constitute a coupled system. The system of equations used by these authors is known as the 16-moment approximation (e.g., Demars & Schunk 1979) is a fluid-type model including transport effects such as heat flow and viscosity and even Coulomb collisions between interacting species. The main purpose of the above authors was to reproduce as good as possible the characteristics of the coronal plasma by including ad-hoc heating functions supposed to mimic the local deposit of energy due to plasma waves. If suitably chosen the heating functions can reproduce the temperature profile and temperature anisotropies which observations suggest to prevail in the solar corona (e.g., Esser et al. 1999). In many respects, our model is much more limited than the above multi-moment models which include most of the ingredients (e.g. radiation, plasma heating through waves, collisions, etc.). However, these models are fundamentally fluid models and many ingredients are not self-consistent. Our model is kinetic and fully self-consistent, but neither waves nor radiation and not even the lower layers of the corona are taken into account. In addition, because of computational limitations, we use an artificially low proton to electron mass ratio, and an exceedingly high coronal temperature, so that transposition of our results to the case of the real Sun must be done critically. One substantial difference between our results and the results from the mentioned multi-moment models is that we find that the electron heat flux in the corona is one order of magnitude larger than the classical value predicted by the Spitzer-Härm formula (Spitzer & Härm 1953), whereas the multi-moment models find it to be of the same order. Of course, both models are subject to their own limitations so that the question of whether the heat flux in the solar wind is classical or not appears to remain an open question.
Even though the physical parameters characterizing the wind simulated in the following section do not correspond exactly to those observed for the Sun, we shall refer to the simulated wind as the solar wind and to the central star as the Sun.
Details of the simulation model have been given in two previous papers
(Pantellini 2000; Landi & Pantellini 2001)
and shall not be repeated here to full extent.
The model is spatially one dimensional, i.e. all fields
depend on the heliocentric distance r only.
An equal number of protons and
electrons are allowed to move freely in the domain
,
where r0 is the solar
radius and
is the outer boundary of the system
located several solar radii beyond the sonic
point. The equations of motion are those of
a particle of mass m and charge q in a central gravitational field
produced by a star of mass M and a radial,
charge neutralizing electric field,
,
i.e.
The physical state of the solar corona at heliocentric distance
r0 (the lower boundary of our simulation domain)
is characterized by the dimensionless parameter
defined as
(
is the Boltzmann constant)
![]() |
(3) |
The equations of motion Eqs. (1) and (2) are
integrated for N protons, and an N electrons
in the radial distance range between
r=r0 and
.
The number N
is determined by the requirement of the collision frequency
of an electron in the system (near r=r0) being roughly equal to the
Fokker-Planck collision frequency of a
plasma with an electron number density
,
which is a typical figure in the solar corona.
The so calculated number N turns out to be of the order 103
(Landi & Pantellini 2001).
Each time a proton or an electron hits the boundary at r=r0
it is injected back into the system according to a non
drifting isotropic Maxwellian velocity distributions
with a temperature T0.
Given that the protons reaching the top boundary at
are generally either supersonic or nearly supersonic most of
it must be re-injected into the system at r=r0.
On the other hand, electrons reaching the
boundary are
injected back into the system either at r=r0 or at
depending on what is needed to make the electron flux
to be equal to the proton flux (zero charge current condition).
The injection method ensures that there are always
N protons and N electrons in the system.
The velocity distribution of the electrons injected at the
top boundary is chosen to be a drifting bi-Maxwellian
with radial and perpendicular
temperatures equal to the radial and perpendicular temperature of the
outgoing electrons. Finally, the drift velocity of the electrons
at
is taken to be equal to the drift velocity of the protons measured
at
.
In case of a supersonic wind this is just
the average velocity of the protons escaping from the top of the system.
In a subsonic wind some protons have to be re-injected into the
system from the top. The temperature and the
number of the re-injected protons is then adjusted iteratively until a coherent
solution is obtained, in a manner similar to the one used
by Landi & Pantellini (2001) to simulate the
static corona.
The electric field profile is adjusted iteratively, during an
initialization phase, until zero charge flux and local charge neutrality
is achieved in all points of the system.
The number density n and the collisionality
of the simulated plasma is dependent on the number
of particles N used in the model.
Figure 1 shows the results of 4 simulations which only differ
in the number N of simulated particles, i.e.
N=400, 784, 1600, 6400, the
N=400 run being the one with the most tenuous (i.e. less collisional)
atmosphere. The corresponding number densities at the base of the
system are
and 13.4, respectively.
The differences between the 4 simulations are substantial
in many respects. The most evident difference is that the wind acceleration is
much more efficient in the high density case, even though the thermal
Knudsen number
(where
is the electron-proton collisional mean free path
based on the electron-proton rate of momentum exchange
(Landi & Pantellini 2001, Appendix B), and where
is the radial
thermal electron velocity)
is much smaller than unity for all runs, ranging from
10-3, for the densest case, to 10-2 for the most tenuous case.
From the figure it appears that the two more tenuous cases do not even
become supersonic with respect to radial proton thermal velocity
(which coincides with the fluid isothermal sound speed when
). The curves
on the bottom panel of Fig. 1
illustrate the effect of collisions on the proton potential energy
.
The latter results from the sum of the gravitational potential and a
charge neutralizing electrostatic potential
:
![]() |
Figure 1:
Flow velocity, Mach number and proton potential energy profiles
for four different values of the number of particles N in the system.
From bottom to top the velocity profiles correspond to
N=400, 784, 1600, 6400, respectively.
In all runs ![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
Open with DEXTER |
In summary, the main consequence of increasing the plasma
density beyond some threshold appears to be a reduction
of the potential barrier the protons
have to overcome in order to escape to infinity, accompanied by
the formation of a local maximum in the
profile.
As we shall see below the formation of the maximum in the
proton potential energy is intimately related to the
existence of both an outward directed, and radially
decreasing heat flux, and a radially decreasing temperature profile.
Both contribute in strengthening the outward directed electric field
,
thus facilitating the extraction of the protons.
In this context we shall remember that if the plasma was static
(impermeable boundaries) with equal electron and proton temperatures,
the charge neutralizing potential
would be the
celebrated Pannekoek-Rosseland potential (e.g., Rosseland 1924)
![]() |
Let us now address the question of the wind acceleration mechanism.
In order to do so we write the energy equation for the species j
for the case of a steady state and purely radial wind.
Indeed, the second moment of Boltzmann's equation (e.g., Endeve & Leer 2001)
leads to the following expression
![]() |
Figure 2:
Proton and electron energy profiles obtained
by evaluating Eq. (6)
for the N=6400 simulation. Note how both,
![]() ![]() ![]() |
Open with DEXTER |
![]() |
(11) |
![]() |
Figure 3: Relative importance of each term in Eq. (10) for the most dense case N=6400 (top panel) and the most tenuous case N=400 (lower panel). The dense case supports a transonic wind (the vertical line in the top panel gives the position of the sonic point r*) while the tenuous case does not. From the comparison of the two figures it appears that the wind acceleration is primarily driven by the heat flux term q/(nv). |
Open with DEXTER |
![]() |
Figure 4: Radial and perpendicular temperature profiles for 4 different values of the N. Note the log-log axis of the top two panels. |
Open with DEXTER |
Figure 4 shows that the radial electron temperature profile
is essentially insensitive to the density
while
is not. Indeed, in the collisionless limit the parallel and
perpendicular temperatures are independent of each other and one should observe
in order to satisfy to the angular
momentum conservation law of individual particles (cf. Eq. (2)).
As a consequence, in the
rigorously collisionless case,
should decrease by a factor
502 from bottom to top of the simulation domain.
Collisions do significantly contribute in limiting this bottom to top
perpendicular electron temperature gradient which ranges from 10 to 30
depending on the value of N.
On the other hand, for all four cases the parallel temperature
only drops by a factor 3 from bottom to top, leading to strong temperature
anisotropies
(bottom panel in Fig. 4). This is not particularly surprising as
in the collisionless limit the parallel temperature
of a plasma plunged in a potential field is constant as long as
the velocity distribution function is close to Maxwellian.
We can now summarize the wind acceleration scenario from a
kinetic point of view. The natural decrease of the temperature
with distance (essentially due to the
fact that in the collisionless limit
)
implies the existence of a radial heat flux
predominantly conveyed by the electrons.
The heat flux is transfered from the electrons
to the protons which become accelerated
in the outward direction.
Since the momentum of the
wind is mainly carried by the heavy protons (rather than by
the light electrons) the plasma as a whole becomes accelerated
in this way.
This mechanism must be particularly efficient in the
region located inside the spherical shell
(location of the maximum of
)
where the protons have to climb uphill in order to escape
from the potential energy well (cf. Fig. 1).
As already stated above, collisions contribute in increasing
the electric field strength. Since the electric field is directed
outward, increasing the electric field favors the extraction of the
protons from the gravitational well by reducing the height of the
maximum in the proton potential
.
The fluid estimate of the macroscopic electric field
for a
spherically symmetric electro-proton plasma can be obtained
by differentiating Eq. (6) for the electrons.
Neglecting the small terms proportional to
and taking advantage of the fact that
is approximately constant (in particular if one
compares it to the
profile shown
in Fig. 2) we then obtain
It is instructive to apply Eq. (12) to a
a weakly collisional system. In such a case the parallel temperature
is roughly constant and the perpendicular
temperature
.
This leads to
the collisionless approximation
![]() |
Figure 5:
Parallel electron velocity distribution functions
(solid lines) in arbitrary units at two different positions in the system,
near the base and in the supersonic region. Shown are results of the
N=6400 simulation. The dashed lines represent Maxwellian
distributions for which the first three moments (density, mean velocity
and temperature) are those of the measured distributions.
Velocities are normalized to the electron thermal velocity
![]() |
Open with DEXTER |
For the N=6400 simulation the parallel electron velocity
distribution functions at the base of the system at r=r0,
and in the supersonic region at r=45 r0 are
shown in Fig. 5.
Since the collisional
mean free path
is proportional to v4
high energy electrons are nearly unaffected by collisions
on their journey through the system. Thus, the velocity
distribution of the high energy electrons flowing downward is the imprint
of the upper boundary condition whereas the velocity distribution of the
high energy electrons flowing upward is the imprint of the
lower boundary at r=r0. On the other hand, the
low energy electrons, which populate
the core of the velocity distribution function, are strongly
affected by collisions. As a consequence, at low velocities,
the electron velocity distributions are approximately
isotropic Maxwellians which do not carry any heat flux. Instead, the
heat flux is carried by the high energy electrons which are
responsible for the asymmetry of the distribution function.
As illustrated by the profiles in Fig. 5
the heat flux is due to a deficiency of downward
flowing high energy electrons near the lower boundary
and to an excess of upward flowing particles in the
supersonic region.
![]() |
Figure 6:
Electron heat flux profile (solid line) and thermal Knudsen number
![]() ![]() ![]() ![]() |
Open with DEXTER |
Figure 6 shows the heat flux and the thermal Knudsen number KT measured in the N=6400 run.
One observes that while
it remains approximately
constant in the supersonic region above the
level,
KT grows steeply in the subsonic region,
where it increases from 10-3 to 10-2.
Since
in whole simulation domain, it is not
particularly surprising that the heat flux closely follows a
dependence (dot-dash curve) as in the case of the
classical Spitzer-Härm heat conduction formula (Spitzer & Härm 1953).
This conclusion is misleading, since, despite the
smallness of the Knudsen number, the heat flux
is strongly non classical.
As already pointed out by several authors in the past
the classical heat conduction
formulation breaks down either because the heat flux intensity exceeds
a value of the order
10-2q0 (Gray & Kilkenny 1980),
or because the Knudsen
number is larger than 10-3 (Shoub 1983),
or because the flow velocity is a significant fraction of the sound
speed (Hollweg 1974; Alexander 1993), or
even because the electric field in the system
is of the order of the Dreicer field
(Scudder 1996b).
Indeed, in the N=6400 simulation the electric field at the sonic point is
,
reaching
an intensity
in the N=1600 case. Concerning the
heat flux intensity, Fig. 6 shows that
it is small enough for the
the low heat flux intensity condition established by
Gray & Kilkenny (1980) to
be satisfied. One can therefore conclude that
the heat flux intensity is low enough for the plasma
to be capable to support a Spitzer-Härm flux.
Let us now address the problem of the heat flux
in a flowing, and weakly collisional plasma. As discussed by
Hollweg (1974) and Alexander (1993)
a non negligible fraction of the of the energy is carried
by a collisionless term of the form
where
is a positive numerical factor of order unity
(note that the electron temperature has been supposed to be isotropic
by these authors).
Given that collisions are still relatively important in our simulations
we make the ansatz that the observed electron heat flux
is made of the sum of a classical (collisional) term
(e.g., Braginskii 1965)
and a collisionless term
![]() |
Figure 7:
Electron heat flux calculated using the collisionless approximation
![]() ![]() |
Open with DEXTER |
Given the artificially low mass ratio
in our simulations there is a concern about the
sensitivity of the results on the value of
.
In order to address this question we show the same simulation for two different
values of
in Fig. 8.
The other parameters are
identical for both simulations , i.e.
and N=6400.
In both cases the formation of a transonic wind occurs,
in association with the formation of a
maximum in the proton potential. However, the maximum's amplitude is
substantially higher, and less peaked, in the low mass ratio simulation.
The discrepancy is likely due to the fact that in the high mass
ratio case the scattering of the electrons in velocity space
by the protons is more efficient than in the low mass ratio case.
Indeed, the temperature ratio
reaches a value of 3 at the upper boundary in the
case (cf. Fig. 4) and a value of 4 in the
case. As a result the absolute value of second term on the right hand side
of Eq. (9) is significantly
smaller in the high mass ratio case than in the low mass ratio case.
Since the sign of this term is negative it contributes in reducing the
the strength of the overall positive electric field.
From Fig. 4 one may argue that
a similar argument applies to the
observation that the electric field strength increases with increasing
plasma density (i.e. with increasing collisionality) as does effectively
show Fig. 1.
![]() |
Figure 8:
Dependence of the Mach number and the proton potential energy
on the proton to electron mass ratio for the simulation with N=6400. Even though the qualitative behavior is similar it
appears that a high mass ratio is in favor of a stronger
acceleration of the wind. The definitions for the Mach number
and the normalizing energy ![]() |
Open with DEXTER |
Extrapolating these observations
to
and N=6400 one therefore expects the maximum
of the proton potential to drop to an even lower level.
The peak is expected to be at least as marked as for the
case.
From self-consistent kinetic simulation of a solar type wind
we find that, unless an efficient isotropization mechanism
for the electron velocity distribution (e.g. wave-particle interaction),
or some source of suprathermal electron distributions are
invoked (e.g. shock produced), the formation of a transonic wind
is only compatible with a sufficiently
high collisionality in the vicinity of the sonic point
r=r*. In oder words, the coronal density must exceed
a threshold density for the wind acceleration to
be sufficiently strong to become supersonic.
Given the admittedly oversimplifications in our model, combined
with the fact that the parameters we use are rather unrealistic
(excessively high coronal temperature and low
ratio) makes it impossible for us to specify
an upper limit for the thermal Knudsen number at the sonic point r*.
We do merely show that the number cannot be arbitrarily large,
the limiting value most likely being of the order unity, or less.
As already stated, non Maxwellian boundary conditions, feeding
an excess of suprathermal particles into the system (e.g. kappa
distributions) may help overcoming the low Knudsen number condition.
However, the existence of non thermal particle distributions
rises the question of their origin. We chose not to
address this question and assume Maxwellian boundary conditions
which have the advantage of not requiring the introduction
of additional ad hoc parameters into the model.
Thus, in the absence of any electron scattering mechanism
(apart from collisions),
and unless special boundary conditions are imposed
at the base of the simulation domain,
collisions appear to be the essential ingredient for
the wind to become accelerated to supersonic velocities,
mainly because collisions are necessary to convert
the electron heat flux into bulk energy of the plasma.
The enthalpy gradient does also contribute to the
acceleration of the wind. However, the acceleration
associated with the radially decreasing
enthalpy is found to be weakly dependent on the
plasma collisionality and doesn't seem to be the
discriminating factor in the acceleration of the wind
to supersonic velocities.
In simulations where a transonic wind forms we find that the proton
potential has a maximum near (but above) the sonic point.
Typical values of the electric field
near the sonic point r* are found to be of the order
of Dreicer's field, or larger. The presence of such
strong electric field intensities may contribute in making
the electron heat flux depart from
the classical Spitzer-Härm formula (which requires
the electric field being much weaker than Dreicer)
but the main reason
for the observed heat flux to depart from the Spitzer-Härm
prediction is due to the presence of a strong
"non collisional'' heat flux
.
The latter appears to contribute significantly to
the total electron heat flux, even in the region
where the wind velocity is much smaller than the sound speed.
We are aware of the fact that extrapolating the above results to the
"real'' Sun is a perilous exercise. However, we do
not expect the qualitative behavior
of a system with real coronal temperature and real proton to electron
mass-ratio to behave in a substantially different way from the high density
case discussed in this paper. In particular, increasing the
proton to electron mass ratio from 400 to 1836 implies a factor
two increase in the electron thermal velocity, only. Given that
the electron thermal velocity at the base of the system is already
one order of magnitude larger than the escape velocity for the
case, not much difference is expected
in a system with twice this thermal velocity. The skeptical
reader may also argue that using a realistic mass ratio would
substantially modify the transport properties of the plasma.
This is certainly true, but we expect the
modifications to be small, essentially because neither the classical
electron heat flux
nor the collisionless heat flux
depend on the mass ratio, at least as long as
.
We expect the excessively high coronal temperature used in our simulations
to be a more crucial limitation in the process of transposing the above results to the
solar case just because the proton thermal velocity and the escape velocity
are of the same order, i.e.
.
Indeed, the coronal temperature of the real Sun
is roughly 3 times smaller that the value we use here. This
implies a factor
difference for the order unity
quantity
.
The impact is certainly non negligible from
a quantitative point of view but there aren't any reasons for
us to believe that the qualitative aspects of our results
do not apply to the solar case. For instance, whether or not the
collisionless heat flux near the solar sonic point
is really one order of magnitude stronger that the classical
Spitzer-Härm flux remains an open question since this finding
discords with recent results from multi-moment models.
Acknowledgements
We thank the referee for the detailed comments which very much helped us during the revision of the paper.