Development of a novel DC micro-discharge code on dsmcFoam+for plasma eﬀect in DSMC ﬂow proﬁle simulations and its application to a micro-plasma system

Recent advances in the miniaturisation of plasma-based microelectromechanical systems have permitted their application in the space sector, such as ﬂow control of ﬂying vehicles and the potential usage in high precision attitude and orbit control systems for CubeSats. To determine its contribution in pressure, velocity, and temperature proﬁles for thrust analysis on micro-ﬂuidic channels in space/vacuum environments, a novel DC micro-discharge code has been developed in the dsmcFoam+ solver in order to simulate the ionic eﬀects on the gas ﬂow during the breakdown voltage in microgaps. The rareﬁed gas ﬂow was solved using the direct simulation Monte Carlo (DSMC) method, with the Fowler–Nordheim equation to calculate the ionic current density during the vacuum ﬁeld emissions eﬀect for the micro-plasma simulation. The incorporation of the Fowler–Nordheim model avoids the addition of an exponential increase in computational power requirements for a DSMC simulation. Instead of introducing new electrons and ions as equivalent particles into the system, the existing DSMC gas particles in the cells are accelerated. Through this, a computationally lightweight alternative for micro-scale plasma simulations util-ising the dsmcFoam+ solver is achieved, which can be tailored to the thrust requirements of a particular mission.


Introduction
In recent years, micro-electromechanical systems (MEMS) have seen exponential growth in their employment as flow control devices and in on-board attitude and orbit determination systems (AOCS) for CubSats (Wie et al., 2014). Applications of plasma micro-scale AOCS thrusters yield nano-Newton to micro-Newton scale thrust, with potential applications on space applications in the near future (Zhuang et al., 2014;Krejci et al., 2019;Demmons et al., 2022;Levchenko et al., 2018). Rarefied gas dynamics must be taken into account since the thrusters are firing into the vacuum of space. Therefore, the Boltzmann equation is required to provide a solution for the simulation of cases where molecular collision rates reduce, making the continuum fluid assumption in the Navier Stokes equations invalid due to the degree of rarefaction (Ebrahimi and Roohi, 2017).
Flow rarefaction is defined by Knudsen numbers Kn > 10 in a free molecular regime, where inter-molecule collisions are not an important factor in the transfer of mass, momentum, and energy (White et al., 2018;Capon et al., 2017). Therefore, rarefied flow simulations based on real scenarios involve the transition Knudsen number regime (10 < Kn < 0.1), where a solution to the Boltzmann equation must be sought, but inter-molecular collisions are still an important process, which makes this equation difficult and numerically expensive to solve. Methods such as direct simulation Monte Carlo (DSMC), particle-in-cell (PIC) and test particle Monte Carlo (TPMC; Capon et al., 2017) are widely used to obtain solutions to problems in the transition Knudsen number regime.
A significant problem is the inclusion of plasma in simulations due to the method's limitations and the high processing power requirements, which is a numerically expensive situation. An electron-gas and ion-gas particle interaction model needs to be simplified in order to replicate plasma phenomena in a more computationally efficient manner. There are two key capabilities that are not considered by the majority of the available methods; one is a phenomenon known as the quantum effect on microgaps, based on an increase in the electric field intensity by an enhancement factor when a voltage is applied to microgaps between electrodes. The second key is the macro view of the gas flow after electrons, gas particles, and ions' interactions. Herein lies the novelty of this work, which includes a method for simulating plasma that takes electric field enhancement and macroscopic flow analysis into account.
Particle-in-cell methods are used to model the movement/trajectory of charged particles within an electric field (Auweter-Kurtz et al., 2005) for plasma simulations. Unfortunately, to model the effect of plasma within the flow, the methods usually do not consider the interactions between the particles (Auweter- Kurtz et al., 2005) and the surface (Capon et al., 2017). Some solvers, such as DEMO-CRITUS (Lapenta, 2011), PicUp3D (Forest et al., 2006) and CPIC (Delzanno et al., 2013) are used to model the charging of satellites and instrument calibration (Capon et al., 2017).
DSMC simulations can reproduce the gas-surface interactions by modelling a single atom/molecule as a single particle that is representative of a large number of real gas atoms. Available DSMC solvers include dsmcFoam+ (White et al., 2018), MONACO (Padilla, 2010), DAC (Padilla, 2010), and SPARTA (Plimpton and Gallis, 2015). For applicability, dsmcFoam+ fulfils a similar role to other alternatives, but with the advantage of being an open-source (as is SPARTA) and implemented within OpenFOAM (Weller et al., 1998) (Open Source Field Operation And Manipulation), providing an open and extensible C++ based software package containing a wide range of libraries (White et al., 2018).
Combining the PIC capabilities to model the charged particles with the DSMC method for gas-surface interactions results in a hybrid PIC-DSMC method. Such PIC-DSMC codes are primarily focused on the research and design of plasma thrusters (Capon et al., 2017), such as PICLas (Auweter-Kurtz et al., 2005;Reschke et al., 2016), and DRACO (Brieda et al., 2004). Additionally, PIC-DSMC solvers are commonly used to model orbiting objects in ionised environment, such as low Earth orbit, an example of which is the pdFoam solver (Capon et al., 2017) based in OpenFOAM. Unfortunately, PICLas (Auweter-Kurtz et al., 2005) and DRACO (Brieda et al., 2004) do not focus on the micro-scale effects on plasma, avoiding the variations and enhancements of the electric field due to by the microgap between the electrodes. Furthermore, electron temperature (T e ) variation is essential in micro-plasma modelling to obtain an accurate electron density distribution; the Draco solver considers T e to be constant (Brieda et al., 2004). Additionally, PIC-DSMC methods for modelling micro-scale plasma thrusters results in an exponential increase in computational resources when they are used to properly model number density by the difference between the densities of the charged particles (electrons) and the gas/neutral particles. To illustrate the challenge, each simulator particle in pdFoam represents a large number of real gas atoms/molecules. Therefore, for an inlet flow between 1 kPa to 100 kPa, the particle density range of the gas is between 10 23 m À3 to 10 25 m À3 . An estimated range of electron density can be defined between 10 11 m À3 to 10 19 m À3 according to the electrical properties of the material, size and the applied current-voltage. PIC-DSMC methods require setting the gas and electron densities as initial and boundary conditions. It is recommended to have a minimum of twenty particles per cell to minimise the statistical scatter and ensure correct collision statistics when using dsmcFoam+ (Hassanvand et al., 2018). During the initialisation, some of the cells will contain just the electron density; consequently, with an electron density of 10 19 m À3 , each equivalent particle at least needs to represent 10 17 real atoms/molecules to fulfil the minimum of twenty particles requirement. Therefore, the inlet flows of 1 kPa and 100 kPa need 10 5 and 10 7 equivalent particles in each cell, assuming the highest electron particle density value. Thus, the more the electron density decreases, the number of equivalent particles exponentially increases up to a maximum of 10 15 particles per cell. With this assumption, modelling a case of 1000 cells will contain roughly a minimum of 10 8 equivalent particles, demanding extremely high computational resources. To get a point of comparison, White et al. (White et al., 2018) executed parallel tests with the dsmcFoam+ solver to determine its performance and scalability by running millions of cells with a total of 2 Â 10 7 ; 1 Â 10 8 ; 2 Â 10 8 , and 4 Â 10 8 equivalent particles using 960 cores on ARCHER-the UK's national supercomputing service. In contrast, running a case for 1000 cells on PIC-DSMC method would require up to 960 cores. PIC-DSMC quickly becomes an expensive and impractical method for the analysis of micro-thrusters.
Due to the flexibility and open-source nature of the solver, a dsmcFoam+ extension was developed to optimise the modelling of the DC micro-plasma effect on flow. The code combines the particle-surface interactions from dsmcFoam+ with the numerical solution of the Fowler-Nordheim equation to model the electrical charge interactions and ionisation rates, allowing for the simulation of micro-plasma effects.
Consequently, a micro-plasma systems with two inlets and a single outlet orifice were used to demonstrate the application of the DC micro-discharge code in the simulation and analysis of the system. Three electrode configurations were determined; (a) 1 anode/2 cathodes, (b) 2 anodes/2 cathodes, and (c) 1 anode extended/2 cathodes, and placed near the outlet for the purpose of analysing the different variations of flow profile caused by plasma.

Theory and Background
In the breakdown of a gas, usually termed Townsend's breakdown, there are two primary mechanisms that contribute to a significant rise in carriers-gaseous charge production through electron impact; the process of ionization the process of the cathode charge production through secondary emission (Go and Pohlman, 2010).

Townsend's first ionization coefficient, a
The coefficient a describes the number of electrons produced by a single electron moving 1 cm in the electric field direction during electron-neutral interactions (Kontaratos, 1965). The second electron produced by the first collision will repeat the same phenomena, impacting another neutral particle to create a third and subsequent accumulative electrons. Therefore, the number of electrons in the gap increases exponentially and this is known as electron avalanche. Through these collisions, ionization occurs at the moment the neutral particle specie releases an electron, changing its charge (Matejčik et al., 2015). Townsend's first ionization coefficient is expressed as where p is pressure, E is the electric field, and A and B are coefficients specific to the gas species. The rate of ion-electron pair production from a single electron after the subsequent avalanche on a distance x across the gap d is denoted n and determined as (Kontaratos, 1965;Wijsman, 1949), Consequently, the electron number density n e given by a current flowing in the gap after the avalanche effect is then given by (Meek and Craggs, 1953), where n e 0 is the initial electron number density.

Townsend's second ionization coefficient (c SEE )
Secondary electron emission c SEE is described as the ratio of the secondary electron flux emitted from the cathode surface to the influx of the positive ion bombardment (Meek and Craggs, 1953;Venkattraman and Alexeenko, 2012) c SEE e ad À 1 Depending on the application, c SEE can be treated as a constant (Go and Venkattraman, 2014). For a consstant discharge, the electron density produced by positive ions bombardment, considering secondary electron emission, is defined as (Meek and Craggs, 1953)

Field emission coefficient (c FE )
Unfortunately, the previously presented equations (which are based on the Paschen curve) improperly describe the breakdown voltage and do not take into account the enhancement effect of the electric field in microgaps. In this micron scale regime, electron field emission plays a significant role in the breakdown phenomena. The Fowler-Nordheim equation takes into consideration the ion-enhanced field emission given by the effects of the electrodes' surface conditions for the breakdown voltage. Consequnetly, Boyle and Kisliuk (Boyle and Kisliuk, 1955) took this approach and defined the field emission coefficient c FE as the ratio of field emission current incident to the ion current on the electric field E given by (Go and Venkattraman, 2014) where K is a constant, with an approximate value of 10 7 (Go and Pohlman, 2010), b is the field enhancement factor (typically for ideal polished metal surfaces this is equal to 50 (Moore et al., 2013;Go and Pohlman, 2010)), / is the work function of the electrode's material, and v y ð Þ is a corrective function. It is insufficient just to consider the coefficients a and c SEE to estimate the electron density for an applied current-voltage in a micrometre gap separation between the electrodes. Furthermore, the incorporation of the c FE coefficient in Eq. (5) is needed to consider the FEE, resulting in the following expression (Go and Venkattraman, 2014) This expression describes the rate of increase in the electron current density for an applied voltage in microgaps.

Ion density, n i
The ion density n i at the cathode is related to the total current/electron density in the gap (Venkattraman and Alexeenko, 2012). Each additional electron generated from the initial current comes with an ion, forming an electronion pair. Thus, the following relation was made Establishing a relation from Eq. (7) with Eq. (8), the ion density can also be expressed as 2.5. Corrective function v y ð Þ According to Miller (1966), v y ð Þ in Eq. (6) is a corrective function from the image force effect on the emitted electrons expressed by elliptical integrals in form of where v y ð Þ is expressed as (Miller, 1966) To facilitate the computation of the corrective function, a polynomial curve fitting equation were calculated for v y ð Þ. Reference values were taken from Miller's table I (Miller, 1966) between values of 0 to 1 of v y ð Þ from Eq. (11), with intervals of 0.01. The following expressions were obtained: v y ð Þ ¼ À0:7519y 2 À 0:2744y þ 1:0161: 2.6. Applied current density For a given applied current I the following relation needs to be taken into consideration (Godse and Bakshi, 2009): where n e 0 is the number of free charges per unit volume that flows through a volume segment Ad; q is the amount of charge on each carrier over a specific amount of time Dt.
For the development of the dsmcFoam+ code, Dt is the iterative time value that represents the time-step size.

Electric field
For a given applied voltage V in a microgap distance d, the electric field is calculated according to the number and location of the electrodes. In this work, it is important to highlight the consideration of plane electrode plates for all calculations of E.
For the DC micro-discharge code to execute and demonstrate its application (See Fig. 1), the electric field for the three different micro-plasma electrode configurations must also be defined. Therefore, electric field lines and values of E are needed for (a) 1 anode/2 cathodes, (b) 2 anodes/2 cathodes, and (c) 1 anode extended/2 cathodes configurations.
For parallel electrodes, one anode and one cathode, E is uniform along the x-direction of the gap, as shown in Fig. 2, caused by the same direction and average magnitude of the electric field lines. The electric field can be determined by the relation (Chabay and Sherwood, 2015) Eq. (14) describes an ideal case with two parallel plates. If the electrodes are set with a displacement between them affecting their relative orientation, or an additional electrode plate is include in the case, variations in the quantity, direction and magnitude of E at each point on the space are introduced, according to how close or far an anode segment is from the cathode. This can no longer be described by Eq. (14) and the modified electric field E n must be determined by the vector summation of the fields from the individual source charges q j (Chabay and Sherwood, 2015;Jahn, 2006), where Q is the total charge to the source given by the sum of each individual charge on the field, e 0 is the dielectric constant of the medium, r is the distance between the charges andr is the direction denoted as a unit vector. As illustrated in Figs. 3 and 4, having more than two electrode (anode-cathode) plates, as in the presented micro-plasma systems, creates an extra electric field. Where the charge Q represents the total negative charge from its source (cathode), denoted as Q cathoden to facilitate its identification. The DC micro-discharge code developed here accelerates the positively charged ions. As consequence, it is required to consider the calculation of the electric field from the point of view of the positive charges, where the ions are being attracted towards the cathode.
An additional cathode plate is introduced to the system at the same height as the first cathode and an anode is placed between them on the opposite side of the microgap, as illustrated in Fig. 3. In this Figure, the red line represents the vector sum of the electric fields E 1 and E 2 on a positive charge resulting in the total field E experienced at this specific location. Then, to calculate the components E i and E j of the total field E the following expressions can be used: Additionally, if two pairs of parallel anode-cathode plates are placed near each other with a separation in the x-axis, as shown in Fig. 4, each parallel electrode pair experiences two electric fields; (1) the uniform field due to its own flux of charges (Eq. (14)),and (2) the proximate field produced by its own electrode neighbour plates, described in Eqs. (16) and (17).
Thus, the total electric field experienced in a cell located inside the electrodes region can be obtained from

Ionic density distribution
A Maxwell-Boltzmann distribution is applied to determine the allocation and the specific electron-ion number of particles for each cell in the microgap region of the mesh.
The distribution of the microscopic thermal velocities of the electrons is calculated from the following function (Chen et al., 1984; Moisan and Pelletier, 2012): where f v ð Þ is the probability density for finding a particle with a velocity v; m e is the electron mass, k B is the Boltzmann constant, and T e is the local electron temperature.
The Maxwell-Boltzmann distribution curve is calculated for each cell column by integrating along the gap, as shown in Fig. 5, and is a function of the temperature and the velocity of the electron.
The calculation of the electron temperature in each column T e;i is possible once the average electric field E ij is obtained from the column of cells where the electron is located (Eq. (18)) by the relation between the electric potential energy U E and the average kinetic energy per particle K E . Therefore, where q is the charge of the particle exposed to an electric field E ij (Chabay and Sherwood, 2015). Additionally, the average kinetic energy per electron K E is given by (Chen et al., 1984), The change in the kinetic energy of the particle, in this case an electron, is equal to the change in its potential energy from the law of conservation of energy. Therefore the change in the potential of the electron in an electric field can be expressed as (Chabay and Sherwood, 2015), Solving for T e , the temperature of the electron is To solve Eq. (19), once the temperature of the electron is obtained, the average thermal speed of the Maxwell-Boltzmann distribution can be calculated by (Chen et al., 1984;Katz, 2015), Additionally, the most probable speed of the distribution can be obtained as (Chen et al., 1984;Katz, 2015)   Therefore, for the assignation of a velocity from the Maxwell-Boltzmann distribution along a column of cells, the width of the curve needs to be obtained and segmented by the number of cells in the column. For this, the probability distribution, Eq. (19), is a function of the thermal velocity v, consequently, as v is the independent variable, the width of the curve can be found by calculating the approximation of the maximum thermal velocity v e;i;max , given by, where v e;i;max represents the total width of the distribution curve, C is the confidence interval set as 0.998, an acceptable value allowing the width range to capture the distribution shape in the cell column. The closer C is to unity, it will cover more the width of the distribution curve tending to inifinity. Then the distribution width is segmented by the total number of cells n j;total in the column, as follows, Each cell n j in the column, corresponding thermal velocity v i;j , contains the cumulative sum of the segment velocities v i;seg from its predecessor cells, starting from the lowest to the highest cell position of the column, As result, the probability of the Maxwell-Boltzmann distribution for each cell is given by the function f v ij À Á from Eq. (19). In order to use this probability as a distribution factor to calculate the electron density in a cell, f v ij À Á must be normalised through division by the probability distribution function of the average thermal speed, f v e;i;ave ð Þ, obtained from Eqs. (19) and (24): The distribution factor d f v ð Þ ij has the purpose to take f v e;i;ave ð Þ as a factor 1.0, and initially coinciding it with the applied current density from Eq. (13). Thus, the average electron density of a column with the applied current density per unit volume is defined as, Therefore, to calculate the electron density after the breakdown voltage in a specific cell at i; j ð Þ, the distribution factor d f v ð Þ ij is incorporated in Eq. (7), resulting in the following expression: Consequently, by the relation of electron-ion pairs, described in Section 2.4, from Eq. (9), the ion density assigned to a specific cell is given by: For the DC micro-discharge code, n i;ij represents the density of gas particles to be accelerated in a cell located at i; j ð Þ in order to simulate the micro-plasma effect on the flow.

Particle acceleration
Once the ion density and its distribution are obtained, the energy of the positive ions must be calculated to determine the acceleration of gas particles to be applied in dsmcFoam+.
The ion's kinetic energy K E;i is obtained by (Kim et al., 2006), where E ij is the electric field experienced in the cell, n gas Dt ð Þ is the gas particle density as a function of the time step and r is the collision cross-section between the electrons and gas particles.
It is important to emphasise the discretisation of both particle densities n e 0 and n gas as a function of the timestep Dt for the algorithm purpose. The first variable is already discretised solving n e in Eq. (13), but for the case of n gas , it follows that where N A the Avogadro constant, M mass the molar mass, T gas the gas temperature, and q the gas density. Additionally, K E;i must consider the amount of energy transferred to a DSMC particle. Each DSMC particle represents a large number of real gas atoms/molecules, F N . If the ion particle density per unit volume n i;ij in the cell is less than F N ; K E;i needs to be multiplied by a transfer factor C P DSMC , the purpose of which is to assign the proportional amount of energy to a DSMC particle: Consequently, the ion energy K E;i assigned to a DSMC particle when n i;ij < F N , is defined as Therefore, if n i;ij < F N then just one DSMC particle in the cell at i; j ð Þ will receive the ion energy. Afterwards, the kinetic energy is converted into the average thermal Then the average thermal speed is converted to the acceleration a ! i þ on the DSMC particle by diving v i þ ;ave by the characteristic flow length, which is d, the gap between the electrodes (Su et al., 2017) If v i þ ;ave is a function of the electric field E i;j , consequently a ! i þ is a function of E i;j . Hence, the vector components x and y of a ! i þ are defined by the unit vectorr ij given from the distance r between the particle and its direction to the electric field obtained in Section 3.1.
Subsequently, the velocity vector u ! i;j is calculated by where Dt is the timestep set for each iteration. Finally, the DSMC particle is accelerated by adding the calculated u ! i;j velocity to the current velocity u ! 0;i;j , resulting in the following velocity, This is an iterative process, whereby u ! 0;i;j is updated in each cell when it is located between the defined electrodes area inside the mesh and experiences an ionization ratio E ij .

Leapfrog method
A leapfrog method is implemented in the code to avoid the non-physical increase of energy of the DSMC particles given by the numerical solutions during the acceleration. Since n i;ij ( F N is in the majority of the cases, only a small number of DSMC particles are accelerated. Therefore, Eq. 39 is split into two parts; half of the velocity is added to u !t p;i;j before the movement part of the DSMC algorithm, then the second half after the collisions with a velocity u !tþDt=2 p;i;j (Szabo, 2001). At the end of the timestep, the particle will be moving at u !tþDt p;i;j (Arber et al., 2015;Verboncoeur, 2005),

Time step selection (Dt)
For obtaining solutions independent of the time step in DSMC, the time step should be a fraction of the local mean collision time (Gallis et al., 2008) and the mean residence time (Espinoza et al., 2016). The hard-sphere mean collision time t c , can be calculated as (Rader et al., 2006) where k the hard-sphere molecular mean free path and v 1;prob is the most probable molecular speed. The cell residence time t res is defined as the mean time it takes a particle to traverse a cell (Espinoza et al., 2016), given by the following relation (Mata et al., 2010), where Dx is the cell size and u 1 is the average cell stream velocity of the gas. The DC micro-discharge code requires an additional restriction in the definition of the time step. The code internally calculates the electron kinetics and density based on micro-discharge theory, allowing the simulation of plasma effects on flow by the discharge phenomena in dsmcFoam+ without directly simulating the electron particles. This is made possible by accelerating the DSMC particles after the electron kinetics timing, such as the electron transport and electron-gas particle collision. Therefore, to define a Dt range for an accurate solution, values were taken from References (Mobli, 2014 andPark et al., 2019) based on the contrast between the characteristic time scale Dt char of the plasma effect with the numerical analysis of the integrated time step Dt int on Fig. 6. Therefore, to capture the discharge phenomena without the addition of electron particles in the mesh, for an accurate solution, the time step must be greater than the characteristic time Dt char of the electron kinetics, but lower than the local mean collision time and the mean residence time, i.e. 1 Â 10 À9 s 6 Dt Dt < t res Dt < t c 8 > < > :

Code validation
Cases from References (Moore et al., 2012 andFu et al., 2018b) have been utilised to verify the accuracy of the DC micro-discharge model's prediction of ion density and distribution.

Validation case 1
A 15 lm gap DC breakdown discharge based on a 1D PIC-DSMC simulation case was investigated in Reference (Moore et al., 2012). The flow of charged particles across the gap between electrode plates may be seen using a 250 by 250 cell grid with a temporal resolution of 2.66 ns. A timestep at 2 Â 10 À14 s was used, with an applied voltage of 215 V in atmospheric conditions consisting of 78.8% N2 and 21.2% O2 at 300 K. The cathode is iron, the field enhancement factor (due to surface roughness) was taken as 50, and the secondary yield was set as 0.026. Similarly, a DC micro-discharge code case with 250 cells and a cell size of 0.06 lm was used, taking the boundary conditions from the PIC-DSMC case, with the exception of the timestep being set as 2.66 ns, fulfilling the requirements for D t according to Section 3.5, an applied current of 1 lA and a work function for a 1 lm Â 1 lm iron cathode of 4.5 eV (Michaelson, 1977).
The ion and electron densities were measured in each cell across the microgap and the results are shown in Fig. 7. Dashed and solid lines represent the results obtained from the dsmcFoam+ code and Reference (Moore et al., 2012), respectively.
The code particle densities results show continuous values given by the calculation of the Maxwell-Boltzmann distribution from Eq. (19), meanwhile the 1D PIC-DSMC simulation presents the discretised solution of the densities by inserting electrons directly in the domain to interact with the gas particles. Due to the random nature of the particle behaviour (Moisan and Pelletier, 2012), breakdown (Wijsman, 1949), and field emission (Rokhlenko, 2011), it is challenging to acquire precise findings when comparing with alternative methods. The results from dsmcFoam+ show a N2 þ and O2 þ ion density peaks of 1:34 Â 10 21 m À3 and 3:60 Â 10 20 m À3 , respectively, both at 5.4 l m. The 1D PIC-DSMC solution predicts the highest density values of 2:2 Â 10 21 m À3 for N2 þ at 4.7 lm and 5:03 Â 10 20 m À3 for O2 þ at 5.5 lm. For the electron density, the PIC-DSMC results show an initial a decrease and a steeper upward tilt to the density curvature trend before reaching 5.5 lm along the gap. According to Reference (Moore et al., 2012) this occurs as a result of the Auger neutralization, where a quick stream across the gap creates ions, slowing down the acceleration during the process towards the cathode. However, the same stream neutralises some of the created ions, decreasing the number of electrons by recombination. The DC micro-discharge code calculates the total electron density from the applied current and the ionization of the gas species without any reduction in density by neutralisation interactions between charged particles and has a continuous curvature shape. After reaching a peak value at around 5.5 lm, both methods predict similar electron densities.

Validation case 2
In Reference (Fu et al., 2018b), a two-dimensional fluid model was used to analyse the Townsend breakdown DC voltages in microgaps based on Paschen's curves. A 450 lm height hemi-ellipsoidal protrusion with a 200 lm radius was set on the 1000 lm radius cathode's surface, meanwhile the anode remains as a 1000 lm radius plane. A 150 V DC voltage is applied with a secondary emission coefficient fixed at 0.1. The working gas is argon at 300 K and 500 torr.
To compare the DC micro-discharge code with the model from Reference (Fu et al., 2018b), the ion density was measured from the minimum point distance of 50 lm between the electrodes. The same boundary conditions were set as in the electrode plane plates case with additional considerations to fulfil the code's execution requirement; the electrode material was set as Au metal plates with a work function of 5.1 eV (Michaelson, 1977) and an applied current of 0.35 mA. The gas constants A and B were set as 12 and 180 respectively (Chen, 2016). Therefore, the size of each cell in the mesh was set to 1 Â 1 Â 1 lm, and the timestep was set to 1 ns, which is less than the t res value of 2.37 ns.
As shown in Fig. 8, the Ar þ ion densities from both methods have similar values. Therefore, mismatches or major changes in the ion density distribution between 0lm and 20lm along the Z-axis, as well as minor deviations between 40lm and 50lm are displayed. This is due to the fact that the values of the electric field enhancement will differ depending on the method utilised; for example, the electric field enhancement between the protrusion tip and the opposing electrode is described by an exponential distribution in Reference (Fu et al., 2018b). Meanwhile, the DC micro-discharge code takes into account the electric field with the enhancement factor as an average value along the distance between the electrodes (Z-direction for this case), where the ion allocation is a function of the Maxwell-Boltzmann distribution (Eq. (19)). Where the distribution of each method varies according to the shape of the function.

Electron temperature validation
Correction factors were required to take the theoretical electron temperature T e closer to a more reliable values T e ref satisfying the cases displayed on Table 1. Allowing for a more precise distribution curvature shape between the cell columns in the mesh for the DC micro-discharge code.
To correct the electron temperature, the most probable speed v e;i;prob from Eq. 25 is obtained. Then, the electron energy requires to be recalculated by solving the kinetic energy K E as a function of v e;i;prob and the addition of the correction factors C eÀn C p , i.e.
Thus, the recalculated electron temperature is given by, where K 0 E is the corrected kinetic energy of the electron, C eÀn is the correction factor of the transferred energy, and C p a correction factor in function of the pressure, Both factors provide a closer approach to a more realistic electron temperature by taking into account the energy exchange between electrons and neutral gas particles.
The factor C eÀn was developed based on the transferred energy factor C m during a collision between a sphere of mass m and another sphere of mass M (Franz, 2009), The adaptation of C m is based on considering each sphere to have mass in proportion to the number density of the species it represents. Therefore, one sphere represents the electrons with mass m e multiplied by the electron number density n e 0 from the applied current. Similarly, a sphere can represent the gas species with mass m gas multiplied by the number density n gas .

Ion temperature validation
By contrasting the ion energy from several sources, the DC micro-discharge code verification was possible ( Table 2). The dependence of the thermal velocity (v iþ;ave ) as a function of the ion's energy ensures a proper accelera- Fig. 8. Ar þ density distribution in a 50 lm gap Reference (Fu et al., 2018b) and the DC micro-discharge code at 500 t.orr with 150 V. Table 1 Comparison between corrected electron temperature T eC and reference values T eref with their respective boundaries conditions. Showing a close approximation between T eC and T eref . Timesteps (Dt) of 1 ns were set on cases where it was specified in the reference, denoted with (*).

Type
Gase P (torr) tion ( a ! i;j ) value on gas particles for a valid flow profile, see Eq. (37). Solutions from K E;i (Eq. (36)) without correction factors provides irrationally high energy values compared with conventional agreements (Kushner, 2004;Lin and Wang, 2015)[ (Choi et al., 2007)- (Kim et al., 2009)]. Three factors influence the proper calculation of K E;i values; (1) the stochasticity of the methods to describe and define the ion behaviour and its parameters, (2) the adaptation of discretised n e 0 , and n gas in Eq. 36, and (3) the linear increase of the electric field in the function.
1. To approximate the ion's kinetic energy values closer to conventional agreements (Kushner, 2004;Lin and Wang, 2015) [ (Choi et al., 2007;Delzanno et al., 2013;Demmons et al., 2022;Ebrahimi and Roohi, 2017;Espinoza et al., 2016;Forest et al., 2006;Franz, 2009;Fu et al., 2018a;Fu et al., 2018b;Gallis et al., 2008;Go and Venkattraman, 2014;Go and Pohlman, 2010;Godse and Bakshi, 2009;Guskov et al., 1965;Hassanvand et al., 2018;Hilal Kurt and Salamov, 2020;Jahn, 2006;Katz, 2015;Kim et al., 2006;Kim et al., 2009)], it is required to consider the stochasticity of the ion's energy by calculating the probability of its velocity by the Maxwell-Boltzmann probability distribution, where m i is the ion's mass and the temperature T i can be calculated in similar manner as Eq. (23) but focusing on the ion energy instead, where v i;prob is defined as, 2. K E;i and r are continuous equations, therefore the decrease in n e 0 and n gas by being a function of Dt for discretisation affects the consideration of the continuity of the equations. Thus, discretised n e 0 and n gas are taken as continuous values to calculate K E;i . The resultant energy is increased by low densities in a continuous equation, where particles reduce the energy loss by having fewer collisions with other particles in their surroundings. To fix this problem, a dimensionless correction factor in the form of a collision cross section was implemented, f r ð Þ ¼ 3:689 log 10 r ð Þ 2 þ 161:8log 10 r þ 1774:3 10 : 3. For the correction of the linear increase of the electric field E to a more realistic exponential increase, the electric field correction factor f E ð Þ was set as, f E ð Þ ¼ 2:11 Â 10 À9 E 1:653 1 Â 10 5 : Consequently, the corrected ion's kinetic energy is, where K 0 E;i replaces the kinetic energy from Eq. (33) with closer values to conventional agreements shown in Table 2. Table 2 shows five microgap DC discharge cases based in different works on micro-plasma effect analysis, providing the ion temperature in electronvolts (eV), making possible to contrast them with the temperature results from the DC micro-discharge code for dsmcFoam+. Corrected ion temperature T ic , from K 0 E;i found from the code, shows a closer approximation to the reference temperatures, T i ref .

Simulation cases
A 1 lm deep and 50 lm wide microchannel with two inlets and a single outlet is considered, as shown in Fig. 9. Half symmetry and the mesh was constructed using blockMesh in OpenFOAM. A 2.5 Â 2.5 Â 1 lm cell size was used, and there were a total of 4612 cells in the mesh. For all situations, the timestep Dt was adjusted to 1 ns to satisfy the restriction imposed by Eq. (45). With a cell residence time t res of 5 ns and an average collision time t c of 5.07 ns. Considering a N2 particle number density of 2:41 Â 10 20 m À3 with a k of 2:14 Â 10 À3 m and a mass of 4:64 Â 10 À26 kg. Temperature of the N2 gas was set at 300 K and an estimated speed of 422.53 m/s was computed.
Three plane-electrode configurations were set; (1) one anode/two cathodes, (2) two anodes/two cathodes, and (3) one anode (extended)/two cathodes. To simplify the naming of each configuration, they will be referred as con- Table 2 Comparison between corrected ion temperature T iC and reference values T iref with their respective boundaries conditions. Timesteps (Dt) of 1 ns were set on cases where the timestep was not specified in the reference, denoted with (*).

Type
Gas figuration 1, 2, and 3, respectively. In each configuration, two cathodes were place on the top surface of the microchannel and the location of the anodes are varied in each case. The electrode surface area was set as 35 Â 1 lm, with the exception of configuration 2, where the anode surface was defined as 105 Â 1 lm. The work function was taken as 4.5 eV as several different materials are near this value such as pristine gold (4.6), and a SEE coefficient of 0.02 to make a more general assumption on the selection of the electrode material.
The working gas was nitrogen at 300 K and three inlet pressures were tested; 10 kPa, 50 kPa, and 100 kPa. To simulate vacuum conditions an outlet pressure of 1 Pa was established, with a 200 lm wide external expansion region outside the efflux section.
Applied DC currents of 10 lA, 50 l,A and 100 lA were set for each tested voltage of 100 V, 400 V, 800 V and 1200 V for each inlet pressure case. There were therefore a total of 108 simulation cases run on a HP EliteDesk 800 G5 SFF Core i7-8700 3.2 with with 32 GB of RAM using one logical thread for each case.

Results and discussions
For the investigation of problems involving the flow of rare gases, the dsmcFoam+ employs an explicit timestepping solution with stochastic molecular collisions (White et al., 2018). Instead of trying to solve Newton's equations of motion for a very high number of individual atoms/molecules, an evolving simulation mimics the physics of a real gas in this way. Since dsmcFoam+ is a time-evolving solver, the outcomes of the DC microdischarge code likewise rely on the time dependency. To demonstrate the capabilities of the micro-discharge code, Figs. 10 shows the evolution of the N2 flow velocity profile over time from 25 ns to 325 ns, before the plasma cases become steady, after the application of high voltage at 800 V with 50lA across three electrode configurations. Separated by a 50lm gap and an inlet pressure of 50 kPa in a 1 Pa pressure environment.
The results obtained from the DC micro-discharge code in dsmcFoam+ allows the visualization of pressure, velocity, and temperature for each case, as shown in Figs. 10 and 11 for time-evolution intervals and steady state flow, respectively. Additionally, dsmcFoam+ provides the mass density, number density, rotational, translational, and vibrational temperatures for each case.
The plasma cases were performed after a continuous flow of 300 ns with a timestep of 1 ns, which is crucial to note. Then the DC micro-discharge code was activated and run until it reached a desirable interval, or up to 400ls for the analysis of the plasma steady flow cases. Three main parameters; (1) highest temperature, (2) thrust, and (3) mean velocity, were taken from the steady state simulation results to evaluate the performance and feasibility of each of the three electrode configurations along with different current-voltage (I-V) values to be employed in plasma-based micro-thruster devices. The thrust was obtained from the cross-sectional area of the outlet, meanwhile, the highest temperature is the maximum value registered at any point inside the microchannel. This is critical to determine the feasibility of the electrode configuration and the operational range for the electrical parameters (I-V) in order to avoid values near to or higher than the sublimation point of the electrodes/thruster materials or any temperature range that can involve a risk for the nearby components.
Due to the emphasis on the macroscopic description of the flow, it is crucial to highlight the several orders of magnitude discrepancies between temperatures in Fig. 12. In addition, it is important to emphasize the capability of the code to simulate theoretical cases beyond real scenarios where the temperature, pressure, and velocity rates can be greater than the material limits. Provided by the increase of I-V above the breakdown voltage, the DC micro-discharge code is not restricted by the mechanical and structural properties of the materials. Such cases can be observed at 100 kPa and 100 lA at 1200 V in Fig. 12-C, where the highest temperatures in flow reach 7:1 Â 10 5 K; 3:4 Â 10 8 K, and 1:5 Â 10 7 K for configurations 1, 2 and 3, respectively. This is a result of the voltage being 600 V above the breakdown voltage (Massarczyk et al., 2017) concentrated in a small surface area of 35 lm 2 . The temperature values are beyond the sublimation point at 1 kPa of common electrode materials, such as gold (1,663 K), silver (1,285 K), copper (1,512 K), tungsten (3,475 K) (Schmidt, 2013), where in a realistic scenario the material will be evaporated before it gets close to 1200 V under a vacuum environment.
The drastic rise in temperature at 1200 V on Fig. 12 is caused by the exponential rise of ion densities induced by the cascade effect inside microgaps where the electric field got a value of 2:4 Â 10 7 V=m, reaching non-real temperature values by increasing the voltage to double or more from the breakdown voltage point, which is around 350 V for a 50lm at 100 KPa (Fu et al., 2018a) given by Eqs. (6) and (7). Where an increased level of accelerated Fig. 11. Steady DC micro-discharge simulation results for pressure, velocity magnitude, and temperature for N2 flow at an inlet pressure of 50 kPa and an applied current of 50 lA at 800 V for the three different electrode configurations. particles with high energy (Eqs. (36)) by the electric potential raises the temperature to theoretical levels. Same case for 50 kPa at 50lm, where the breakdown is around 550 V (Hilal Kurt and Salamov, 2020), and 10 kPa at 350 V, where the voltage was set to 800 V in all cases. The highest temperature T max and thrust F values rise simultaneously as the voltage and current increase in most simulation cases, as shown in Fig. 12, as a result of the increase in potential energy by I-V. Lower temperature and thrust rates were commonly found at 100 V to 400 V, as the potential is not great enough to reach the breakdown voltage for a sustainable plasma. Presenting a steady 305 K in the three configurations, followed by 410 nN and 820 nN on thrust for 50 kPa and 100 kPa cases, respectively. These T max and F values are equal to the results from the gas flow without any electrodes, performing as similar as a cold-gas micro-thruster (Ranjan et al., 2018).
Additionally, at current values of 10 lA from 800 to 1200 V on the three configurations, the results show constant values for the range between 100 V and 400 V with a temperature value of 305 K and thrust values of 410 nN and 820 nN for 50 kPa and 100 kPa, respectively. For the 50 kPa cases at 1200 V with 10 lA, the results show a slight increase in temperature up to 309.7 K, 336.56 K, 326.8 K, and thrusts of 436 nN, 487 nN, 467 nN for configurations 1, 2, and 3, respectively. (1) one anode/two cathodes, (2) one anode (extended)/two cathodes, and (3) two anodes/two cathodes for 100 V, 400 V, 800 V, and 1200 V. Results are split by nitrogen inlet pressures of 10 kPa, 50 kPa, and 100 kPa. The applied currents of 10 lA, 50 lA, and 100 lA are distinguished by purple, turquoise, and orange lines, respectively.
Similarly, for the 100 kPa cases where the temperature reaches 306.1 K, 310.5 K, 310.1 K, and the thrust 836 nN, 852 nN, and 871 nN for configurations 1, 2, and 3, respectively. Therefore, the reason for no or a small increase of thrust at 10 lA between 800 V to 1200 V is given by the low electrical current, where the electron-ion densities(10 19 m À3 ) are not enough to sustain the breakdown even it is reached the requested voltage (600 V (Massarczyk et al., 2017)).
The electrons/ion potential energies are dissipated by a much larger neutral gas number density, an approximated 10 25 m À3 for 50 kPa to 100 kPa, and 10 24 m À3 for 10 kPa. Contrary, with current values of 50 lA and 100 lA, the electron-ion densities (10 20 m À3 ) become enough to cause an acceleration effect on the flow.
The breakdown voltage does not depend just on the pressure, gap distance and voltage as the Paschen's curve (Massarczyk et al., 2017) and Fowler-Nordheim equation (Go and Pohlman, 2010) define; it also depends on the magnitude of the applied electrical current, as shown in the simulation results. Thus, according to the simulation results, the maximum difference between the electron number density and the gas number density to provide enough charged particles with an applied potential to accelerate the flow in the microchannel becomes approximately 10 5 m À3 . Consequently, the smaller this difference becomes between the electron and neutral number densities, the greater the magnitude of T max and F.
The highest temperature values from Fig. 12 are commonly exhibited on configuration 3, followed by configuration 2, with the lowest values found for configuration 1, along with 400 V to 1200 V from 10 lA, 50lA, and 100 lA. Contrastingly, configuration 1 shows the highest thrust values, followed by configuration 2 and, with the lowest values, configuration 3. There is therefore an inverse relation between T max , and F provided by the difference in the streamlines and the vorticity inside the microchannels.
Vortex formation, induced by the displacement of the anode far of the perpendicular direction from the outlet, accelerates the flow by ion-gas particle collisions towards the cathodes. Vorticies are shown in Fig. 13-B and -C at coordinates (6:5 Â 10 À5 m; 3 Â 10 À5 m) and (17:5 Â 10 À5 m; 3 Â 10 À5 m). There are also two vortices located at the centre of the microchannel near the outlet with an approximated radius of 1:71 Â 10 À5 m for configuration 2, and 1:18 Â 10 À5 m for configuration 3. These vortices at the centre of the microchannel limit the flow towards the outlet, leaving an open area for configurations 2 and 3 of 8.6 lm and 10.8 lm, respectively, between the wall and the vortex.
The larger this open space becomes, the more volumetric flow will be transported to the outlet to be exhausted, increasing the velocity, heat dissipation, and consequently the thrust. Thus, the lack of vortex formation in configuration 1 (see Fig. 13-A) allows the flow to travel directly from the inlet to the outlet without restriction, consequently presenting the lowest temperature and the highest thrust.

Conclusion
The DC micro-discharge code implemented in dsmcFoam+ solver is an alternative method to model DC micro-plasma effects on flow with the primary focus for the research and design of plasma thrusters. The code avoids adding electrons as particles by accelerating an existing number of equivalent gas particles by calculating a numerical solution of the Fowler-Nordheim equation to emulate the ionization, reducing the computational requirements.
The applicability of the DC microdischarge code was tested by simulating three electrode configurations; (1) one anode/two cathodes, (2) two anodes/two cathodes, and (3) one anode extended/two cathodes. Inlet pressures were set as 10 kPa, 50 kPa, and 100 kPa, with an applied Fig. 13. Velocity flow profiles at steady state on 50 lm gap at 800 V with an applied current of 100 lA. The nitrogen inlet and outlet pressures were set as 100 kPa and 1 kPa; Three configurations are shown: (a) 1 anode/2 cathodes, (b) 2 anodes/2 cathodes, and (c) 1 anode extended/2 cathodes. Vortex formation is present in (b) and (c), with a distance between the wall and the vortex of 84 lm and 108 lm, respectively. voltage of 100 V, 400 V, 800 V, and 1200 V, as well as an applied current of 10 lA, 50 lA, and 100 lA in a 50 lm gap. Configuration 1 shows the most effective electrode pattern, with the lowest increase in temperature between the configurations, and also the highest velocities and thrust values.
During the demonstration of the code application, it was identified that there must be a maximum difference between the electron and neutral number densities of 10 5 m À3 between the electrodes in order to provide enough charged particles to accelerate the flow. More work needs to be done to analyse the specific amount of current density required to contribute to the change in the acceleration of the flow. In addition, the code requires a performance test, contrasting it with other PIC-DSMC methods as future work to show its performance on the consumption of computational resources. The DC micro-discharge code is a promising tool for the prototyping simulation of plasmabased micro-propulsion systems for nano-and picosatellites applications.

Declaration of Competing Interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: [Hassam Israel Guevara Jelid reports financial support was provided by National Council on Science and Technology. Hassam Israel Guevara Jelid reports financial support was provided by Fundacion ESRU.]