2.1. Background Theory and Techniques
A simple way to simulate the evolution of a set of interacting species is to apply Equation (
2) to each reaction (
1) for a time-step
, so that the change in each species is:
For each species
i, the gains and losses are added up to give the total gain
and total loss
so the new density
of species
i after time
is:
The development of the densities of all species can be simulated by iterative application of Equation (
6) over the required time, but the magnitude of
is limited by the requirement that the density must not go negative and so this “explicit” formula is unsuitable for simulation of systems over long time intervals.
An alternative “semi-implicit” method [
15] is justified by the fact that, through the time interval
, the loss rate
is proportional to the instantaneous density
. Hence, an approximation is being made irrespective of whether the value of
L is approximated by being proportional to
or
. Thus, substituting the final density for the initial density:
and rearranging:
leads to:
As the new density cannot be negative, Equation (
9) avoids the main problem with Equation (
6) and so allows much longer values of
to be used.
There are more sophisticated time-step algorithms (e.g., [
15,
16]), but the requirement here is to be able to discriminate between the error due to the time-step algorithm and any error inherent in the approximation of the electron energy spectrum. For example, the Gauss–Seidel method is quoted to have an error of 1% [
16]. As Equation (
9) would be expected to have no error in the limit of very small time steps, the errors due to the simulation method can be separated from the errors due to the time-step calculation.
To allow the simulation to run to equilibrium with minimum computation time, an adaptive time step
is required to implement Equation (
9) efficiently [
3]. The initial value of
is set very small (e.g.,
s) and is then successively increased as:
where
is the minimum time interval for any of the constituents
i to fall to zero in the next time step
at the current rate of change, and
f is a fraction that acts to reduce the error in the calculations.
The rate constant for collisions between an electron and a gas molecule in a unit volume is
, where
v is the electron’s velocity and
is the cross section or probability for the interaction. Thus, the rate constant as a function of energy is:
It therefore follows that the rate constant
for all transitions of electrons starting in a range [
] is:
where
is the electron energy distribution function (EEDF). For a Maxwell–Boltzmann distribution:
where
T is the electron temperature, and
is Boltzmann’s constant [
17].
The cross section
for a superelastic collision, for electron impact energy
E, can be determined from the inelastic cross section using [
18,
19]:
where
is the threshold energy of the inelastic excitation.
Electrons impacting atoms and molecules will gain or lose energy by elastic scattering. Published cross sections [
20] for elastic electron scattering by OH are only for electron energy above 1 eV. Thus, as an approximation, the formula for the elastic electron energy transfer rate for O of Banks [
21] is used, i.e.,
where
is the energy transfer rate (eV cm
s
),
is the temperature of a Maxwellian distribution of electrons,
T is the temperature of the O atoms,
is the electron density (cm
), and
n(O) is the O density.
2.2. Proposed Method
To incorporate the electron reactions (
4) into a time-step simulation, the proposed approach is to split the electron energy range
into
N subranges R
–R
, with R
set as the range
.
In place of the single reaction
, a series of reactions:
is entered into the list of chemical reactions, with the number density of electrons in each energy range R
treated in the same way as the density of a chemical species. (For example, in the current implementation, the number of electrons in the range [0,
] is stored in a variable R
0. While a logarithmic distribution of subranges is desirable, a linear distribution is considered first for simplicity before proceeding to the complications added by a logarithmic distribution.
Consider a case where electrons lose or gain energy
in elastic and superelastic collisions with gas molecules. In
Figure 1a, the case where the transition energy
is less than the size of the energy subranges is considered for two energy subranges with boundaries at
,
and
and midpoints at
and
. Applying Equation (
12) with
(as the electron–energy distribution is varying), the rate constant
for transitions where the electron crosses from R
to R
is:
The rate constant
, for cases where the electron energy remains in the same subrange while the excited species is produced is similarly:
If
were constant, then energy would be conserved because the energy lost by electrons that remain in the same subrange would be offset by the higher energy implied for those that transition to the next subrange, i.e.,
However, as
varies across the subrange, energy will not be conserved so a correction is necessary. This is made by solving Equation (
19) for a modified value of
:
In the case of
Figure 1b, where
is greater than the subrange size, all electrons transition to a lower subrange and there is no physical value of
. However, a notional value of
can be calculated to maintain conservation of energy:
While
is an unphysical quantity and can be negative, it implements conservation of energy by correcting for the approximation that
is constant over the subrange. Applying a similar analysis for the case where
is much larger than the subrange size leads to a general equation:
for inelastic collisions and
for superelastic collisions.
As electrons may be initially introduced at high energy, but then proceed to lose energy through a series of collision processes down to a very low energy, it is desirable to make the energy subrange boundaries on a logarithmic scale to keep the total number of subranges
to a viable minimum for the calculations while allowing for adequate resolution at low electron energies. Equations (
22) and (
23) are applicable, even where electrons in one subrange can transition to several lower subranges.
Figure 1.
Electron energy losses (arrows) in relation to electron–energy subranges: (a) For an energy transition that is less than the subrange size , electrons with initial energies in the shaded area cross the boundary to range R, while the others remain within range R. (b) For greater than the subrange size , all electrons transit to a lower subrange, with those in the shaded region going to R and the others to R.
Figure 1.
Electron energy losses (arrows) in relation to electron–energy subranges: (a) For an energy transition that is less than the subrange size , electrons with initial energies in the shaded area cross the boundary to range R, while the others remain within range R. (b) For greater than the subrange size , all electrons transit to a lower subrange, with those in the shaded region going to R and the others to R.
Aspects of the model are illustrated in
Figure 2, for inelastic and superelastic collisions of electrons with OH molecules in the lowest and first vibrational level of the ground state. It shows rate constants for individual energies and average rate constants for the simulation where the electron energy range is divided into 30 logarithmically-spaced energy subranges between 0.01 eV and 20 eV, labelled R
–R
, with electrons below 0.01 eV assigned to the subrange R
.
Inelastic electron–impact cross sections for the 0→1 vibrational excitation in OH, calculated using the method of Riahi et al. [
22], are shown by a solid curve. This case was chosen as it provides a smoothly varying curve, so that the accuracy of the simulation can be evaluated with minimum contribution from any structure in the cross sections. The superelastic cross sections for the
transition, calculated using Equation (
14), are shown by the dashed curve. Using these cross sections in Equation (
11), rate constants for the
excitation (
eV) are shown for 2500 logarithmically-spaced initial energies (between 0.005 and 20 eV) by horizontal green lines drawn between the initial and final electron energies at the appropriate vertical position. Rate constants for the superelastic deexcitations are plotted similarly in purple. Transitions ending in the lower-energy ranges all start from within a small energy range near the threshold, so, in applying Equation (
12) to calculate the averaged rate constants, the steps in energy must be sufficiently small.
Equations (
17), (
18), (
22) and (
23) were applied to calculate the rate constants for transitions between pairs of subranges. These are represented in
Figure 2 by left-pointing arrows for inelastic collisions and by right-pointing arrows for superelastic collisions, drawn from the centre of the initial subrange to the centre of the final subrange, with the vertical position showing the value of the average rate constant. At low electron energies, the average rate constants are much lower than the individual ones because, particularly near the threshold, only a fraction of the possible transitions starting within the higher energy subrange end in each lower-energy subrange. Thus, these are better described as “effective” rate constants. This is an example of the way in which details of the cross section within an energy range are embedded in the rate constants. At higher energies, where the transition energy is much less than the width of the energy subranges, the effective rate constants are again lower because only transitions that cross a subrange boundary are included. The effective rate constants for collisions that leave the electron in the same subrange are plotted as octagons for the inelastic collisions and squares for the superelastic collisions.
In some cases, the corrected rate constants produced by Equations (
22) and (
23) can be negative. While this is unphysical, it makes a necessary correction to maintain conservation of energy. The negative values are indicated by filled symbols in
Figure 2.
As examples, the reactions added to the model for inelastic collisions starting in subranges R
and R
are specified as:
A proof-of-concept test of this proposed method, to investigate validity and accuracy, is to simulate the injection of electrons of one energy into a gas of OH molecules that are initially all in the lowest level of the ground state and to consider only excitation to the first vibrational level. While this is unphysical, as it does not consider the electron–electron interactions that would dominate in this case, the rationale is to make a stringent test in a case where a theoretical equilibrium can be calculated. Initially, the electrons will lose energy in exciting the OH molecules, then regain some of it via superelastic collisions, until an equilibrium is reached. The total energy of the system (of electrons plus excited OH molecules) should be equal to the total input energy, while the electrons should have a Maxwell–Boltzmann distribution as given in Equation (
13).
Consider that
electrons with a total energy
are mixed with
OH molecules in the lowest level of the ground state. Assuming that at equilibrium the vibrational temperature of the gas and the electron temperature are the same, this temperature
T can be found by solution of the equation:
In order to compare results for logarithmic and linear spacing of the subranges, where the average energy of the electrons
(i.e., centre of subrange R
) is different between the two,
is chosen so that the equilibrium gas energy
is the same:
The electrons can also gain or lose energy in elastic collisions with the gas molecules. To calculate rate constants for elastic collisions, Equation (
15) for electron scattering from O is used, making the further approximation of applying the characteristic temperature
of a Maxwellian distribution of electrons to that of a single electron. This gives the change of energy of an electron with energy
E (eV) for unit gas and electron densities as:
where
. As
is a small fraction of the electron subranges, the rate for transfer from subrange R
to R
is divided by the number of individual collisions required to transfer a total energy of
, where
is the midpoint of the energy subrange R
. Thus, the rate constant for elastic collisions that transfer energy from range R
to range R
is:
To verify this method, the simulated energy of the electrons can be compared with that predicted in an iteration of the total electron energy
:
where
is calculated using Equation (
15).