Astron. Astrophys. 333, 1069-1081 (1998)

3. Model equations

For the description of the transport in plasma a (stationary) five-moment-approximation of Boltzmann's equation will be used (see e. g. Schunk, 1975). In the energy balance heating and radiative cooling are included as parameterized functions of temperature and density. The moment expansion is closed by assuming the pressure (of an ideal gas) to be isotropic, , and that the classical formulation of the heat flux proportional to the temperature gradient can be used (Sect. 3.3). In Sect. 3.4 and 3.5 numerical solutions of the equations are described and the boundary conditions are formulated.

3.1. Multi-fluid transport equations

In the formulation of Marsch et al. (1995) the stationary equations of continuity and momentum for a species j read as

Here and denote the particle density and the velocity of the species j. The sources and sinks for the particle flux density, , are due to ionization and recombination with the respective rates .

In this paper the indices k, j, and are used. Whenever k and j are found, an interaction of the the type is present, i. e. the particles stay the same, e. g. as in elastic collisions. j and k represent different elements. j and k are used also in the case of resonant charge exchange, e. g. , because before and after the interaction the same particles are present, just with exchanged identity. In a reactive interaction (e. g. ionization), the particles before and after the interaction are not the same, but they can change their state of excitation or ionization. Because they are from the same element, the same letter is used, with or without a prime: j and . In contrast to Marsch et al. (1996) we do not use j and because j and are more general: e. g. the continuity equation (1) is valid for both, neutrals and ions. If we would have used j and , we would have to distinguish between the case of neutrals and ions, because there is no symmetry between the ionization and recombination rate. Additionally j and can also be used to describe excitation and de-excitation processes.

The sound speed and the ion-acoustic speed are given by and respectively, with Boltzmann's constant , the temperature and the atomic mass . In the present case of equal temperatures, , these two speeds are equal, . denotes the charge number of the species j and is the gravitational acceleration.

The exchange of momentum between the species is due to elastic collisions and ionization/recombination with the respective rates and (see Sect. 3.2). The influence of the magnetic field can be formulated in an analogous way with a "magnetic" frequency . Here is the gyro-frequency with the elementary charge e and the magnetic field strength B. The direction of the magnetic field is given by . Both the elastic collision frequencies and the "magnetic" frequencies obey the symmetry relations and .

The term in (2) describes the momentum that a particle which is created or destroyed (by ionization/recombination) adds or subtracts to the momentum of the species (see e. g. Geiss & Bürgi 1986). This was neglected in the model of Marsch et al. (1995) because the ionization/recombination rates are some orders of magnitude smaller than the momentum transport due to resonant charge exchange; compare (9) and Table 2 with Table 1. But as it can be seen immediately from (2) this term can become important if the density of either the neutral or the ionized species becomes very low. Thus this term is considered in all numerical solutions.

Table 1. Collisional rates in a hydrogen-helium mixture consisting of neutral and (first) ionized components. and denote the temperature in K and the particle density in m^-3 respectively (from Geiss & Bürgi, 1986 and von Steiger & Geiss, 1989).

Table 2. Boundary conditions used for the hydrogen-helium mixture. The conditions marked with an asterisk, *, are not needed for a pure hydrogen gas. In the calculations the given values are used, unless stated otherwise.

In the derivation of the momentum equation (2) the mass of the electrons was assumed to be much smaller than the mass of the ions, , and quasi-neutrality and zero-current was presumed,

The (ambipolar) electric field is given by , a relation which has been exploited in the derivation of (2).

As outlined in Sect. 2.2 the temperature is assumed to be the same for all species, . For this only one energy equation has to be solved. Following Schunk (1975) the energy equation for the electrons is given by

The sources and sinks in the energy balance are due to

namely elastic collisions with the heavy ions, heating and cooling. Often also the heat flux is comprehended as a source or sink of energy. These processes are discussed in more detail in Sect. 3.3.

3.2. Ionization, recombination and collisions

In the chromosphere the material becomes (first) ionized, and thus in the particle dynamics ionization and recombination play an important role. Additionally elastic collisions are of importance, particularly in the nearly neutral regions of the atmosphere where diffusion may be of relevance.

Ionization

In the atmospheric layer considered here the most important process is photoionization, for which the rate for hydrogen is about (von Steiger & Geiss 1989). This process overwhelms ionization due to electron collisions at chromospheric temperatures of K by some orders of magnitude: following Lotz (1967) the corresponding rate is .

With the photoionization rate and the typical diffusion velocity as mentioned in Sect. 2.1 one can define an ionization length

which is of the order of 50 km. Thus the thickness of the ionization layer of hydrogen is somewhat smaller than the gravitational scale height in the chromosphere (300 km), but not negligibly small.

As in models for the earth's ionosphere (see e. g. Banks & Kockarts, 1973), the chromosphere is assumed to be irradiated from above by UV-photons originating in the transition region and corona. In the considered layer these photons are absorbed in the ionization process. As a consequence, the ionization rate will vary with depth in the atmosphere. The below presented description of this variation of the ionization rate follows the one given in the analytical model of Peter & Marsch (1997). For a further discussion of the assumptions leading to (7) see their paper. The change of the ionizing radiative flux over the distance is proportional to the flux itself and the density n of the absorbing material: . If the cross section for photoionization is assumed to be nearly constant over the relevant wavelength band, then one can apply the same relation to the ionization rate: if s is the vertical coordinate, the respective photoionization rates of hydrogen and helium are determined by

where the indices n and i stand for the neutral and ionized components of the respective elements. This corresponds to the processes as found in the earth's ionosphere (see e. g. the textbook of Banks & Kockarts 1973). The latter equation for the ionization rate has to be solved together with the transport equations (1), (2) and (4).

The cross sections for photoionization as needed in the description of the ionization rate (7) are taken from Vernazza et al. (1981). The mean values (see Peter & Marsch 1997) can be estimated as

This treatment of the ionization rate and photon flux is just an approximation: on the one hand the cross sections and the photon flux depend on the wavelength, on the other hand the exited states of the hydrogen atom are neglected. Thus in a complete atmospheric model one would have to solve not only two continuity and momentum equations in the case of a pure hydrogen gas, but 26 if e. g. a 12-level hydrogen atom is described (two for every level and two for the protons). Additionally one would have to solve properly the rate equation for the photon flux as a function of wavelength.

This has been done e. g. by Vernazza et al. (1981) or Fontenla et al. (1990). But they solved these equations just for the single-fluid hydrostatic case, i. e. for zero velocities, and thus did not have to solve the momentum equations. In contrast, in the present paper the emphasis is on the diffusion and flow of the material. To avoid the problem of combining the radiative transport with the plasma dynamics we use the practical approach described above, where we consider a "mean" or "effective" hydrogen atom. A more complete model including radiative transport would probably result in a (slightly) different profile of the ionization rate with depth, but this should not change the main results for the plasma dynamics.

Recombination

The respective rates for the recombination are simply taken from von Steiger & Geiss (1989).

Here denotes the temperature in K and the particle density in m^-3.

Elastic collisions

As this paper deals with a mixture of hydrogen and helium, each being either neutral or (singly) ionized, the following processes in the elastic collisions have to be considered: hard sphere collisions (between neutrals), induced dipole interaction (ions and neutrals) and Coulomb-collisions (ions). The respective collisional rates are derived in various textbooks, e. g. in Burgers (1969). For the collisions between the neutral and ionized particles of hydrogen as well as helium the resonant charge exchange (RCE) is the most important process. The respective cross sections can be found in Banks (1966).

In this paper the (parameterized) collisional rates as given in Geiss & Bürgi (1986) and von Steiger & Geiss (1989) are used. These are listed in Table 1.

Energy sources and sinks

The energy sources and sinks are of great importance for the energy balance in the lower atmosphere. The most important processes are mechanical heating and radiative losses. The loss of kinetic energy (of the electrons) due to collisional ionization can be neglected in the lower chromosphere, because it is not efficient at those low temperatures. However, in a model which includes the transition region this process would have to be taken into account.

Heating

The exact determination of the heating rate is an unresolved problem. For that reason mostly parameterizations are used.

Following Ulmschneider & Kalkofen (1977) the heating in the lower chromosphere is due to the damping of shocks. Priest (1982) gave some heuristic arguments that in this case the heating rate behaves like and is exponentially damped. This result corresponds to the behaviour of H as described by Rosner et al. (1978).

With parameter values taken from the references mentioned above, giving a damping length of about 3000 km and an energy flux of at the bottom of the chromosphere, the resulting heating rate is

Once more, and are the density in and the temperature in K respectively.

Radiative cooling

Besides heat conduction radiative cooling is the most important energy loss process. Although the radiation is emitted by the atoms and ions, this process is a loss mechanism for the energy of the electrons, since the excitation of the atoms and ions is due to electron collisions and thus the required energy is taken from the kinetic energy of the electrons.

If the radiation is not treated self-consistently, an approximate description is required. For this purpose the radiative loss function L is often assumed to be a power law of the temperature. A great variety of work has been done on this subject.

The description in the present paper will follow the work of Peres et al. (1982). For chromospheric temperatures they give the following parameterization:

As before, the indices ₁₆ and ₄ indicate that the density has to be taken in units of and the temperature in K.

It should be noted that this treatment is a little problematic, because the main assumption going into (11) is to describe an optically thin plasma. But this may not be true at low temperatures. However, Kuin & Poland (1991) showed that for the optically thick plasma of a loop a similar (though corrected) radiative loss function can be used.

Heat flux

Instead of describing the heat flux vector by a heat flux moment equation, it is assumed to be proportional to the temperature gradient, i. e.

Following Geiss & Bürgi (1986) this long known relation can be derived by assuming a subsonic flow, an isotropic pressure and time stationarity. In this case the heat flux equation for the electrons can be reduced to

with the dimensionless number

For a plasma consisting only of neutral and singly ionized particles, as considered in the present paper, the quasi-neutrality renders a constant, .

Using the value for the collisional rate as given in Table 1, the simplified heat flux equation (13) leads directly to the relation (12) with the thermal conductivity along the magnetic field being

This relation together with (12) is also known as Spitzer's law.

It is just interesting to note that this conductivity in the chromosphere at K is of the same order as in the air on the earth at room temperature!

In this paper the conductivity across the magnetic field is not considered, because it is some orders of magnitude smaller than . Given the strong magnetic fields of 10 - 100 Gauss, the electrons are strongly magnetized and despite of collisions constrained to move along field lines.

3.3. Formulation of the transport equations appropriate for numerical treatment

In this section the transport equations (1), (2) and (4) will be formulated for one spatial dimension along the magnetic field (see Sect. 2.1) and in a way appropriate for the numerical treatment. As the aim is to use standard numerical routines solving a system of differential equations of first order, the equations have to be written in the form

where stands for the densities, velocities and temperatures of the respective elements.

Here as in the following the prime, , denotes the derivative with respect to the vertical coordinate s. Concerning the numerical treatment it is better to use the particle flux densities instead of the velocities, where is the component of the velocity along the coordinate s. This is because for very low densities the flux remains finite while the velocity can become extremely high due to numerical errors.

First of all the energy equation (4), which is of second order, will be divided into two equations of first order. For this purpose a new variable , with from (15), is introduced,

Because and can be written in terms of the variables but not their derivatives (see below), the energy equations thus attain the required form (16).

Before re-formulating the continuity and momentum equations the following abbreviations are introduced:

The derivative of the squared sound speed is given by . All these abbreviations are functions of the variables itself but not of their derivatives!

Now the equations of continuity and momentum can be written in compact form as

Because of the quasi-neutrality condition (3) this system of equations can always be written in the standard form (16).

In the case of a pure hydrogen plasma the system (19) and (20) reads

It is straight forward to obtain the system for the hydrogen-helium mixture. The corresponding eight equations are just a little longer.

It is easily proven, that for vanishing helium components, i. e. the densities and fluxes of He and He are set to zero, this system results in the system (21) for a pure hydrogen gas.

It should be noted that in both systems, in (21) and (22), no use was made of the conservation of flux. It would have been easy to replace one of the two continuity equations for each element respectively. But as one aim of this model is to determine the total flux for every element, the full system is solved, while no absolute number of the flux is given as a boundary condition (see Sect. 3.5).

For the complete model the ionization equations (7), one for each element, and the energy equations (17), (18) have to be solved together with the continuity and momentum equations (21) or (22). This means that in the case of a pure hydrogen gas seven first order differential equations have to be solved. In the case of a hydrogen-helium mixture this number is twelve.

As this system has the explicit form (16), numerical standard routines to solve a system of ordinary differential equations can be applied. In the present case the routine D02RAF from the NAG-library was used. This routine allows the supply of an approximate first-guess solution, automatic addition of grid points and a flexible formulation of the boundary conditions. To solve the differential equations it uses a deferred correction technique and a Newton iteration.

3.4. Boundary conditions

To solve the above mentioned first order differential equations, in the case of a pure hydrogen gas seven, or for a hydrogen-helium mixture twelve boundary conditions are needed. These are chosen from a physical point of view.

The lower boundary of the considered chromospheric layer is placed at about 8000 K, well above the temperature minimum.

Here the particle density is assumed to be of the order of some . This value is the often used, e. g. in Vernazza et al. (1981). The abundance of helium at the bottom is assumed to be the same as in the photosphere; this value is known quite exactly (e. g. Anders & Grevesse, 1989). The degree of ionization at the bottom is calculated from the Saha equilibrium.

At least the ionized and the neutral component of every element should enter the layer from below with the same velocity, which is due to their tight "collisional" coupling in the form of charge exchange or ionization-recombination balance. Please note, that this does not mean that the different elements must have the same velocity. There can (but does not have to) be diffusion at the bottom of the chromosphere.

The upper boundary of the layer is placed about 1000 km above the bottom. The numerical calculations have shown that for a thicker layer the results remain unchanged.

At the top the ionization rate is fixed. This means that the photon flux is given at the top. Here the flux values calculated by von Steiger & Geiss (1989) from a solar spectrum are used.

The relative velocity of the neutral and ionized component, e. g. , is assumed to be constant at the top. This was used instead of the stronger condition of equal velocities, e. g. , because of numerical reasons: in the case of a bad first-guess solution the solving routine was much more stable when using the weaker condition. But in the end the numerical solution results in invariably equal velocities at the top, e. g. (see Fig. 3).

The heat flux at the top of the layer, where the temperature reaches about K, is chosen according to Hansteen et al. (1993), leading to about 0.01 W/m² by using (15). This value is somewhat arbitrary, as it is not well known, but only calculated in models and not measured directly. However, the numerical calculations have shown that the result does not depend on the exact value of the heat flux given at the top. If it is chosen "not correctly", this results in a boundary layer for the temperature at the top, i. e. the temperature rises or drops rapidly on the last km. For this reason the heat flux is chosen to give a smooth temperature profile at the top.

In Table 2 the boundary conditions are summarized and the respective values as used in the present paper are given.

It should be noted that no absolute value for the particle flux or the velocity is given as a boundary condition or imposed as a constant of integration. The only absolute values given as a boundary condition are densities, temperature and ionization rates. Especially no absolute values of the velocities are given. Thus the particle flux through the layer will result from the boundary conditions in terms of ionization rates on the top and densities at the bottom of the layer.

© European Southern Observatory (ESO) 1998

Online publication: April 28, 1998

helpdesk@link.springer.de