Monte Carlo analysis of field-dependent electron avalanche coefficients in nitrogen at atmospheric pressure Free

Gas discharge physics and the evolution of a plasma in the presence of high electric fields are important topics and have been discussed in numerous studies.^1–7 Interest in this area is driven by the many plasma applications that include anisotropic etching,^8,9 thin-film deposition,¹⁰ pollution control involving ozone depletion substances and volatile organic compounds,¹¹ biomedical treatments and plasma medicine,^12–14 plasma-assisted combustion,¹⁵ pumping of gas lasers,¹⁶ thrusters for space,^17,18 fusion technology,¹⁹ plasma electronics,²⁰ etc. Most of the plasma-based phenomena are crucially dependent on the space-time evolution and energy distribution of the electrons. Hence, it becomes important to understand and predict these characteristics, from which appropriate transport and rate coefficients can be obtained.

Modeling for charged particle transport in non-equilibrium plasmas is a very mature field, and there have already been excellent reviews and discussions on this subject ranging from approaches that solve the Boltzmann equation^21,22 to kinetic particle-in-cell analyses.²³ Early work by Kumar²⁴ introduced the ideas based on continuum theory, while Viehland, Mason, and collaborators first formulated solutions of the Boltzmann equation for ions²⁵ that later became the basis for solving the electron problem.²⁶ Another important development in the time-dependent theory for charged-particle swarms driven by electro-magnetic fields was the work of Ness.²⁷ Kinetic Monte Carlo and particle-ion-cell (PIC) modelling approaches gained popularity as they are able to solve complex problems with much fewer assumptions or approximations and lend themselves more easily to the inclusion of complex boundaries. The roots of the PIC method go back to the self-consistent calculations of Buneman²⁸ and Dawson.²⁹ Much of the theory was later worked out by Hockney and Eastwood³⁰ and Birdsall and Langdon³¹ and discussed in more recent reviews.^23,32

The behavior and characteristics of non-equilibrium plasmas are strongly influenced by ionization,³³ especially at high E/N values. In this context then, the appropriate distribution of energies between the primary and secondary particles becomes important. The first reported analysis on energy partition following impact ionization was by Tzeng and Kunhardt.³⁴ However, comparisons with available experimental data were not attempted. The other important point is that values of the ionization coefficient (α) have not been calculated at high E/N levels and carefully compared with experiments. Here, we perform Monte Carlo simulations to probe the field-dependent behavior of the ionization coefficient for a wide range of E/N values (including fields as high as 4000 Td), with a focus on the energy partition between the primary and secondary electrons. Nitrogen gas at atmospheric pressure is used as a simple illustrative example, and our results are carefully compared with available experimental data.

The theory of the ionization of atoms by electron impact was developed by Wannier,³⁵ who treated the motion of escaping electrons on the basis of classical mechanics. Later, Vinkalns and Gailitis,^36,37 as well as Rau,³⁸ used Wannier's technique to investigate the angular distributions of outgoing electrons. It was shown that, due to electron repulsion, the two electrons escape in nearly opposite directions from the residual positively charged core to minimize system energy. The angular distribution was close to a Gaussian (with a mean of π–radians) and did not depend much on the energy, while the width of the angular distribution was a monotonically decreasing function of the energy.³⁹ Here, the ideas discussed in the above approaches to impact ionization of neutral molecules by energetic electrons are used for a Monte Carlo implementation, with details given below.

Instantaneous collision for an incident particle of energy E_inc is assumed, leading to the formation of a scattered primary particle of energy E′ and a secondary electron of energy E_sec. Thus, E′+ E_sec = E_inc − E_ion, where E_ion represents the ionization threshold energy. Furthermore, without any loss in generality, the incident electron is chosen to have a velocity along the z-direction. If Ψ represents an angle between the scattered primary electron and the secondary electron, one has v′·v_sec/[|v′| |v_sec|] = cos(Ψ), where v′ and v_sec are the velocity vectors of the scattered primary and secondary electrons. In addition, conservation of momentum also holds and leads to the relations: v_xsec = −v_x′ and v_ysec = −v_y′. As a result, one obtains the following relation for the angle Ψ:

Finally, if the partition of energies between E′ and E_sec is given by the ratio R, then E′= (E_inc − E_ion) (1−R) and E_sec = (E_inc − E_ion) R. This partition ratio can be chosen as a variable. In keeping with the results by Gailitis³⁷ that showed the energy partitioning function to be almost constant, R was taken to be a random number in the range [0,1]. Support for this near uniform energy partition has also come more recently from quantum-mechanical calculations by Shakeshaft,⁴⁰ though the close agreement with results from classical dynamics is also remarkable.^41,42

The distribution f(Ψ) for the scattered angle Ψ between the primary and secondary electrons [in Eq. (1)] was chosen to be Gaussian as suggested in the literature.^38,39 More specifically, f(Ψ) = K exp[−(Ψ−π)²/(2σ²)], for σ = 3.5 E^¼ with E = E_inc − E_ion denoting the energy in eV.³⁸ For numerical Monte Carlo implementations, the angle Ψ can be obtained by the standard technique of generating two random numbers r₁ and r₂ each in the interval [0,1]. Thus, Ψ  =  π + σ sin[2π r₁ (−2 Ln{r₂})^½]. Finally, based on momentum conservation, one has

which leads to

Given E′ and E_sec, one can solve Eq. (3) for v_z′, which can then be substituted into Eq. (1) to yield the v_zsec value. Finally, upon picking a value for v_x′ ≤ [2E′/m − v_z′²]^½, i.e., v_x′ = [2E′/m − v_z′²]^½r₁ based on a random number r₁ ε [0,1], the value of v_y′ is realized as: v_y′ = [2E′/m − v_z′² − v_x′²]^½. Velocities for the secondary electron are obtained from v_xsec = −v_x′ and v_ysec = −v_y.

Here, the Monte Carlo approach has been used for computations of the electron ionization coefficient. Details of the numerical implementation have already been reported and discussed by our group elsewhere.^43,44 The standard scattering mechanisms for both elastic and inelastic collisions between electrons and nitrogen were followed,^23,45 and the cross-sections were taken from the study by Itikawa et al.⁴⁶ Finally, the first Townsend coefficient α was calculated from the MC simulations as the ratio of the ionization rate S(E) and the electron drift velocity (v_dr), averaged over the entire ensemble of particles at each time step for different energy partitions. In order to probe the role of energy portioning, the following three cases were analyzed: (i) assignment of zero energy to the secondary particle, (ii) a random distribution of energies between the primary and secondary particles, and (iii) equal energy sharing between the primary and secondary electrons.

Monte Carlo simulations were carried out for nitrogen gas at atmospheric pressure to determine the field-dependence of the first Townsend coefficient α. The results were first obtained for the mean electron energy as a function of the applied external electric field. Figure 1 shows the mean energy with the coefficient R taken to be a uniformly distributed random number in the [0,1] range. Given the interest here in field-assisted impact ionization, relatively high values of the field F were chosen. The dependence is roughly parabolic as may be expected from a crude calculation based on the equation of motion: mdv/dt = qF − mv/τ_m, where m denotes the electron mass, v the velocity, and τ_m the energy-dependent momentum relaxation time. This leads to the mean kinetic energy: ⟨mv²/2⟩ ∼ E² ⟨τ_m⟩².

The results for the electron drift velocity were also obtained and are given next, again for randomly distributed partition coefficient R. Figure 2 shows the field dependence of the drift velocity. Also shown in the figure are the experimental data points reported in the literature for nitrogen. The values were taken from the studies by Hasegawa et al.,⁴⁷ Wedding and co-workers,⁴⁸ Roznerski,⁴⁹ Fletcher and Reid,⁵⁰ Yousfi et al.,⁵¹ and Blevin and Kelly.⁵² The Monte Carlo results seem to be in very good agreement with the data, at least over the lower E/N values. Measurements have not been performed at very high fields, and it is possible that simulation values predicted at large field strengths could differ from measured data. This could arise, for example, because experiments typically probe bulk drift velocities, while Monte Carlo simulations fold in both the drift and diffusive components. The latter can be non-negligible if strong impact ionization creates successively larger electron densities, leading to a diffusive flux that opposes field driven drift. More details on this phenomenon have recently been discussed by Vass et al.⁵³ and Dujko and co-workers.⁵⁴ In any event, Fig. 2 does suggest that the numerical scheme incorporating secondary electron generation, based on a uniformly randomized distribution of available energy between the primary and secondary electrons, roughly produces results in line with experiments. Although satisfying, values of the electron impact ionization coefficient and their comparison with available data at high electric fields constitute the real test of our formulation and the dependence on energy partition.

So finally, field-dependent values of the first Townsend coefficient were obtained from the Monte Carlo simulations and directly compared to experiments. The focus was on determining the role of energy partitioning, and the following different assignments were chosen: (a) the secondary particle was designated zero energy, (ii) a random distribution of energies between the primary and secondary particles, and (iii) equal energy sharing between the primary and secondary electrons. In the simulations, electric field values ranged to values as high as 4000 Td, starting from 1000 Td. The Monte Carlo results are shown in Fig. 3, along with comparisons to reported experimental data. The older data by Felsenthal and Proud⁵⁵ are seen to fall somewhat on the higher side. However, reports by Hays et al.,⁵⁶ Haydon and Williams,⁵⁷ Legler,⁵⁸ Yousfi,⁵¹ and Pitchford and Phelps⁵⁹ are all close together. The Legler and Pitchford data points match the simulations especially well. The Felsenthal and Proud result, which is much older, perhaps deserves some comment. Figure 8 of their report⁵⁵ for nitrogen was used, and paired values for the electric field F and Pτ were obtained. Here, P denotes the gas pressure and τ the formation time.⁵⁶ Expression for the Pτ product is provided in their Eq. (7), which was used to extract the ionization coefficient α. The numerator in their equation was taken to be “Ln(10⁸)” as indicated in their Appendix I, and the drift velocity was also provided by Felsenthal and Proud. However, there may be several issues with their approach. For example, the numerator in each case, which is at different points along the E/N curve, is really not known. There could also be some reservations regarding their field-dependent drift velocities. Also, electron transport properties are especially difficult to quantify experimentally at high fields because they depend on initial and experimental conditions, and issues arising from runaway particles may not be precisely determined.

Our predicted values with random energy partitioning also do compare well with the plot obtained by Hays et al.⁵⁶ Figure 3 shows that with secondary particles assigned zero starting energy, the field-dependent ionization values were predicted to be the lowest. Using equal energy sharing, on the other hand, yielded higher ionization values. These outcomes are all qualitatively expected based on the inherent physics. For example, a zero energy assignment to a secondary particle simply increases its “dead-time” before it can cross the ionization threshold and become capable of producing an ionizing event. Although the primary electron has more energy in this partitioning scheme, it is also more prone to stronger scattering given the monotonic dependence of this process on energy. The net outcome is a reduction in the ability of such a primary-secondary pair to strongly fuel impact ionization events. Equal sharing of excess energy can be expected to have just the opposite effect and produce the most conducive environment for ionization, a possibility supported by the results of Fig. 3. The random portioning scheme produces results between the two scenarios discussed above. It may be mentioned that partitioning based on the differential cross-sections had been proposed in the literature,³¹ with data from Opal, Peterson, and Beaty (OPB) being used for the differential cross-sections.⁶⁰ However, the randomization is probably a better option for Monte Carlo calculations given its inherent simplicity and the convenience of not requiring additional differential cross-section data.

For completeness and convenience, the simulation result of Fig. 3 was curve-fit to yield an analytic expression for the field-dependent electron impact ionization coefficient. The relation α(F) = 1.4764 × 10⁴ exp(−3 × 10⁵/F) resulted, with the field F in units of V/cm and ionization coefficient α given in units of cm⁻¹. The mathematical form α(F) = K₁ exp(−K₂/F^β) was chosen based on discussions reported by Meek and Craggs.⁶¹ This field dependence is very similar to the phenomenological relationship by Shockley, Wolff, and Baraff originally discussed in the context of impact ionization in semiconductors. For low-to-negligible phonon scattering (as would be the case in gases), β = 1 and so the form α(F) = K₁ exp(−K₂/F) results.⁶² For completeness, it may be mentioned that at very high fields in semiconductors where phonon and inter-valley scattering dominates, the value of β is 2,⁶³ while at intermediate fields, one has 1 < β < 2.⁶⁴ 

In conclusion, the behavior and characteristics of non-equilibrium plasmas are known to be strongly influenced by electron ionization, especially at high E/N values. Both primary and secondary electrons have to be considered, and in this context, the appropriate distribution of energies between the two particles becomes important. Here, several energy partitioning schemes, based on Wannier's classical ionization treatment, were used to probe the evolution of electron swarms in response to electric fields that ranged up to 4000 Td. Our Monte Carlo simulation results for the field-dependent drift velocities matched the experimental data available in the literature. Most of the reported experimental values are at the lower fields. More interestingly, the first Townsend coefficient predicted by the present Monte Carlo procedure was shown to be in very close agreement with reported data for E/N values ranging up to 4000 Td. Of the different energy portioning schemes tested, the random distribution of excess energies between the primary and secondary particles yielded the best comparisons and hence provides a simple approach to such simulations.

This work is based upon the work supported by the Air Force Office of Scientific Research under Award No. FA95501010106. One of us (RPJ) would like to thank K. Ness (Intertek Robotic Laboratories, Australia) for useful discussions.