Inclusion of Electron Interactions by Rate Equations in Chemical Models

: The concept of treating subranges of the electron energy spectrum as species in chemical models is investigated. This is intended to facilitate simple modiﬁcation of chemical models by incorporating the electron interactions as additional rate equations. It is anticipated that this embedding of ﬁne details of the energy dependence of the electron interactions into rate equations will yield an improvement in computational efﬁciency compared to other methods. It will be applicable in situa-tions where the electron density is low enough that the electron interactions with chemical species are signiﬁcant compared to electron–electron interactions. A target application is the simulation of electron processes in the D-region of the Earth’s atmosphere, but it is anticipated that the method would be useful in other areas, including enhancement of Monte Carlo simulation of electron–liquid interactions and simulations of chemical reactions and radical generation induced by electrons and positrons in biomolecular systems. The aim here is to investigate the accuracy and practicality of the method. In particular, energy must be conserved, while the number of subranges should be small to reduce computation time and their distribution should be logarithmic in order to represent processes over a wide range of electron energies. The method is applied here to the interaction by inelastic and superelastic collisions of electrons with a gas of molecules with only one excited vibrational level. While this is unphysical, it allows the method to be validated by checking for accuracy, energy conservation, maintenance of equilibrium and evolution of a Maxwellian electron spectrum.


Introduction
Simulations of atmospheric processes involving chemistry often involve the processing of a set of reactions specified in the form: where A, B, C and D are atoms or molecules and k is the rate constant. The rate r at which A and B interact is given by: where [x] represents the density of species x and the rate constant k has the units L 3 s −1 where "L" is the unit of length in which the density is specified (e.g., cm 3 s −1 ).
For each species i with number density n i , a one-dimensional mass continuity equation can be written: where P i and l i are the production rate and loss probability [1].
The set of continuity equations can be solved iteratively until equilibrium is obtained [2], or applied over a series of time steps for non-equilibrium simulations [1]. Even where equilibrium is required, a time-step simulation is useful to determine the time for equilibrium to be attained [3].
Electron interactions can be (e.g., [4]) incorporated into the set of rate equations in forms such as: where * represents a change in energy and k r , k a and k ex are the average rate constants, for recombination, attachment and excitation, for electrons with a Maxwellian distribution characterised by an electron temperature T e . However, in many atmospheric applications, such as those including auroral, cosmic-ray or VUV input, the electrons are not in thermal equilibrium [5,6]. For example, in the Earth's nighttime mesosphere, electrons are created in ionization produced by cosmic rays and by Lyman-α radiation [7], and they then lose energy in collisions with atmospheric molecules before recombining with positive ions or attaching to molecules to produce negative ions. At the same time, chemical processes produce vibrationally-excited OH [8], with some of this energy being transferred to other species and then to free electrons in superelastic collisions. The time-scale of these processes depends on the density of each constituent. Therefore, calculating the electron density requires a non-equilibrium simulation that includes both the electron-impact processes and chemical processes. Another potential application may be to enhance Monte Carlo simulation of nonequilibrium electron-liquid transport [9]. In cases where there are too many tracks (from ionisation) or too few tracks (from attachment) to simulate, approximate approaches such as "re-scaling" [10] or "weighting" [11] are used and the distribution of the excited species is not considered. It is possible that using a chemical model that includes electron interactions could be applied in such situations. In addition, the chemical model could be applied to determine a non-equilibrium distribution of excited species in a gas prior to applying a Monte Carlo simulation of electron impact.
Ionizing radiation applied to biological materials can produce secondary electrons, radicals (e.g., OH • ) and ions (e.g., H 3 O + ), all of which can cause damage to DNA [12,13]. As the electrons may lose energy in a series of interactions with water and other molecules before interacting with a DNA molecule, a chemical model that incorporates a changing electron energy spectrum should be applicable to simulating damage to DNA by ionizing radiation.
Thus, an aim of this work is to develop a method to calculate the development of the electron energy spectrum for a system where the electrons are interacting with atoms and molecules that are simultaneously interacting with each other. This has previously been addressed in plasma physics as part of a "state-to-state" approach [14] in which excited states of atoms and molecules are treated as independent species. This chemical model is coupled to the electron distribution via the Boltzmann equation, with the reaction rates between electrons and chemical species being recalculated each time the coupling is made. The aim here is to do this with a different approach in which groups of electrons of similar energy are treated as individual chemical species and consequently all interactions are calculated in each time step. An advantage of this method is that the electron interactions can be added into the list of rate equations in existing models, avoiding the work involved in merging different computer codes that are geared to different problems. Another advantage is that some detail (e.g., thresholds and narrow resonances) of the electron interactions can be incorporated into the reaction rates, so they do not have to be recalculated at each time step and so a relatively coarse energy grid can be used to reduce computation time.
The proposed method is to divide the electron energy range into subranges and to treat each subrange in the same way as a chemical species i.e., each reaction, such as those in Equation (4), is replaced by a series of reactions that represent the interaction for a set of different initial and final electron-energy subranges. Due to the large range of energies of different processes (from ionisation to recombination), a logarithmic distribution of the energy ranges is preferable.
A 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, where only vibrational excitation of one level of OH is considered. OH is chosen as it is relevant to two of the examples postulated above.
In Section 2.1, the background theory and computational techniques required to develop and test the proposed method are outlined. In Section 2.2, the method and its computational implementation are outlined. The results of various tests of the method to show its viability and the errors involved are given in Section 3. These results are discussed in Section 4, from which conclusions are drawn and given in Section 5.

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 ∆t, so that the change in each species is: For each species i, the gains and losses are added up to give the total gain G i and total loss L i so the new density n i of species i after time ∆t 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 ∆t 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 ∆t, the loss rate l i is proportional to the instantaneous density n i (t). Hence, an approximation is being made irrespective of whether the value of L i is approximated by being proportional to n i (t) or n i (t + ∆t). 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 ∆t 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 ∆t is required to implement Equation (9) efficiently [3]. The initial value of ∆t is set very small (e.g., 10 −10 s) and is then successively increased as: where ∆t min is the minimum time interval for any of the constituents i to fall to zero in the next time step ∆t j+1 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 vσ, 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 k(E a , E b ) for all transitions of electrons starting in a range [E a , E b ] is: where F(E) is the electron energy distribution function (EEDF). For a Maxwell-Boltzmann distribution: where T is the electron temperature, and k B is Boltzmann's constant [17]. The cross section σ s 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 dU e /dt is the energy transfer rate (eV cm −3 s −1 ), T e is the temperature of a Maxwellian distribution of electrons, T is the temperature of the O atoms, n e is the electron density (cm −3 ), and n(O) is the O density.

Proposed Method
To incorporate the electron reactions (4) into a time-step simulation, the proposed approach is to split the electron energy range [E min , In place of the single reaction A + e k −→ A * + e * , a series of reactions: is entered into the list of chemical reactions, with the number density of electrons in each energy range R i 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, E min ] 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 E 01 , E 12 and E 23 and midpoints at E 1 and E 2 . Applying Equation (12) with F(E) = 1 (as the electron-energy distribution is varying), the rate constant k 21 for transitions where the electron crosses from R 2 to R 1 is: The rate constant k 22 , for cases where the electron energy remains in the same subrange while the excited species is produced is similarly: If σ(E) √ 2E/m 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 σ(E) √ 2E/m 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 k 22 : 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 k 33 . However, a notional value of k 33 can be calculated to maintain conservation of energy: While k 33 is an unphysical quantity and can be negative, it implements conservation of energy by correcting for the approximation that σ(E) √ 2E/m 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 N + 1 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. 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 1 -R 30 , with electrons below 0.01 eV assigned to the subrange R 0 .
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 1 → 0 transition, calculated using Equation (14), are shown by the dashed curve. Using these cross sections in Equation (11), rate constants for the 0 → 1 excitation (ε = 0.443 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.

Rate constants
Effective rate constants Rate constants Effective rate constants 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 N e electrons with a total energy E tot are mixed with N g 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 E N (i.e., centre of subrange R N ) is different between the two, N e is chosen so that the equilibrium gas energy E g 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 T e 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 T e = E/k. As ∆E is a small fraction of the electron subranges, the rate for transfer from subrange R j to R k is divided by the number of individual collisions required to transfer a total energy of E j − E k , where E m is the midpoint of the energy subrange R m . Thus, the rate constant for elastic collisions that transfer energy from range R j to range R k 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 E e : where ∆E is calculated using Equation (15).

Results
The 119 effective-rate equations illustrated by the arrows and symbols in Figure 2 were applied to a gas containing (per cm 3 ) 10 8 OH molecules in the ground state and 562,666 electrons in range R 30 (determined using Equation (26) to given a final OH energy at equilibrium of 9.7 MeV) for a logarithmic distribution of subranges. The simulation was run (with f = 10 −6 ) for 1000 s, giving the results shown in Figure 3 for the OH, electron and total energies as a function of time. The simulated values at every thousandth time step are shown by symbols. In addition, the electron distributions at the beginning and the end of the simulation are shown. The theoretical and simulated equilibrium values are plotted as horizontal lines over each energy interval. It can be seen in Figure 3 that energy is transferred from the electrons to the gas molecules, reaching close to an equilibrium after 700 s. At this stage, both the gas energy and the total energy are lower by about 2%, while the electron energy is close to the expected value. The simulated electron energy distribution is close to the calculated Maxwellian distribution ("Equil. ED") but with a residual high-energy tail.
The results for repeating the simulation with 506,677 electrons and a linear distribution of 31 subranges are shown in Figure 4. The simulated OH energy reaches about ∼1% of the equilibrium value after about 400 s, but the electron energy stabilises at (proportionally) a significantly higher absolute value, resulting in a slight excess in the total energy. The excess electron energy appears as higher values in the number of electrons at <1 eV in the electron distribution.  In Figure 5, the simulated OH energies are plotted as a function of time for linear and logarithmic spacing of 31 subranges, each for the uncorrected rates (Equations (17) and (18)) and for the rates corrected for energy conservation (Equations (22) and (23)). While there is little difference until about 150 s into the simulation, the values in the corrected cases approach the predicted equilibrium while those for the uncorrected cases rapidly gain excess energy. The energy in the corrected linear case is constant between 500 s and 1,000,000 s, while in the logarithmic case it is constant from 700 s to 500,000 s and then declines slightly. In addition, the OH energies for the corrected logarithmic case with the negative rates set to zero are shown. These show excess energy appearing after about 1000 s.
The simulations summarised in Figure 5 were repeated for N = 100, giving the energies plotted in Figure 6 as a function of time. In both the corrected linear and logarithmic cases, the OH energies converge to the predicted value and the equilibrium is then maintained for 1,000,000 s, while for the uncorrected cases the energies still diverge, or go to zero in the case of setting negative rates to zero. The rise time for the logarithmic case is slightly longer than for the linear case.     The simulated electron distributions, at the times shown by the vertical dashed lines in Figure 7, are plotted along with the theoretical equilibrium distribution in Figure 8. The errors (being the difference between the simulated and theoretical distribution as a percentage of the higher value) are plotted as crosses where the simulated values are higher and circles for where they are lower. At 110 s, most of the electrons are still at high energy, but, by 300 s, the Maxwellian distribution has appeared, with a remaining high-energy tail. This tail, which is almost gone by 1000 s, has disappeared entirely by 10 6 s, at which time the rest of the simulated distribution is unchanged since 1000 s, seen clearly in the similarity of the details of small-scale discrepancies. A possible reason for this remaining small-scale structure is an aliasing effect between the fixed value of ε and the changing sizes of the subranges that result in varying proportions of electrons transferring to a lower or higher-energy subrange.  Table 1, which gives the consequent computational parameters (the number of rate equations N eqn , the number of timesteps in the simulation N ts and, as an indicator of the computational load, their product) and the errors at 10 6 s in the electron energy, OH energy and total energy. The percentage errors in the total energy are plotted in Figure 9 as a function of the computational load (N eqn N ts ).  In Figure 10, rate equations based on Equation (28) are added to the model to include elastic scattering. The peak OH energy produced by excitation is reduced by about 10%. The OH energy then declines as it is transferred to the electrons and then to the thermal energy of the gas by elastic scattering. As a check on the approximation embodied in Equation (28), Equation (15) is applied at each time step to calculate the energy loss using the instantaneous electron energy, assuming a Maxwellian distribution. The total energy from the simulation ("Sim. total energy") and that from assuming a standard Maxwellian distribution ("Total energy-Maxwellian") are close, particularly in the first 100 s, indicating that the approximations involved in Equation (28)

Discussion
In Figure 6, it is shown that, with the energy range split into a sufficiently large number of subranges, the simulated OH energy approaches close to the predicted value and is maintained at that value over a long time. In Figure 7, it is shown that the OH, electron and total energies all converge to the predicted values, while Figure 8 shows the simulated electron distribution to be close to the theoretical distribution. Thus, Equations (22) and (23) and their implementation are confirmed to be valid. This is emphasised in Figure 6 by the gross departures from equilibrium in cases where the corrections for conservation of energy are not applied.
Comparison of Figures 5 and 6 shows that, while the rise time is longer for both the linear and logarithmic corrected cases for a smaller number of subranges, the effect is more pronounced in the logarithmic case. As this lag in rise time is not seen for the uncorrected logarithmic case with N = 30, where it is close to the values of the linear case until the predicted equilibrium value is reached, it can be inferred that the lag is associated with the correction for energy conservation and that it is more pronounced in the logarithmic case.
As a major rationale for this method is to increase computational efficiency, it is essential to reduce the number of subranges to a minimal value. This introduces inaccuracies as shown in Figures 3 and 4. In the linear case, the electron energy is proportionally high, presumably due to the major part of the electron distribution being averaged into one subrange. In the logarithmic case, the OH and total energy are both slightly reduced and the rise time is substantially increased. Thus, a choice must be made between using a logarithmic or a linear spacing, depending on whether rise time or accuracy of the electron spectrum is more important. An approach to this may be to run both a linear and logarithmic model, so that the difference in results gives an indication of the magnitude of the errors.
The near-accurate calculation of the OH energies in Figure 4 demonstrates the successful incorporation of the detail of the electron-impact cross sections into a few rate equations for transitions in and out of the one subrange [0.01-0.8] eV, which incorporates almost all the electron-impact cross sections below the peak value. This demonstrates the potential of this method to reduce computation times, relative to recalculating the excitation rates for each iteration of the electron spectrum.
In Figure 9, it is seen that, for N = 100 and f < 10 −6 , very high accuracy in the equilibrium total energy is maintained to 10 6 s. For each value of N, the discrepancy asymptotes to a fixed value with decreasing f , showing that the time-step method is not a source of error in the results for f < 10 −6 . The figure also shows that there is a substantial increase in accuracy in reducing f to 10 −6 , but a further decrease makes little difference to the accuracy while substantially increasing the computational load. Thus, to achieve a desired level of accuracy for this particular application, the strategy would be to use f = 10 −6 and then choose the appropriate number of subranges. The inaccuracy introduced by reducing N is substantially greater than by increasing f . Thus, it appears that an optimum value of f , in a compromise between computation time and accuracy, is 10 −5 . Figure 10 shows that elastic scattering (calculated using values for O atoms) reduces the initial excitation of OH by about 10%. Thus, elastic scattering for all the species present must be included in a simulation.

Conclusions
A method to simulate nonequilibrium interactions of electrons with gas molecules was proposed and tested. In this method, the energy range of the electrons is split into subranges that are then treated in a time-step calculation in the same way as chemical species, so the electron interactions can be incorporated easily into existing simulations without new coding being required. It was found that, in excitation of gas molecules with one vibrationally-excited level, the initial energy of the electrons was transferred to the gas molecules until an equilibrium was reached that, with sufficiently small subranges, was very close to the predicted equilibrium values. This equilibrium was then maintained over a long time (10 6 s), validating the method of calculating the rates for the electron interactions. The simulated electron spectrum was also very close to the predicted Maxwellian distribution. On reducing the number of energy subranges (which is desirable for reduced computation time), two discrepancies became apparent. For a linear spacing of the subranges, the simulated electron distribution had a higher energy, presumably because most of the equilibrium electron distribution was represented by a single subrange. However, the accuracy of the results despite this low resolution confirmed the potential of the method to reduce computation times by incorporating details of the cross-section spectrum into a small number of rate equations. For a logarithmic spacing, the rise in gas energy was slow, and the final gas and total energies were slightly less than the predicted values. Thus, the proposed method is capable of producing accurate results, but the minimum number of subranges, and thus computational efficiency, will need to be assessed for the requirements and circumstances of the particular application.