Three-dimensional first principles simulation of a hydrogen discharge

Townsend discharge theory is commonly used to describe and approximate the ionisation fraction growth rate in the very early phase of plasma initiation in tokamak devices via ohmic breakdown. The prediction of the ionisation fraction growth rate is done most commonly with continuum or kinetic models, which in turn boil down to the relation between the first Townsend’s coefficient α, pressure p and electric field strength E (namely, α/p and E/p). To date there are few computational models that attempt to simulate the ionisation fraction growth rate via explicit modelling of each ionisation event through electron-neutral collisions. This is largely due to the challenge of addressing the exponential growth of charged particles from ionisation processes, combined with the high computational cost of N-body simulation. In this work, a new fully three-dimensional, first-principles model of a Townsend hydrogen discharge is demonstrated and benchmarked against prior experimental findings. These tests also include comparisons of three separate models for the scattering angle and their impact on the obtained α/p and mean electron drift velocity. It is found that isotropic scattering combined with restricting the freed electron’s scattering angle along the incident electron’s velocity vector during ionisation events gives the closest agreement of α/p compared to experimental measurements.


Introduction
In the very early stages of plasma initiation in a tokamak via ohmic breakdown, the main mechanism by which the * ITER is the Nuclear Facility INB no. 174. The views and opinions expressed herein do not necessarily reflect those of the ITER Organization. This publication is provided for scientific purposes only. Its contents should not be considered as commitments from the ITER Organization as a nuclear operator in the frame of the licensing process.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 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 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 superparticle 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 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 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-ofconcept 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.

Townsend avalanche theory
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 was proposed to describe the amplification factor of the current I that flows across the electrodes compared to the photoelectric current I 0 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 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 m −1 Pa −1 and B = 93.6 V m −1 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 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 2 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 Details of the numerical set up will be covered later in section 5.

Particle motion and collision models
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 2 collisions.

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 with m representing the mass of the charged particle, v the charged particle's velocity, q denoting its charge, E and 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 ) .
The updated velocity v t+∆t/2 is then obtained by treating the linear acceleration and angular rotation separately. Denoting v − and v + as the velocity vectors before and after rotation due to magnetic field, provides the updated velocity for the subsequent particle position update for time t + ∆t, given as 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 with ω g being the gyro-frequency, q e denotes the charge of electron and m 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 In this work, the external electric field strength E ext would be a prescribed constant throughout the simulation. The internal electric field E 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.

Electron-hydrogen cross sections
The cross sections of electron-H 2 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 + H 2 → e + H 2 (elastic scattering) • e + H 2 → e + H * 2 (vibrational and rotational level transition) 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 → 1 and J = 0 → 2 are also included, considering the sizeable cross sections in the lower energy range between 10 −1 and 10 1 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 + H 2 → e + H + 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 2 molecules. At the end of a typical simulation run, the estimated number density of H is at least 10 11 times smaller than for H 2 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 O(1 µs), thus it is assumed that the ionisation fraction is negligible.
The cross section for total scattering σ T is included in figure 1, as it plays a role in the computation of electron-H 2 collision probabilities which will be presented in section 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: where σ c (ε i ) denotes the cross section for collision outcome of type c which is dependent on electron's energy ε i , and n(x i ) describes the hydrogen number density at electron's position x i . The distance travelled by any electron is obtained by multiplying the time interval of ∆t with its velocity magnitude ∥v i ∥. In the current implementation, it is assumed that any given successful ionisation event is sufficiently rare such that n(x i ) can be taken constant over time. It is also assumed that the hydrogen molecules are uniformly distributed such that n(x i ) = n. The collision frequency is then 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, ν ′ , 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 ν ′ can be accepted as long as it fulfils ν T (ε i ) describes the collision frequency computed from summing all considered cross sections of electron-hydrogen molecule collisions, Since it is established that ν ′ should always be equal or larger than ν T (ε i ), the resulting collision probability is then P ′ denotes the collision probability calculated with ν ′ and P i refers to the actual collision probabilities computed from the energy carried by electron i, using the total scattering cross section σ T (ε i ). In this set up, a 'null collision' outcome with no real interaction between hydrogen molecule and electron is already included in ν ′ . This would mean that only when a random number A second random number R 2 is then used to determine the actual collision outcome. This implies that choosing an overly large ν ′ 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 ν ′ is optimised as 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 ⩾ 2) per time step as In order to have less than 1% chance for each simulated particle to have more than one collision per time step (r < 0.01), then ∆t must necessarily be chosen to fulfil

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.

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 where R 3 and R 4 are two random numbers in the range of [0, 1]. With the above θ and ϕ, the unit vector of the velocity after collision then is 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 → 2 and 0.516 eV for vibrational level transition of v = 0 → 1. Assuming that the incident electron's velocity magnitude is denoted by v init , the net velocity vector after scattering then is 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 nondissociative 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]).v inc refers to the initial unit vector of the incident electron, whilev ′ 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.
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 Thus an assumption that the newly created ion does not acquire any significant portion of the incident electron's energy (∥v 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 v ′ freed and incident v ′ inc electrons are fully described with Equation (14) is solved for every instance where ionisation occurs, it can be further simplified since m freed = m 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 v ′ inc and v ′ freed would be forward scattered relative to the incident electron's unit vector during ionisation events.

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]: where R 3 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 4 . 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. 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.

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 ,

ξ(ε) changes depending on different interaction potentials and is deduced from
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 ξ.
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.

Parallel N-body avalanche model
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 × 10 5 cm −3 (approximately 1.8 × 10 6 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 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.
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.
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.

Simulation setup
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 2 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).
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 −12 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 −1 Torr −1 and pressure at the range of 0.5 < p < 10.5 Torr. One set of reported parameters is E/p = 156 V cm −1 Torr −1 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 E ext which all simulated charged particles will experience. In all the presented simulation cases, E 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 + 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.
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 −1 Torr −1 , note that the recorded q an is mostly from electrons. It can be observed that the recorded charges settle to a steady fluctuation about a constant mean value.
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 with q b.s. describing the number of back scattered electrons recorded at the cathode. It should be noted that the spread of the recorded q 0 reduces as E/p increases. Higher E/p can be achieved by either increasing E ext or reducing gas pressure p.
Having larger E 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. The mean charges at the anode Q an and the mean net photoemitted electrons Q 0 are computed (along with their standard deviations), followed by the first Townsend coefficient α obtained through The uncertainty of α is then calculated through propagation of errors.

Results
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 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. Recall that the choice of time step size has an upper bound given in equation (10). This implies that a pre-filled pressure of 2.5 Torr (0.33 kPa) would require a maximum time step size of O(1 ps) in order to keep the probability of more than one collision per time step per particle below 1%. Influence of electric field strength on time step size is also considered. However, the corresponding limit is much larger than the suggested picosecond scale.
In order to choose the proper time step size for the simulation, a set of tests is done with the initial choice of 1 ps, 0.4 ps, 0.25 ps and finally 0.1 ps. The convergence tests are done at E/p = 400 V cm −1 Torr −1 with pressure of 2.5 Torr, using the presented random scatter model. A comparison of the obtained α/p for each of the chosen time step sizes is done and the result is tabulated in table 1.
The tabulated results shows that the obtained mean values are within 1 standard deviation of each other. This suggested that the choices of time steps sizes are all acceptable. However, the variation of obtained mean α/p values as the step sizes get smaller do show a trend of convergence. In light of that, the lower time step size choices of 0.25 and 0.1 ps are preferred. For the following simulation cases, 0.25 ps is chosen to help reduce the total time required to complete the simulations.

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 −1 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 −1 Torr −1 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 −1 Torr −1 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 b.s. in equation (18), thus giving a narrower spread of q 0 as E/p increases. Similar trend of reduction in the uncertainty of recorded q 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 −1 Torr −1 , this could be attributed to the neglected hydrogen dissociations that produce hydrogen atoms, Figure 12. Comparison of our results for the mean electron drift velocity using all three models to reported values from: Simulation data obtained from Saelee et al [31], derived data from Schlumbohm [32] and Lisovskiy et al [33], measured data from Blevin et al [34], and Roznerski et al [35]. which in turn provide another channel for additional impact ionisation.
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 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 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 DE from presented scattering angle models and the various prior works [31][32][33][34][35]. It can be observed that the models for calculating the scattering angle by Vahedi and Okhrimovskyy give good  figure 4. The random scatter model gives an overall lower estimation of the mean electron drift velocity V 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 ⩽ 200 V cm −1 Torr −1 , but it is seen that the agreement improves at 400 and 600 V cm −1 Torr −1 . The comparison between the three presented models shows that both energy dependent scattering angle calculations give better approximation of V DE .

Discussion
In our simulations, electrodes are only defined as boundaries where charge collection occurs. Intricate interaction between charged particles and electrodes taking place in actual experiments, such as charged particles sheath formation in close proximity to electrodes or secondary electron creation through ion collision on cathodes, are not captured. With that in mind, the current study chooses to benchmark the proposed model with experimental work focusing on measuring the first Townsend coefficient α. As such, the presented set of benchmark simulations are performed with the shortest separation distance d derived from Rose's work as well as using the lower ranges of operating pressure.
A number of modifications to the numerical setup are required to extend the pressure and electric field range of validity for α measurements in electrode plate experiments. Chief among them is the modelling of the electric field via potential difference over a defined separation distance rather than a prescribed constant E ext . Electric potential changes across the electrodes at longer separation distance would then be modelled correctly.
However, the current model is more suitable in the simulation of tokamak ohmic breakdown since the electric field experienced by charged particles are induced by the current ramp up in the central solenoid rather than potential difference between electrodes. Establishing a method to predict the rise of electron population via electron-neutral collisions, as demonstrated in this work, is one of the components required to capture the very early phases of tokamak plasma initiation where the electron number density is well below 10 3 m −3 . Since the obtained α presented in this work has good agreement with reported values in high E/p ranges, it fits the ohmic breakdown phase of tokamak devices such as NSTX at E/p ≈ 500 V cm −1 Torr −1 [1] or ITER at a predicted E/p ≈ 400 V cm −1 Torr −1 [18].
On a separate note, Townsend's work [17] reported different values of current amplification with experiments at the same E/p ratio with differing operating pressure. This trend is noted in the simulation as well, but not explicitly reported in studies conducted by different authors [36]. More work could be done to study the effect on obtained α at different pressure settings while maintaining a constant E/p ratio as well as constant electrode separation d.

Conclusion
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 DE . In order to improve the approximation in both α/p and V 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 −3 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.