Three-dimensional first principles simulation of a hydrogen discharge

In the very early stages of plasma initiation in a tokamak via ohmic breakdown, the main mechanism by which the charged particle population increases can be described by the Townsend discharge process [1]. This is triggered when certain conditions are met, i.e. the combination of electric field strength and gas pressure is such that free electron population growth becomes exponential with respect to their trajectory length. Tokamak breakdown theory is often a simplified one-dimensional (1D) description of a Townsend discharge in an actual tokamak setting. As such, this causes an overestimation of electrons' collisional drag force since they only travel in one direction, rather than a prescribed scatter angle distribution in the direction perpendicular to principle external field. A 1D model would be unsuitable for approximation of runaway electrons as well, as it is heavily dependent on electrons' interaction with complex three-dimensional (3D) field geometry. For these reasons, there is a need for fully 3D models that can simulate Townsend avalanches in the presence of a magnetic field. This work represents the first step in an effort to create such a model, where we first benchmark the code using the well-established discharge process in an electrode plate experiment. Ultimately our goal is to then apply the model to tokamak devices, expanding on the works of 0D model simulation by Jiang et al [2] and 2D simulation by Yoo et al [3]. In order to simulate the discharge process, several key components such as numerical representation of charged particles, determination of electron-neutral collision likelihoods and outcomes, scattering angle calculations for electrons experiencing collision events and numerical method to handle the space-charge effects needs to be considered.

To date there have been numerous computational models proposed for the simulation of Townsend discharges. Early attempts in the 1960s by Ward [4] (which was further examined by Börsch-Supan et al [5]) involve computing the growth of electron currents at the cathode by solving the differential form of Gauss's law while assuming a 1D setting for the discharge tube. This method approximates the measured current over time via an iterative finite difference method with a given initial first Townsend's coefficient α. As such, there is no explicit representation of electrons nor the need to consider collision events. During the same period, early works of Monte Carlo collision (MCC) calculations were proposed [6, 7]. Electrons were explicitly tracked and the cross sections considered. However, the simulated domain remained 1D and the tracked number of electrons was in the order of thousands. A sufficient sample size required to approximate α and mean electron drift velocity $V_\textrm{DE}$ was obtained through repetition of the simulation. Due to the small population of charged particles, space-charge effects were not considered. More recent works of MCC methods [8] are capable of handling a larger charged particle population through the description of super-particles. However, loss of details in the actual collision frequency of individual charged particle is unavoidable since a super-particle has difficulties in tracking the evolving energy distribution of the particles it represents over time.

Selection of the electron-neutral cross sections also requires careful consideration. Balance between computation efficiency and accuracy of the approximation is important for cases which involve a high number density of simulated particles. Benchmarking a computational model with existing experimental results provides a quantitative measure of approximation accuracy, which is an important step to make informed choices in regard to which cross sections are important.

In the application of MCC in Townsend discharge simulation, the choice of scattering angle calculation influences the obtained α and $V_\textrm{DE}$, as well as the velocity distribution of electrons (aside from the electric field and neutral gas pressure), which in turn affects the collision frequency of electrons. Surendra et al described the scattering angle based on a look-alike form of screened Coulomb scattering [9], while Okhrimovskyy et al derived an expression via the first Born approximation scattering theory [10]. Comparison between the models and previously established experiments is essential to understand how the different nature of anisotropy in the resulting scattering angle affects the ionisation rate.

As the charged particle number continues to rise during an avalanche, Coulomb forces will become more apparent and eventually need to be accounted for. One of the established methods to compute the internal fields in discharge simulations is the combined particle-in-cell and Monte Carlo collision (PIC-MCC) method [9, 11]. Another notable method in the simulation of Townsend discharge is a hybrid model of Monte Carlo simulation and Boltzmann equation solver [12], where the the collision outcome is determined by collision operators and the internal electric field is computed by solving the continuity equation as a function of position. Both the described PIC method and hybrid model work well with simple simulated domain geometries such as cylinders, since the created regular grid can easily represent the physical boundaries. It is less computationally efficient, however, to create the same grid structure for a toroidal geometry. There will be a large number of cells where particles should not travel, and additional work is required to ensure the grid structure respects the physical boundary. Thus, a mesh free method which circumvents the meshing process is presented in our work and it is capable of solving the internal field regardless of the simulated geometry.

While there are numerous computational models and methods for discharge between electrodes in simulations, their application with hydrogen as gaseous medium is lacking. Ward's work focused mainly on the breakdown in air [13], while Tkachev et al applied the PIC method in combination with known α and analytical solutions to predict electron growth in helium [14, 15]. Vahedi et al applied the method of PIC-MCC to simulate discharges in oxygen and argon [8]. A more recent work by Jovanović et al applies the MCC method in the simulation of methane discharge [16]. Additional work in simulating hydrogen discharge between electrodes is required to quantitatively measure the accuracy of the obtained α and $V_\textrm{DE}$, before applying such a model to tokamak plasma initiation scenarios.

This work presents an alternative combination of methods to simulate the Townsend discharge and derive the first Townsend coefficient (α) in hydrogen. Firstly, it is a proof-of-concept of a complete discharge simulation with all charged particles (electrons and positive ions that arise from ionisation events) simulated and tracked in three dimensions in both space and momentum. The internal electric field is then computed by employing a parallel tree algorithm to compute the electrostatic field, which is helpful as it is less constrained by the geometry the charged particles are occupying. Additionally, a new scattering angle model is developed for the MCC method and compared to the works by Vahedi et al and Okhrimovskyy et al. The accuracy in reproducing the obtained α in low pressure conditions and high E/p range will be beneficial in the simulation of the ionisation growth rate of very early phases of tokamak ohmic breakdown, which falls in the Townsend discharge regime.

The paper's structure is as follows: We start with the basics of Townsend avalanche theory described in section 2. This is then followed by detailed descriptions of the particle motion as well as the implemented collision models in section 3. The framework of our parallel tree code for resolution of charged particle potentials is briefly discussed in section 4, which is followed by a detailed description of the numerical experiment setup presented in this work in section 5. Results and discussions that follow are then presented in section 6.

Experimental work performed by J S Townsend on the growth of the electrical current between parallel electrode plates due to impact ionisation formed the basis of Townsend discharge theory. His work focused on finding the relationship between the exponential rise in free electrons through ionisation and the influence of variables such as neutral gas pressure and applied electric field [17]. As a result of impact ionisation, the expression

$$\begin{equation} \frac{I}{I_0} = \frac{e^{\alpha d}}{1 - \gamma (e^{\alpha d} - 1)} \end{equation} \tag{ 1 }$$

was proposed to describe the amplification factor of the current I that flows across the electrodes compared to the photoelectric current I₀ at the cathode. Here, d denotes the separation between the two electrodes, α(E, p) defines the number of ionisations per unit length and γ is a coefficient to include secondary electron emissions due to positive ion impact at the cathode.

An analytical expression for α can be obtained by considering the electrons that have reached the ionisation energy, their corresponding mean free path at a given operating pressure as well as the probability of collision with the help of known cross sections. This gives the expression of

$$\begin{equation} \alpha = Ap \,\exp \left(-\frac{Bp}{E}\right), \end{equation} \tag{ 2 }$$

with A and B constants which are dependent on the ionisation energy of the neutral target and its cross sections. For the work presented here, the neutral target is purely hydrogen molecules, so $A = 3.83~\textrm{m}^{-1}\,\textrm{Pa}^{-1}$ and $B = 93.6\ \textrm{V m}^{-1}\, \textrm{Pa}^{-1}$ would normally be used to predict the current growth [18].

Electrons in a neutral gaseous medium will gain energy due to the external electric field and over time potentially cause ionisation of the neutrals during collisions. The collisions between electrons and neutral targets will in turn prevent an arbitrary rise in electron kinetic energy, resulting in a mean drift velocity ($V_\textrm{DE}$) in parallel with the accelerating field. The acceleration experienced by charged particles scales proportionally with the magnitude of the external electric field, hence it is expected that the mean drift velocity scales with E/p.

Considering plasma initiation via ohmic breakdown in a tokamak setting, ionisation is mostly due to electron—H₂ impacts. Secondary ionisation mechanisms through positive ions impacting vessel boundaries are not considered. In this work, an experimental setup that involves a small electrode separation d and operates below breakdown voltage is deliberately chosen for comparison with the numerical simulations. The γ parameter is negligible under such conditions [19], thus simplifying equation (1) to

$$\begin{equation} \frac{I}{I_0} = e^{\alpha d}. \end{equation} \tag{ 3 }$$

Details of the numerical set up will be covered later in section 5.

This section provides the description of the implemented equation of motion integrator, as well as the various components that concern the events of electron-H₂ collisions.

3.1. Velocity update scheme

In order to prepare for the application of this code in scenarios with magnetic fields, a numerical integrator that has a suitable bound on the energy conservation error when both electric and magnetic fields are present has to be selected. The integrator should also avoid unbounded error accumulation over an arbitrarily large number of time steps. The Boris scheme [20] is implemented for the current work. Although it does not meet the criteria of being a symplectic integrator, it does conserve phase space volume as well as exhibiting a global bounded energy error [21].

The equation of motion of charged particles under the influence of both electric and magnetic fields is given by

$$\begin{equation} m \frac{\textrm{d}\mathbf{v}}{\textrm{d}t} = q(\mathbf{E} + \mathbf{v} \times \mathbf{B}) \end{equation} \tag{ 4 }$$

with m representing the mass of the charged particle, $\mathbf{v}$ the charged particle's velocity, q denoting its charge, $\mathbf{E}$ and $\mathbf{B}$ each representing the experienced electric and magnetic field, respectively. Numerical evaluation of the acceleration is performed by first discretising equation (4) in time t with central differences, which gives

$$\begin{equation*} \frac{\mathbf{v}_{t + \frac{\Delta t}{2}} - \mathbf{v}_{t - \frac{\Delta t}{2}}}{\Delta t} = \frac{q}{m} \left(\mathbf{E} + \frac{\mathbf{v}_{t + \frac{\Delta t}{2}} + \mathbf{v}_{t - \frac{\Delta t}{2}}}{2} \times \mathbf{B}\right). \end{equation*}$$

The updated velocity $\mathbf{v}_{t + \Delta t/2}$ is then obtained by treating the linear acceleration and angular rotation separately. Denoting $\mathbf{v}^-$ and $\mathbf{v}^+$ as the velocity vectors before and after rotation due to magnetic field,

$$\begin{equation*} \begin{split} \mathbf{v}^{-} & = \mathbf{v}_{t - \frac{\Delta t}{2}} + \frac{q\mathbf{E}}{m} \frac{\Delta t}{2}\\ \frac{\mathbf{v}^+ - \mathbf{v}^-}{\Delta t} & = \frac{q}{2{m}} (\mathbf{v}^+ + \mathbf{v}^-) \times \mathbf{B}\\ \mathbf{v}_{t + \frac{\Delta t}{2}} & = \mathbf{v}^{+} + \frac{q\mathbf{E}}{m} \frac{\Delta t}{2} \end{split} \end{equation*}$$

provides the updated velocity for the subsequent particle position update for time t + Δt, given as

$$\begin{equation} \textbf{x}(t + \Delta t) = \textbf{x}(t) + \textbf{v}_{t + \frac{\Delta t}{2}} \Delta t. \end{equation} \tag{ 5 }$$

There is a time step size limit associated with the Boris integration scheme in case the gyro-motion of the fast moving electrons are to be resolved. The limit is expressed as

$$\begin{equation*} \Delta t \lt \frac{2}{\{\omega_\textrm{g}\}_{\textrm{max}}} = \frac{2m_\textrm{e}}{q_\textrm{e} |B|_{\textrm{max}}}, \end{equation*}$$

with $\omega_\textrm{g}$ being the gyro-frequency, $q_\textrm{e}$ denotes the charge of electron and $m_\textrm{e}$ describes the electron rest mass. However, this limit is not considered in our work since the magnetic field is absent. The time step size is instead bounded by the collision probability in the Monte Carlo method, more details will be discussed in section 6.1. The influence of electrostatic potentials of other charged particles is also considered, such that the electric field vector is

$$\begin{equation} \textbf{E} = \textbf{E}(t) = \textbf{E}_\textrm{ext} + \textbf{E}_\textrm{int}(t). \end{equation} \tag{ 6 }$$

In this work, the external electric field strength $\textbf{E}_\textrm{ext}$ would be a prescribed constant throughout the simulation. The internal electric field $\textbf{E}_\textrm{int}(t)$ due to the spatial distribution of charged particles is solved at every time step through the use of the PEPC solver, a parallel tree code implementation of the Barnes-Hut algorithm [22]—see section 4.

3.2. Electron-hydrogen cross sections

The cross sections of electron-H₂ molecule collisions used in the current work are based mainly on the reported values from Yoon et al [23]. They are a compilation of decades of experimental work and recommended cross sectional data for an electron energy ranging from 0.001 to 1000 eV as shown in figure 1. For the current work, the following reactions are included for the simulated collision outcomes:

• $e + \textrm{H}_2 \rightarrow e + \textrm{H}_2$ (elastic scattering)
• $e + \textrm{H}_2 \rightarrow e + \textrm{H}^*_2$ (vibrational and rotational level transition)
• $e + \textrm{H}_2 \rightarrow 2e + \textrm{H}^+_2$ (non-dissociative ionisation)
• $e + \textrm{H}_2 \rightarrow 2e + \textrm{H} + \textrm{H}^+$ (dissociative ionisation)

Zoom In Zoom Out Reset image size

The included collision cross sections are chosen to focus primarily on the ionisation of zero-point ground state of hydrogen molecules, where its rovibrational level v and J is both at 0. The cross sections for rovibrational level transition for $v = 0 \rightarrow 1$ and $J = 0 \rightarrow 2$ are also included, considering the sizeable cross sections in the lower energy range between 10⁻¹ and 10¹ eV. Electrons that are involved in collisions of this manner will have their energy subtracted via equation (13), specifically Δε of 44.1 meV for rotational level transition and 0.516 eV for vibrational level transition of hydrogen molecule. It should be noted that the prefilled neutral gas molecules are not explicitly simulated in general, in an effort to save computational resources.

Note that the $e + \textrm{H}_2 \rightarrow e + \textrm{H} + \textrm{H}$ cross section is absent in this work. A deliberate choice is made not to include this, since its overall contribution to α via ionisation of hydrogen atoms is negligible due to its low number density compared to H₂ molecules. At the end of a typical simulation run, the estimated number density of H is at least 10¹¹ times smaller than for H₂ at 273.15 K and 333.31 Pa. The hydrogen molecule number density is also assumed to be constant throughout the simulation. This assumption is made since the simulated currents reach an equilibrium at a time scale of $\mathcal{O}(1\,\mu$s), thus it is assumed that the ionisation fraction is negligible.

The cross section for total scattering $\sigma_\textrm{T}$ is included in figure 1, as it plays a role in the computation of electron—H₂ collision probabilities which will be presented in section 3.3.

3.3. Collision probability

As noted earlier, the target hydrogen molecules are treated as a background species, only the electrons and the subsequently created ions are individually simulated. Since the background species is not treated explicitly, a collision event is then modelled by a collision probability based on the average number density of hydrogen molecules, the electron velocity, the time step size, and collision cross sections. This results in the following expression for the collision probability of the ith electron with a hydrogen molecule:

$$\begin{equation} P_{c,i} = 1 - \exp\left(-\sigma_c(\varepsilon_i) \,n(\textbf{x}_i) \,\|\textbf{v}_i\| \Delta t \right), \end{equation} \tag{ 7 }$$

where $\sigma_c(\varepsilon_i)$ denotes the cross section for collision outcome of type c which is dependent on electron's energy ε_i , and $n(\textbf{x}_i)$ describes the hydrogen number density at electron's position $\textbf{x}_i$. The distance travelled by any electron is obtained by multiplying the time interval of Δt with its velocity magnitude $\|\textbf{v}_i\|$. In the current implementation, it is assumed that any given successful ionisation event is sufficiently rare such that $n(\textbf{x}_i)$ can be taken constant over time. It is also assumed that the hydrogen molecules are uniformly distributed such that $n(\textbf{x}_i) = n$. The collision frequency is then

$$\begin{equation} \nu_c(\varepsilon_i) = \sigma_c(\varepsilon_i) \,n \,\|\textbf{v}_i\|, \end{equation} \tag{ 8 }$$

which is solely a function of electron energy.

In order to avoid potentially costly computation of all individual collision frequencies ν_c for every particle at every time step, Skullerud proposed the null collision method which involves defining a constant collision frequency, $\nu^{^{\prime}}$, that is large enough to encompass all possible collision outcomes over the full range of electron energy [24] in a first step. As such, any measure or definition of $\nu^{^{\prime}}$ can be accepted as long as it fulfils

$$\begin{equation*} \nu^{^{\prime}} \geqslant \nu_\textrm{T}(\varepsilon_i),\ \forall \varepsilon_i. \end{equation*}$$

$\nu_\textrm{T}(\varepsilon_i)$ describes the collision frequency computed from summing all considered cross sections of electron-hydrogen molecule collisions,

$$\begin{equation*} \nu_\textrm{T}(\varepsilon_i) = n \|\textbf{v}_i\| \sum_c \sigma_c(\varepsilon_i) = n \|\textbf{v}_i\| \sigma_\textrm{T}(\varepsilon_i). \end{equation*}$$

Since it is established that $\nu^{^{\prime}}$ should always be equal or larger than $\nu_\textrm{T}(\varepsilon_i)$, the resulting collision probability is then

$$\begin{equation} \begin{split} 1 - \exp(-\nu^{^{\prime}} \Delta t) &\geqslant 1 - \exp(-\nu_\textrm{T}(\varepsilon_i) \Delta t)\\ P^{^{\prime}} &\geqslant P_i, \ \forall i. \end{split} \end{equation} \tag{ 9 }$$

$P^{^{\prime}}$ denotes the collision probability calculated with $\nu^{^{\prime}}$ and P_i refers to the actual collision probabilities computed from the energy carried by electron i, using the total scattering cross section $\sigma_\textrm{T}(\varepsilon_i)$. In this set up, a 'null collision' outcome with no real interaction between hydrogen molecule and electron is already included in $\nu^{^{\prime}}$. This would mean that only when a random number R₁ = [0, 1] fulfills $R_1 \leqslant P^{^{\prime}}$, proper computation of the each available $\sigma_c(\varepsilon_i)$ is done. A second random number R₂ is then used to determine the actual collision outcome.

This implies that choosing an overly large $\nu^{^{\prime}}$ will not impose a penalty to the computational cost, it would be the same as not applying this simplification in the first place. However, if the value of $\nu^{^{\prime}}$ is optimised as

$$\begin{equation*} \nu^{^{\prime}} = n\ \max_\varepsilon(\sigma_T(\varepsilon) \|\textbf{v}\|), \end{equation*}$$

the resulting algorithm will be most efficient.

The choice of time step size is important to ensure that only one instance of collision occurs for every charged particle in each simulation step. Vahedi et al [8] gave an expression for the probability r of two or more collisions ($k \geqslant 2$) per time step as

$$\begin{equation*} r \sim \sum^\infty_{k = 2} P^k_{i} = \frac{P^2_{i}}{1 - P_{i}}. \end{equation*}$$

In order to have less than 1% chance for each simulated particle to have more than one collision per time step ($r \lt 0.01$), then Δt must necessarily be chosen to fulfil

$$\begin{equation} \sigma_\textrm{T}(\varepsilon_i) \,n \,\|\textbf{v}_i\| \Delta t \leqslant 0.1. \end{equation} \tag{ 10 }$$

3.4. Scattering angle

In a quasi-neutral plasma, Coulomb collisions dominate in the collisional drag force experienced by electrons. However, in the early stages of Townsend discharge, presence of ions is not prominent within the medium. Most of the collision events experienced by the electrons are between electron and neutral targets and thus deserve special attention. Since the differential cross sections related to electron-hydrogen elastic scattering over the wide energy range of interest is lacking, several methods of computing the scatter angle are considered in this work. This includes a new method of scatter-angle calculation that is compared to two others by Vahedi et al [8] and Okhrimovskyy et al [10], respectively.

The interaction between simulated charged particles is handled by a tree code algorithm called Pretty Efficient Parallel Coulomb (PEPC) solver in its force summation subroutine (as we will describe in section 4). As such, a transition from a elastic collision dominated phase to a primarily Coulomb collision phase in the ionisation breakdown simulation can be handled seamlessly.

3.4.1. Random scatter model.

A new method is presented here to resolve the scattering angle of electrons, as well as the resulting energy and momentum partition during collision events. Two distinct treatments of the post collision scattering angle are implemented, specifically for non-ionising and ionising collision events. In the case of non-ionising collisions, the electron scattering angles for elastic scattering and vibrational-rotational level transition collision outcomes are determined in a random, isotropic manner. This is chosen since the neutral targets are not explicitly simulated, thus actual scattering angles cannot be computed. In this work, φ is defined as the polar angle with a range of [0, π] and θ is the azimuthal angle ranged [0, 2π]. In order to ensure that the velocity vector is uniformly distributed over a sphere, θ and φ are calculated as

$$\begin{equation} \begin{split} \phi & = \arccos(1 - 2 R_3),\\ \theta & = 2 \pi R_4 \end{split} \end{equation} \tag{ 11 }$$

where R₃ and R₄ are two random numbers in the range of [0, 1]. With the above θ and φ, the unit vector of the velocity after collision then is

$$\begin{equation} \hat{\mathbf{v}} = \begin{bmatrix} \sin(\phi)\cos(\theta)\\ \sin(\phi)\sin(\theta)\\ \cos(\phi)\end{bmatrix}. \end{equation} \tag{ 12 }$$

Once the randomised unit vector is obtained, the incident energy of electrons before impact will be conserved after the collision for the elastic scattering case. In the case of rovibrational level transitions, the transition energy Δε has to be deducted from the initial energy carried by the electrons before collision. The values amount to 0.0441 eV for rotational level transition of $J = 0 \rightarrow 2$ and 0.516 eV for vibrational level transition of $v = 0 \rightarrow 1$. Assuming that the incident electron's velocity magnitude is denoted by $v_\textrm{init}$, the net velocity vector after scattering then is

$$\begin{equation} \mathbf{v}_\textrm{final} = \sqrt{v_\textrm{init}^2 - \frac{2 \Delta \varepsilon}{m_e}} \cdot \hat{\mathbf{v}}. \end{equation} \tag{ 13 }$$

The treatment of velocities after an ionisation event is markedly different from elastic collisions. During an ionisation event, the ionisation energy is first deducted from incident energy via equation (13)), with Δε = 15.426 eV for non-dissociative ionisation and Δε = 18.075 eV for dissociative ionisation. The remaining energy is then split between the incident and newly freed electrons by satisfying conservation of both momentum and energy. To achieve this, we start with the freed electron's scattering direction determined randomly, but with the restriction that the angle of scattering is always over a hemisphere which is perpendicular to the incident electron vector shown in figure 2 (φ = [0, π/2]). $\hat{\mathbf{v}}_\textrm{inc}$ refers to the initial unit vector of the incident electron, while $\hat{\mathbf{v}}^{^{\prime}}_\textrm{freed}$ denotes the unit vector of the freed electron after ionisation. The incident electron's direction after the impact is then calculated in order to fulfil the conservation conditions.

Zoom In Zoom Out Reset image size

With the hydrogen molecule about three orders of magnitude heavier than an electron, the resulting energy transferred from electron to hydrogen molecule is expressed by

$$\begin{equation*} E_{\textrm{kin,}\textrm{H}_2} = \frac{4m_\textrm{e}m_{\textrm{H}_2}}{(m_\textrm{e} + m_{\textrm{H}_2})^2}E_{\textrm{kin,e}}. \end{equation*}$$

Thus an assumption that the newly created ion does not acquire any significant portion of the incident electron's energy ($\|\mathbf{v}_\textrm{ion}\| = 0$) can be made. The set of equations to be solved in order to obtain the incident electron's post collision velocity reduces to a two-body elastic collision problem.

As such, the assignment of velocities for freed $\mathbf{v}^{^{\prime}}_\textrm{freed}$ and incident $\mathbf{v}^{^{\prime}}_\textrm{inc}$ electrons are fully described with

$$\begin{equation*} \begin{split} \mathbf{v}^{^{\prime}}_\textrm{inc} & = \|\mathbf{v}^{^{\prime}}_\textrm{inc}\|\hat{\mathbf{v}}^{^{\prime}}_\textrm{inc} \\ \mathbf{v}^{^{\prime}}_\textrm{freed} & = \|\mathbf{v}^{^{\prime}}_\textrm{freed}\|\hat{\mathbf{v}}^{^{\prime}}_\textrm{freed}, \end{split} \end{equation*}$$

where

$$\begin{equation} \begin{split} \cos \theta^{^{\prime}} & = \frac{\mathbf{v}^{^{\prime}}_\textrm{freed} \cdot \mathbf{v}_\textrm{inc}}{\|\mathbf{v}^{^{\prime}}_\textrm{freed}\|~\|\mathbf{v}_\textrm{inc}\|}\\ \|\mathbf{v}^{^{\prime}}_\textrm{freed}\| & = \frac{2m_\textrm{inc}}{m_\textrm{freed} + m_\textrm{inc}} \|\mathbf{v}_\textrm{inc}\| \cos \theta^{^{\prime}}\\ \mathbf{v}^{^{\prime}}_\textrm{inc} & = \mathbf{v}_\textrm{inc} - \frac{m_\textrm{freed}}{m_\textrm{inc}}\mathbf{v}^{^{\prime}}_\textrm{freed}. \end {split} \end{equation} \tag{ 14 }$$

Equation (14) is solved for every instance where ionisation occurs, it can be further simplified since $m_\textrm{freed} = m_\textrm{inc}$. Note that the hydrogen atoms created through dissociative ionisation are not explicitly simulated in favour of reduced computational cost. The implication of equation (14) in combination with the forward scattering restriction is that both $\mathbf{v}^{^{\prime}}_\textrm{inc}$ and $\mathbf{v}^{^{\prime}}_\textrm{freed}$ would be forward scattered relative to the incident electron's unit vector during ionisation events.

3.4.2. Vahedi et al model.

In this case, calculation of the loss of energy experienced by the incident electron is still done by equation (13). However, the scattering angle experienced by the incident electron is now dependent on its energy [8]:

$$\begin{equation} \cos \phi = \frac{2 + \varepsilon - 2(1 + \varepsilon)^{R_3}}{\varepsilon} \end{equation} \tag{ 15 }$$

where R₃ denotes a random number of the range [0, 1] and ε refers to the incident electron's energy. Another random number is used to determine the azimuthal angle θ = 2πR₄. We show the dependence of the scattering angle on the random number for a choice of electron energies in figure 3. It is clear that for this scattering model, the electrons are more likely to have scattering angles that are less than π/2 (forward scatter) as ε increases. At an energy of 0.01 eV, there is an almost equal chance of forward or backward scattering.

Zoom In Zoom Out Reset image size

The scattering angle distribution in equation (15) is also applied to the incident electrons during ionisation, while the freed electrons' post-collision velocities are resolved with equation (14). As the electrons' scattering angle is not restricted to the forward hemisphere, there will be electrons that move in the opposite direction to its progenitor after ionisation occurred. Thus, the drag force experienced by the electrons is higher when compared to the random scatter model. Note also, that in Vahedi's original work only the energy is conserved while momentum assignment is randomised.

3.4.3. Okhrimovskyy et al model.

Another scattering model for the calculation of scattering angle φ of electrons follows the method proposed by Okhrimovskyy et al [10]. In his work, one finds that the scattering angle is represented by

$$\begin{equation} \cos \phi = 1 - \frac{2R_3(1-\xi)}{1 + \xi(1-2R_3)}, \qquad R_3 = [0,1]. \end{equation} \tag{ 16 }$$

ξ(ε) changes depending on different interaction potentials and is deduced from

$$\begin{equation} \frac{\sigma_m}{\sigma} = \frac{1- \xi}{2 \xi^2} \Big((1+\xi) \ln \frac{1+\xi}{1-\xi} - 2\xi \Big), \end{equation} \tag{ 17 }$$

where σ_m and σ denote momentum transfer cross section and integral cross section of electron hydrogen collisions.

The result of ξ in turn affects the resulting scattering angle distribution (obtained by solving equation (16)) shown in figure 4, where the scattering angle is more inclined to back scatter the more negative the value for ξ becomes. The opposite is true for positive values of ξ.

Zoom In Zoom Out Reset image size

Note that the resulting scattering angle distribution in figure 4 is markedly different from figure 3. This method produces enhancements in both forward and backward scattering limits depending on the incident electron's energy. For electron energies below 3.8 eV, there is a higher likelihood of backscatter, peaking at approximately 1 eV. Once the electron energy is above the threshold of 4 eV, it quickly becomes highly likely to scatter forward instead, as seen in the case of 8 and 12 eV, respectively. This scattering angle distribution is also applied to both ionising and non-ionising collisions.

While the simulation of hydrogen discharge does not involve magnetic fields, there is strong motivation to employ the use of a mesh free potential solver as a preparation for the latter application in tokamak settings. One of the drawbacks of meshed algorithms like PIC codes is the issue of additional numerical diffusion, observed in anisotropic problems such as magnetised plasmas [25]. This condition is triggered when the mesh is not aligned to the magnetic field's principle direction. A mesh free N-body solver does not suffer from the same difficulty and is thus better suited to tackle simulations that involve complex field geometry. Another merit of the mesh free method is the flexibility it provides in capturing the domain geometry. There is no actual mesh boundary for particle positions, thus avoiding the tedious process of dynamically re-meshing the domain as the simulation continues.

The mesh free solver we use is PEPC, the PEPC-solver developed at the Jülich Supercomputing Centre [22]. Depending on the simulation setup of the Townsend discharge simulation, the number density of simulated particles can reach up to $7.6 \times 10^{5} \textrm{cm}^{-3}$ (approximately 1.8 × 10⁶ particles). PEPC provides the framework to compute the resulting collective Coulomb forces. It adopts an implementation strategy proposed by Warren and Salmon [26] called hashed oct tree algorithm, using a Hilbert space filling curve rather than the original Z order curve for mapping particle coordinates to unique keys. This is done to avoid complicated communication between processors and heavy duplication of local tree data during the tree construction and traversal phase of Barnes–Hut algorithm. As a result, the computational complexity of an N-body problem is now reduced to the order of $\mathcal{O}(n \log n)$ in PEPC. The structure of PEPC is divided into three parts [27].

The tree construction and traversal part of the Barnes–Hut algorithm are built into the tree code routine, shown in the middle of figure 5. This part of PEPC contains most of the parallel communication required for the tree-code, as well as the handling of load-balancing between processors. Various utilities such as initialisation of the MPI environment and parallel input/output operations are also placed here. The interaction specific modules serve as the 'backend' for PEPC, which handle the particle data structure as well as providing subroutines that compute multipole moments of grouped particles [28]. It also includes subroutines and functions for the force calculation during tree traversals. The final part to PEPC is the 'frontend' which includes the implementation of components discussed in section 3. Since the framework for the Barnes–Hut algorithm is mostly coded in the other two parts, this makes development of numerical studies for n-body physics relatively self-contained in the frontend itself.

Zoom In Zoom Out Reset image size

The overall computational workflow performed in the simulation is shown in figure 6. The aforementioned cross sections are first tabulated so that the collision frequency for each electron can be computed during the time step loop. The operating pressure of the experiment is also given as an input during initialisation phase to determine the background neutral number density n, and subsequently the corresponding collision frequency in equation (8). The internal electric field is computed via the tree code algorithm at each time step. The external electric field is then added to the internal field for each simulated particle before performing the velocity update via the Boris scheme described in section 3.1.

Zoom In Zoom Out Reset image size

At the start of every time stepping loop, the list that contains all simulated particles is sorted and filtered for any particle that has come into contact with the electrodes. The remaining particles are then accelerated and translated before collision checks are performed. During the collision checks, the outcome for every simulated electron is then calculated. Any newly created charged particle is added to the particle list so that their influence on the internal electric field can be evaluated in the preparation for the new time step with equation (6). The loop repeats itself until the prescribed total loop is reached.

The experimental setup with two parallel electrodes driven by a constant electric field and the absence of a magnetic field [17, 19, 29] are used as benchmark simulations that measure the α values in hydrogen gas medium. As explained in section 2, this study focuses on discharge via the first ionisation mechanism (electron—H₂ impact ionisation), as investigated experimentally by Rose [19]. The numerical setup is therefore designed to mimic the experimental geometry and the numerically obtained α are compared against his reported values. There is a continuous stream of free electrons that are photoemitted by directing ultraviolet radiation towards the cathode. These freed electrons are then accelerated in between the electrodes and the resulting current is measured. Since the electrode separation d is defined as one of the setup parameters, the first Townsend coefficient α can be obtained via equation (3).

Zoom In Zoom Out Reset image size

In our model setup, the position of the cathode is assumed to be at z = 0, while the anode is simply a numerical offset d in the −z direction, prescribed at the start of the simulation. A preset number of electrons to be 'added' at z =−0.01 for every time step is used to emulate the photoemitted current from the cathode due to ultraviolet radiation. Tagashira et al [29] reported a photoemitted current at the order of 5 × 10⁻¹² A, due to the nature of time discretisation associated with numerical simulation, it would mean that the number of photoemitted electrons per time step at a resolution of picoseconds is less than one. Instead, a series of convergence tests found that a prescribed 20 electrons per time step is sufficient to reach statistical stability for the obtained mean value of first Townsend coefficient α. Their x and y coordinates are resolved by a process of rejection sampling to ensure a uniform spatial distribution of electron source from the cathode. The diameter D of the electrodes is set to 5 cm, similar to the experimental setup used by Rose [19]. Considering the fact that there will be at least four random numbers required for each simulated electron for every time step, a pseudo random number generator that has sufficiently long period and is highly parallelisable has to be chosen. As such, every instance of random number used in this work is generated through 'Random123', a counter-based random number generator [30].

Determination of the electrode separation d used in the simulation is derived from reported values of Rose's experimental work, with E/p value ranging from 15 to 1000 V cm⁻¹ Torr⁻¹ and pressure at the range of $0.5\lt p\lt10.5$ Torr. One set of reported parameters is E/p = 156 V cm⁻¹ Torr⁻¹ with voltage across electrodes measured at 50 V and the pressure of 2.7 Torr. As such, the electrode separation d worked out to be approximately 0.12 cm. All simulations presented in this work are performed at a temperature of 0^∘ which then determines the neutral number density. The electric field created by potential difference between separated electrodes are simulated by prescribing a constant $\mathbf{E}_\textrm{ext}$ which all simulated charged particles will experience. In all the presented simulation cases, $\mathbf{E}_\textrm{ext}$ is parallel to the +z direction.

Each of the added electrons near the cathode is assumed to carry an initial energy of 1 eV, travelling in the −z direction. Any particles that travelled to either above z = 0 or below z =−d are removed from the simulation, while charges 'lost' to the electrodes in each time step are recorded. Figure 8 shows the typical spatial distribution of the simulated electrons during a simulation run. Figure 8(a) shows that there are lighter spots distributed over the electrode region, they are areas which electrons are sparse. It is worth mentioning that H⁺ and H$_\textrm{2}^+$ ions are not shown in this figure to avoid cluttering. Of interest in figure 8(b) are the electrons that are scattered away from the electrode regions as they near the anode. This can be attributed to the electron-neutral elastic scattering and Coulomb scattering effects experienced by the electrons as they travel between the electrodes. These electrons will not be counted toward the measure of α value.

Zoom In Zoom Out Reset image size

At every time step, sum of all the charges carried by charged particles (electrons, H⁺ and H$_2^{+}$) that reached the anode is recorded. Figure 9 shows the recorded sum of the charges (measured with the unit of 1e) at the anode under E/p of 400 V cm⁻¹ Torr⁻¹, note that the recorded $q_\textrm{an}$ is mostly from electrons. It can be observed that the recorded charges settle to a steady fluctuation about a constant mean value.

Zoom In Zoom Out Reset image size

It was shown in section 3.4 that the implemented scattering angle models allow back scattering. As such, rather than the expected 20 photoemitted electrons at every time step as described earlier, the measured net photoemitted electrons q₀ also fluctuates as shown in figure 10. The expression for q₀ at every time step is

$$\begin{equation} q_0 = 20 - q_{\textrm{b.s.}}, \end{equation} \tag{ 18 }$$

with $q_{\textrm{b.s.}}$ describing the number of back scattered electrons recorded at the cathode. It should be noted that the spread of the recorded q₀ reduces as E/p increases. Higher E/p can be achieved by either increasing $\mathbf{E}_\textrm{ext}$ or reducing gas pressure p. Having larger $\mathbf{E}_\textrm{ext}$ will cause larger downwards acceleration (away from the cathode at z = 0) for electrons. As such, electrons are less likely to reach the cathode after back scattering events. On the other hand, lowering the pressure p effectively reduces the collision frequency between electrons and neutrals, thus allowing longer electron mean-free-path.

Zoom In Zoom Out Reset image size

The mean charges at the anode $Q_\textrm{an}$ and the mean net photoemitted electrons Q₀ are computed (along with their standard deviations), followed by the first Townsend coefficient α obtained through

$$\begin{equation} \alpha = \frac{1}{d} \ln(I/I_0) = \frac{1}{d} \ln(Q_\textrm{an}/Q_0). \end{equation} \tag{ 19 }$$

The uncertainty of α is then calculated through propagation of errors.

The explicit representation of electrons enables a detailed study of their spatial distribution, as well as the influence of cross sections and scattering angle calculations on the energy at any chosen time step. Figure 11 shows the electron number density and averaged energy over distance for E/p = 400 V cm⁻¹ Torr⁻¹.

Zoom In Zoom Out Reset image size

An interesting observation can be made regarding the high electron number density near the cathode, which diminishes rather abruptly over a short distance. This is explained by the lack of low energy electrons' mobility as they are initially introduced into the system. As they travel along the accelerating electric field, a spread in the energy distribution occurs due to the probability of collisions with the neutral background (determined by collision frequency ν) and the distribution of the resulting scattering angle. This spread affects the spatial distribution of electrons since they do not move at the same velocity.

It is noted that the electron energy gain over distance is smooth up to approximately 15 eV, which can be attributed to the lack of major energy sink within that energy range. Referring to the cross section plot in figure 1, the only other reactions that can act as energy sink in the energy range below 15 eV are vibrational and rotational level transitions. However, their impact on electron energy is small as the rotational level transition Δε is 44.1 meV and the cross section for vibrational level transition is 1 order of magnitude smaller than elastic scattering at best.

Once the electron's energy threshold goes above the requirement for ionisation, additional cross sections are considered during collision events. The number density and the average electron energy continues to rise as the travel distance increases. There is a small but noticeable drop in number density beginning from a distance of approximately 0.11 cm as well as a rise in average electron energy. This can be attributed to slightly diminishing cross sections above 40 eV, meaning that the electrons encounter less collisions thus more likely to move unimpeded.

6.1. Influence of Δt on obtained α

6.2. Influence of E/p on obtained α

Rose performed the experiment with various combinations of electric field strengths E and pre-filled gas pressures p. In order to compare with those results, appropriate values for electric field strength and pressure are chosen for the numerical simulations. Specifically, the chosen E/p are 60, 100, 200, 400 and 600 V cm⁻¹ Torr−1. All the cases were benchmarked at a plate separation of d = 0.12 cm. Referring to table 2, the simulation severely underestimates the mean α/p at the E/p of 60 V cm⁻¹ Torr⁻¹ for all three implemented scattering angle models. However, good agreement with Rose's reported α/p is observed in the range of 100–400 V cm⁻¹ Torr⁻¹ for the newly proposed random scatter model, while the models by Vahedi and Okhrimovskyy result in lower mean α/p values. Another interesting observation is that the uncertainty of α/p reduces as E/p increases. This is due to diminishing value of $q_\textrm{b.s.}$ in equation (18), thus giving a narrower spread of q₀ as E/p increases. Similar trend of reduction in the uncertainty of recorded $q_\textrm{an}$ is observed as well. The obtained mean α/p are also noticeably lower compared to Rose's reported value at E/p of 600 V cm⁻¹ Torr⁻¹, this could be attributed to the neglected hydrogen dissociations that produce hydrogen atoms, which in turn provide another channel for additional impact ionisation.

Table 2. Comparison of obtained mean α/p values.

	Mean α/p [cm−1 Torr−1]
E/p	Rose	Random	Vahedi	Okhrimovskyy
60	0.50	0.10 ± 1.19	0.09 ± 1.06	0.08 ± 1.30
100	1.26	1.24 ± 0.89	1.07 ± 0.82	0.97 ± 0.97
200	2.55	2.57 ± 0.62	2.18 ± 0.59	2.00 ± 0.67
400	3.60	3.56 ± 0.51	2.94 ± 0.52	2.71 ± 0.57
600	4.10	3.85 ± 0.40	3.14 ± 0.44	2.95 ± 0.47

Note, however, that both Vahedi's and Okhrimovskyy's models were proposed to simulate scattering with larger atoms or molecules (for example, Ar and O$_\textrm{2}$). It is also important to note that Okhrimovskyy's original work did not account for ionisation cases. Subjecting freed electrons to the anisotropic scattering angle distribution in either figures 3 or 4 causes exaggeration of the electrons' energy loss when compared to the forward scattering restriction implemented in the random scatter model. As a result, there are fewer electrons that reach the energy threshold required for further ionisation, which is reflected in the diminished α/p value. This suggests that Vahedi's and Okhrimovskyy's models require different methods to account for energy partitions during ionisation events, other than the classical treatment of energy and momentum conservation presented at equation (14).

In the current work, the mean electron drift velocity $V_\textrm{DE}$ is obtained by taking the last simulated time step of the electron population, followed by averaging all the electrons' velocities. Figure 12 shows the comparison of measured $V_\textrm{DE}$ from presented scattering angle models and the various prior works [31–35]. It can be observed that the models for calculating the scattering angle by Vahedi and Okhrimovskyy give good agreement to other previously reported results. Note that the measured drift velocity using Okhrimovskyy's model gives a slightly lower mean result when compared to Vahedi's model at E/p = 60 V cm⁻¹ Torr⁻¹ whereas it exhibits a higher value above 100 V cm⁻¹ Torr⁻¹. This can be attributed to the extremes exhibited in figure 4.

Zoom In Zoom Out Reset image size

The random scatter model gives an overall lower estimation of the mean electron drift velocity $V_\textrm{DE}$ when compared to either Vahedi's or Okhrimovskyy's model. It is most likely due to the isotropic distribution of scattering angle independent of electron energy as introduced by the random scatter model. This means that there is a higher fraction of the electron population that will be back scattered compared to both Vahedi's and Okhrimovskyy's model, thus introducing more collisional drag to electrons. The random scatter model underestimates reported values of previous works at $E/p \leqslant$ 200 V cm⁻¹ Torr⁻¹, but it is seen that the agreement improves at 400 and 600 V cm⁻¹ Torr⁻¹. The comparison between the three presented models shows that both energy dependent scattering angle calculations give better approximation of $V_\textrm{DE}$.

6.3. Discussion

A 3D first principle simulation of Townsend discharge is presented and benchmarked against prior experimental work. It is found that the proposed random scatter model produces good agreement in obtained α/p values when compared to Rose's reported values. The produced current amplification is a good approximation in situations where the primary means of ionisation is through electron-neutral collisions. In order to properly correct the obtained α/p in settings with higher pressure, additional work would be needed to account for secondary ionisation by positive ions as well as charged particle interaction with electrodes.

Different methods to partition energy and momentum among newly ionised electrons and ions are required, should the scattering angle model of Vahedi et al or Okhrimovskyy et al be implemented. The classical treatment of energy and momentum conservation presented in equation (14) coupled with the back-scattering during ionisation events cause excessive drag force experienced by electrons when the aforementioned models are used, thus underestimating the obtained α/p values. Our work shows that the choice of implemented scattering angle model plays a major role in the simulated Townsend discharge in a full 3D application. It is also observed that while Vahedi's and Okhrimovskyy's models underestimate α/p, they provide a better approximation of the resulting mean electron drift velocities $V_\textrm{DE}$. In order to improve the approximation in both α/p and $V_\textrm{DE}$, a combination of forward scatter restriction presented in section 3.4.1 during ionisation events, with equations (15) or (16) in other scenarios could be studied further.

The proposed random scatter model shows promising signs that it can be used to model the ionisation fraction growth rate during the plasma initiation phase of tokamak ohmic breakdown, especially in low pressure conditions (approximately 10⁻³ Pa) where the main mechanism of ionisation is through Townsend discharge. Future work will focus on the model's application in a realistic 3D toroidal geometry along with the presence of a confining magnetic field configuration.

The authors gratefully acknowledge the computing time granted through JARA, project JZAM04, on the supercomputer JURECA at Forschungszentrum Jülich.