PIC/MCC simulation for a ns-pulsed glow discharge in nitrogen at sub-atmospheric pressure and analysis of its quasi-steady state physics

The present study investigates the physics of a nanosecond-pulsed microplasma operated at a pressure of 200 mbar with the help of a particle-in-cell simulation with Monte Carlo treatment of collision (PIC/MCC) and (semi-)analytical models. The discharge is ignited in a 1 mm gap between two parallel molybdenum electrodes by applying a voltage in the kV-range for several tens to hundreds of nanoseconds. A PIC/MCC simulation code is developed in order to describe the experiment performed under identical conditions. The simulation includes the external electrical circuit to perform ab-initio calculations of the complete experimental setup. Notable features of the PIC/MCC are (i) the adjustment of super-particle weights to reduce the computational load during drastic changes of the plasma density, which reaches values up to a few 1019m−3 , and (ii) the inclusion of dissociative recombination of the N4+ ions, which is the key loss process for the charge carriers in the plasma bulk regions. The current and voltage waveforms obtained from the simulation are compared to the experimentally measured ones and good agreement is found. After ignition, the discharge establishes a quasi-steady state exhibiting spatial features similar to a conventional DC glow discharge. Using the PIC/MCC results, reasonable approximations are identified, which allow the development of various analytical fluid models for the individual plasma regions. These models are able to reproduce the key features of the discharge in agreement with the PIC/MCC results. The simplified models for the different discharge regions can be combined to describe the global behaviour of the discharge and—in a next step—might be used to develop computationally efficient global chemistry models that account for the different power dissipation mechanisms along the discharge gap.


Introduction
During the last decades, near-atmospheric pressure plasmas gained much research interest, not least because of their high potential in numerous applications. Where lowpressure plasmas face strong limitations in the treatment of living tissues and liquids, due to their strict pressure requirements, atmospheric pressure plasmas are a natural alternative. Another important advantage of the higher pressure is the less demanding operation-no expensive vacuum equipment is necessary. On the downside, the much lower mean free path of the electrons increases gas heating. Because a low temperature is critical for many potential applications, several discharge configurations were created to ensure stable operation at low temperatures. A popular approach is the use of a dielectric in so-called dielectric barrier discharges (DBDs) to prevent conduction currents. In this way, sparks or arcs are avoided, which would lead to undesired gas heating and attrition of the electrodes. Many geometries and configurations of DBDs exist [1][2][3][4][5][6], but most share that the individual plasma pulse can only be sustained for short periods of time until the charging of the dielectric significantly reduces the gap voltage [7]. To increase the plasma-on time, while keeping the advantage of a high electron density due to fast high voltage pulses, ns-pulsed DBDs can be combined with a continuous RF-driving voltage [6,8].
An alternative to DBDs are discharges with bare metallic electrodes, which allow DC conduction currents. Examples are pin-to-pin or pin-to-plane geometries, operated with nanosecond or microsecond high-voltage pulses [3,4,9]. While bare metallic electrodes in plane parallel geometry are often associated with radio frequency plasmas [3,4,10,11], it was demonstrated that they can also be operated with nanosecond high-voltage pulses-producing a homogeneous glow discharge [12][13][14][15][16][17]. Similar to DBDs, here approaches are proposed to prolong the effective life time of the discharge by applying a constant lower voltage bias in addition to the transient high voltage pulses [18].
These homogeneous glow (non-DBD) discharges are the focus of the current work. For the investigation, a simple discharge geometry is chosen to allow one-dimensional modelling. With this geometry and the electrical circuit described in section 2, discharge pulses with a repetition rate of several kHz and a duration of several 100 ns can be achieved [16]. Furthermore, during most of the pulse duration the plasma bulk parameters-plasma density, electric field-are constant in time [16,17] and uniform in space over a considerable domain of the system [15]. These constant conditions make the plasma source presented here a valuable tool for the research of fundamental plasma processes, but also for the development and benchmarking of new diagnostics. This work investigates the processes leading to the constant plasma conditions in the so called quasi-DC phase.
The discharge in the quasi-DC phase shows many similarities to the usual DC discharges in the anomalous glow regime, i.e. there are spatial structures analogous to the cathode sheath, the negative glow, Faraday dark space, positive column and, finally, the anode sheath. As it is known, the non-local electron kinetics in the cathode region determine the discharge behaviour [19][20][21]. For this reason, we apply a particle based method (particle-in-cell simulation with Monte Carlo collisions, PIC/MCC), to correctly capture the discharge dynamics. This technique is described in section 3. In section 4 results of the simulation are presented and discussed. Analytical and numerical models are developed for the specific discharge regions, which enable an insight into the basic physics governing the observed phenomena in the various spatial regions and temporal phases of the discharge. These models are based on both, the fluid dynamic picture and the kinetic picture. Separate models are developed for the sheath, the negative glow/Faraday dark space, and the bulk region.
Finally, we would like to note that further validation is performed by extensive optical emission measurements, which are compared with the PIC/MCC results [22]. Excellent agreement is found throughout.

Experimental setup
The discharge is ignited between two plane parallel molybdenum electrodes, which are arranged such that they create a discharge channel of 20 mm length with a 1 mm × 1 mm cross-section. Glass plates ensure a gas flow along the extended axis-i.e. the y-axis-of the discharge channel, as it is shown in figure 1. The gas inlet is positioned in the centre of the discharge, supplying 20 sccm of nitrogen at a pressure of 200 mbar. Under the operational conditions used in this work, the gas temperature was T g ≈ 330 K [16].
An electric circuit (see figure 1) is employed to apply short high-voltage pulses to the cathode with negative polarity, a duration in the order of hundreds of nanoseconds, and a repetition rate of 1 kHz. The repetitive pulsing ensures that remaining charge carriers from the previous discharge pulse ease the ignition of the next. A capacitor (C 0 = 4.7 nF) at the input of a Behlke HTS 81 high voltage switch is charged to 3 kV by a DC power supply (Heinzinger LNC 6000-10neg). When the switch closes (its internal resistance is estimated to be about 12 Ω), the voltage at the electrode increases until the critical value is reached where the gas breakdown occurs. Once the discharge is established, the stray capacitance-estimated to be about C p = 40 pF-is discharged through the plasma, causing an initial peak in the current. After this initial peak, the current through the plasma is maintained by discharging the capacitor C 0 via the series resistor, R s = 255 Ω. An additional parallel resistor, R p = 2 kΩ, is used to ensure that no charges remain in the circuit when the discharge is extinguished. The voltage over the discharge gap is measured by a high voltage probe (Lecroy PPE6kV) and the current through the discharge by a current probe (Tektronix CT-2).

Discharge model and kinetic simulation approach
The plasma formation in the microdischarge is described numerically by an electrostatic 1d3v (one dimension in space, three dimensions in velocity space) particle-in-cell simulation complemented with the Monte Carlo treatment of collision processes, i.e. by the 'PIC/MCC' approach, see e.g. [23][24][25][26][27][28][29]. Being a particle based approach, PIC/MCC is fully capable of capturing kinetic effects and allows the computation of the distribution functions of the charged species.
The PIC/MCC scheme uses a computational grid, which in the simplest case is equidistant with a division ∆x. The potential, the electric field, and the charged particle densities are defined at the grid points. Plasma characteristics, like mean particle velocities, rates of elementary processes, etc, are also computed at these grid points. The simulation is selfconsistent, i.e. the charged particles move under the effect of an electric field that is defined by the potential at the electrodes and the presence of space charge. The approach uses the concept of 'superparticles' that represent a high number of real (charged) particles in the system. Besides space, time is also discretised in the simulation, the particle trajectories are obtained by the numerical integration of their equations of motion. As the method used here is standard in most respects, no detailed description of the simulation is provided here, details can be found in the papers cited above. Further, the core of the code used in this work was benchmarked in [30]. Here, we discuss only the specific features and parameters that are relevant for this work.
The model of the discharge is based on a simplified set of elementary processes in the nitrogen plasma. These processes are, however, expected to reproduce correctly the main physical effects taking place in this system. For the electron-N 2 collisions, we use the same cross-section set as presented in table 1 of [31]. The processes considered include elastic scattering, rotational, vibrational, and electronic excitation, as well as ionisation. To compute the cross-sections at any given energy, we use linear interpolation between the values given in the input data file.
The N + 2 ions created in ionising collisions rapidly convert to N + 4 ions via the three-body reaction N + 2 + N 2 + N 2 → N + 4 + N 2 , which has a rate coefficient of k 2→4 = 5 × 10 −29 cm 6 s −1 [32]. For both ion species we consider only elastic collisions with the background N 2 gas. For N + 2 -N 2 collisions we adopt a differential cross-section that is composed of an isotropic part and a backward scattering part (the scattering angles are to be understood in the centre-of-mass frame), in accordance with [33]. At low ion energy the cross-section converges to the Langevin limit, while at high energies the backward scattering (symmetric charge exchange) becomes dominant. For the N + 4 -N 2 interaction we consider the Langevin cross-section.
The probability of any collision process is given as where n 0 is the gas density, σ t is the total cross-section at the centre-of-mass energy ε, g is the relative velocity and ∆t is the time step for the given species. One of the constraints for the stability and accuracy of the PIC/MCC scheme (see e.g. [34]) requires that P coll is kept at a low value to prevent the occurrence of multiple collisions within a ∆t time step. Due to the relatively high pressure of 200 mbar the plasma is highly collisional and this condition dictates an electron time step of ∆t e = 2 × 10 −14 s. Significantly longer time steps (∆t i = 20∆t e ) are used for the ions ('subcycling'). Other stability and accuracy conditions require that (i) the spatial grid resolves the Debye length, (ii) the time step resolves the plasma frequency, and (iii) most of the particles traverse less than ∆x during ∆t (known as the Courant-Friedrichs-Lewy condition) [35]. The collision probability (1) is calculated in each time step for all the particles (i.e. in each ∆t e for the electrons and in each ∆t i for the ions), based on the total cross-section of collision at the actual energy, ε. For the electrons, the cold gas approximation is adopted, therefore ε equals the electron energy and the relative velocity g equals the electron velocity. For the ions, random collision partners are sampled from the background gas with Maxwell-Boltzmann velocity distribution for the calculation of the actual values of ε and g. Comparison of P coll with a random number, R 01 , having a uniform distribution over the [0, 1) interval decides about the occurrence of the collision. If R 01 < P coll , a collision of the given particle takes place, the type of which is selected based on the values of the cross-sections of the possible reactions with the help of another random number. The collision is executed in the centre of mass frame (see e.g. [34]). All collision processes are assumed to result in isotropic scattering, except for the electron-impact ionisation process, for which the collision algorithm proposed in [36] is utilised.
The discharge cell, as explained in the previous section, is a part of an electrical circuit that contains a number of passive elements and a switch. This circuit has to be considered in the simulations, too. At the start of the simulation, the capacitor C 0 (4.7 nF) is assumed to be fully charged to the source voltage U s = U s,0 and the capacitor C p (for which we take 40 pF in the simulations) is assumed to be completely discharged. When the switch connects C 0 to the series resistor R s , current (I Rp ) flow starts through R p and through the parallel capacitor C p (I Cp ). The voltage on the discharge cell, U d , equals the voltage of the capacitor C p , and I Rp = U d /R p as well as C p dU d dt = I Cp . When the voltage of the latter reaches the breakdown threshold, a current I d starts to flow through the discharge. These currents cause a change of U s according to C 0 dUs dt = I Rp + I Cp + I d . The discharge current I d is composed of the currents carried by the various particles and the displacement current; these quantities are evaluated at the powered electrode. The above equations for the external electrical circuit are solved in conjunction with the Poisson equation.
As already mentioned above, the PIC/MCC scheme uses superparticles that represent several real charged particles expressed by the superparticle weight, W. While for many types of systems W can be specified as a constant, in settings where the density of the charged particles, and consequently the number of superparticles, grows exceedingly during the time scale of interest, some kind of adaptive particle weighting [37][38][39][40][41] has to be applied to improve the local/global accuracy of such simulations. Without such an approach, i.e. using a constant superparticle weight, we would face the problem of either (i) having very bad statistics in the early phase of the pulse, or (ii) reaching a huge number of superparticles in the quasi-stationary phase of the pulse. Therefore, in our code we implement the following approach: • Initialisation: the simulations start by seeding a given number of superparticles (both electrons and N + 4 ions) with a given weight W 0 (here: W 0 = 4167), with uniform spatial distribution within the electrode gap. These quantities define the initial density of the electrons and the ions. • Weight adjustment: whenever the number of electron superparticles grows by a factor of 1 + ν from the moment of the last weight adjustment, or from the beginning of the simulation, we increase the weight W for new particles by a factor of 1 + ω. We emphasise that this weight adjustment affects only the 'new' electrons and ions born in ionisation events in the gas phase and the 'new' electrons created in secondary emission events at the electrode surfaces (see below). The weights of 'old' particles remain unchanged, therefore superparticles with various weights coexist in the system. As an advantage, no particle merging has to be applied. • Ionisation: whenever an electron with a weight W 1 participates in ionisation, and the actual weight for new particles is W, a new electron and a positive ion (N + 2 ) with a weight W will be created only with a probability of P i = W 1 /W. • Secondary electron emission due to ion impact: when an ion with a weight W 1 reaches an electrode, a new secondary electron with weight W is emitted with a probability P s = γW 1 /W, where γ is the secondary electron yield due to ion impact (taken to be 0.04 for both N + 2 and N + 4 ).
The simulations, of which the results will be presented here, are initialised with 600 000 superparticles, which represents a real initial density of n −,0 = n +,0 = 1.25 × 10 11 cm −3 . This value for the initial density is based on order of magnitude estimations and is finally chosen such that the current waveform in the simulation matches with the one in the experiment. In any case, it was found that only the ignition phase is affected by the value of the initial density, but not the steady state. These particles mimic the residual charges that are present in the system after the previous pulse in the experiment. (Repetitive pulses are not described in the simulation due to the very long time between the pulses in the experiment.) Still, the initial density used here should be understood as a rough estimate. The details of the ignition depend on multiple other factors, like the exact voltage rise-time, the spatial distribution of the initial charges and the presence of excited gas molecules. Because the focus of this work is on the quasi-DC phase and the general behaviour during the ignition is well reproduced, these aspects are not investigated here. The parameters for the weight adjustment procedure (explained above) are ν = 0.025 and ω = 0.15. With these settings the number of electron and ion superparticles grows, respectively, by a factor of ∼1.4 and ∼10 during the whole time span of the simulation, i.e. we have about 0.8 million electron and 6.5 million N + 4 ion superparticles at the end of the runs. Due to the efficient ion conversion process, the number of N + 2 ion superparticles is in the order of 10 5 only. The simulations also include the recombination process e + N + 4 → N 2 + N 2 . It is treated in the electron transport part of the simulation, and is added as an electron-impact process with a cross-section σ r [m 2 ] = 7.7 × 10 −38 /ε [J], derived on the basis of the rate constant given in [32]. The probability that an electron participates in this process is computed based on equation (1), with substitutions n 0 → n 4 , σ t → σ r , and g → v e , where n 4 is the local density of N + 4 ions and v e is the velocity of the electron. Whenever a collision occurs in a given spatial (grid) cell, the weights of all electrons and all N + 4 ions in that cell are adjusted in a way that their sum (separately for the electrons and the ions) decreases by W e , the weight of the colliding electron. This approach, while mimicking that an electron and an ion superparticle are lost, actually avoids deleting these particles from the simulation and ensures that the recombination process decreases equally the local densities of the ions and the electrons.
The number of grid points in the simulation is N g = 4000, and a run comprises N t = 2.5 × 10 6 time steps for a time domain of T = 50 ns. A run lasts about two weeks on a Xeon Gold 6132 CPU. For final results the data obtained from 15 independent simulation runs-with different random placements of initial superparticles and different initial random number generator seeds-are averaged. Upon the presentation of the results these averaged data will be used.
The quantities 'measured' in the simulations include the densities, the mean velocities and the current densities of the electrons and ions, the electric field and the potential, as well as the rates of individual elementary processes, space and time resolved. Upon collecting information about these quantities, data are binned for 4000 electron time steps, while the full resolution of the spatial grid (with 4000 points) is retained. As the PIC/MCC scheme uses the leapfrog integration technique, particle positions are known at integer time steps (i.e. at integer multiples of ∆t), while velocities are defined at times shifted by ∆t/2, see, e.g. [34]. Therefore, for a proper calculation of the current densities average values of the velocities of the respective particles at integer multiples of ∆t has to be computed [42]. Besides the quantities mentioned above, the electron impact excitation rate to the N + 2 (B 2 Σ + u ) state (to allow comparison with the measured optical emission from this state) is also computed, and the electron energy probability function is saved with space and time resolution as well. This set of results allows a thorough comparison with the experimental results [22].

Current and voltage: PIC/MCC simulation and measurements
As initial validation of the discharge model, the resulting voltage and current waveforms are compared to measured waveforms from the experiment in figure 2 for U s,0 = 3000 V. While there are some differences at early times, e.g. slightly different voltage rise times-probably caused by non-ideal behaviour of the switch-after the time of the breakdown the agreement for both, current and voltage, is good. This is a first indication that the numerical model can reasonably reproduce the experiment. Further validation is performed in [22] by measurements of the optical emission.
In figure 2, two temporal regimes are notable. First, the very short ignition phase with relevant time scales in the order of a nanosecond. Second, the longer quasi-DC phase with time scales of some 10 ns to several 100 ns. During this phase, the spatial structures of the discharge developed during the ignition stay mostly unchanged.

Ignition phase: PIC/MCC simulation
While the main focus of this work is on the quasi-DC phase, a qualitative discussion of the ignition dynamics still provides interesting insights. Regarding the temporal development of the current, the ignition is divided into three phases (see figure 3(a)).
In the beginning of the first phase (t < 10 ns, figure 3(b)), the initial electrons (which remain from the previous pulse) are pushed to the anode by the applied electric field, resulting in a positive space charge on the cathode side. This process does not lead to any significant current, but the positive space charge slightly increases the amplitude of the electric field closer to the cathode. The secondary electrons emitted due to the initial ions hitting the cathode multiply in this electric field. In contrast to the electrons, the ions produced in these ionisation events stay close to the position, where they were created over the time scale of the pulse. Such, the positive space charge close to the cathode is increased further. Therefore, the electric field is shielded over a shorter distance, leading to a higher field and stronger ionisation. Because of this positive feedback the positive space charge increases more rapidly. This loop results in a cathode directed ionisation wave, which produces an increasing current that-at the cathode-is mostly given by the displacement current due to the increasing field at the cathode.
In the second phase ( figure 3(c)), the ionisation wave comes to within 100 µm from the cathode. The field has reached values, where the ionisation does not increase significantly anymore. Additionally, the spatial extension of the high field region shrinks, decreasing the secondary electron multiplication. These two factors slow down the ionisation wave, so that the field at the cathode increases only linearly with time. The current, which at the cathode is still given by the displacement current up to this point, experiences a plateau.
During the third phase ( figure 3(d)), the ion density at the cathode has increased enough, that the conduction current becomes significant. At the time, when the ion conduction current starts to increase (t ≈ 13 ns), the positive space charge region at the cathode collapses rapidly (compare the field profiles at t = 13.0 ns and 14.5 ns). This collapse is also clearly visible in the emission measurements in [22]. At about t = 17 ns the current reaches its maximum, and consists essentially only of the conduction current of the ions. The displacement current becomes insignificant as a steady cathode sheath develops.
The transition in the current composition at the cathode surface from displacement current (phase (a) and (b)) to ion conduction current also occurs in streamer discharges as they are observed in various dielectric barrier configurations (see for example [43][44][45]).
Finally, after a small overshoot of the current due to the stray capacitances in the setup, the quasi-DC state is reached.

Quasi-DC phase: PIC/MCC simulation and (semi-)analytical models
During the transition from the ignition to the quasi-DC phase distinct spatial regions are formed, which are visible in both, the density profiles of the charged species and the electric field (see figure 4). Please also refer to the structure of the light emission by the discharge presented in the related paper [22]. These spatial regions are very similar to the ones observed in a traditional DC glow discharge, and we will use an analogous nomenclature. In this work the regions are called-from cathode to anode: the cathode sheath (CS) with a very high electric field, vanishing electron and moderate ion densities; the negative glow (NG) with a peak in the plasma density and a vanishing electric field; the Faraday dark space (FDS) where the plasma density transitions with an exponential profile to the value in the plasma bulk (positive column) and the field is  partially reversed; the plasma bulk with uniform densities and field; and finally the anode sheath (AS).

Cathode sheath.
A zoomed view into the cathode sheath region is shown in figure 5 for the charge carrier densities and the electric field. The electric field strength increases from the sheath edge to the electrode in an approximately square root shape and reaches values of up to 5.7 × 10 7 V m −1 (which corresponds to approximately 13 kTd).
From the density profiles of the charge carriers it is clear that the electrons do not contribute to the charge density in the sheath. In contrast to the other parts of the discharge, where the positive charge is mainly given by N + 4 ions, both ion species (N + 2 and N + 4 ) are prevalent in the cathode sheath layer.
The electron current density is very well reproduced by a spatial exponential growth with a constant (average) Townsend coefficient α: where j e,c is the electron current density value at the cathode surface. The space dependent Townsend α coefficient is not an output of the simulation, but can be calculated from the time averaged ionisation rate, R iz via In the bottom of figure 6, it can be seen, that the assumption of a constant α-here, the average value of α(x) over the sheath-is quite reasonable. Also, the numerical value of α = 1.6 × 10 5 m −1 agrees well with the high field limit of the Townsend coefficient available in the literature [46]. The PIC/MCC model described in section 3 is able to provide a very detailed description of the discharge dynamics and-as is seen in the course of this work and in [22]agrees very well with the experiment. But on itself, the simulation returns only numerical values. To obtain results for different discharge parameters or a slightly modified setup, it is therefore necessary to either rerun the simulation for the new parameter set or to extrapolate the results already obtained for other conditions. Because the time needed for a single simulation run is in the order of weeks, an extrapolation with analytical models can save a large amount of computation time. Furthermore, from analytical models physical insights can be extracted that are not accessible through numerical simulations, e.g. the scaling of various quantities with the input parameters. Therefore, in this work analytical descriptions are developed wherever possible and compared to the PIC/MCC results. While the electron kinetics in the cathode sheath are generally non-local due to the high field gradient, the approximation of a constant Townsend ionisation coefficient can be used to derive a fluid model for the ions. With the previously determined constant α (see figure 6) the amplification of the electron flux density in the sheath is given by For a given secondary electron coefficient γ, the distance from the cathode at which the electron current equals the total discharge current is s = ln ((1 + γ)/γ) /α. In the frame of this model s is assumed to be the sheath thickness. Furthermore, the coordinate system is chosen such that ξ = 0 is the position of the sheath-plasma boundary and ξ = s is the position of the cathode. Note that the coordinate system therefore differs from the coordinate system of the PIC/MCC, where x is the distance from the cathode. Further, in the following discussion the electron flux density is defined as positive, as it enters the relevant equations only in the ionisation source terms. There only the amplitude of the electron flux density is essential, not the direction. The electron flux density is then with the values at the boundaries Γ e (0) = I eA and Γ e (s) = I eA where A is the electrode surface. The fluid model for the ions includes-in addition to the production of N + 2 ions by ionisation-the conversion of N + 2 to N + 4 with a frequency of ν 2→4 = k 2→4 n 2 0 calculated from the same rate coefficient k 2→4 = 5 × 10 −29 cm 3 s −1 and gas density n 0 as is used in the PIC/MCC simulation. The corresponding fluid equations are then: with the electric field determined by Poisson's equation Here, m iσ , Γ iσ , u iσ with σ = 2 (for N + 2 ) or σ = 4 (for N + 4 ) are the mass, the flux density and the average velocity of the ions, respectively. E is the electric field, e the elementary charge and ε 0 is the vacuum permittivity. In (8) the contribution of the electrons to the electric field in the cathode sheath is neglected due to their low density.
The momentum transfer collision frequency for N + 4 is assumed to be constant with ν m4 = 1 ns −1 (at the given pressure of 200 mbar). For the N + 2 ions the polarisation collisions have a higher frequency ν m2 = 3.5 ns −1 and ion charge exchange collisions are considered to be important with a mean free path of λ cx = 0.76 µm. These parameters are chosen to be consistent with the cross-section input in the PIC/MCC simulation. Equations (7) and (8) can be normalised to the natural units of the problem: with the normalisations: are compared to the PIC/MCC results in figure 7, for U s,0 = 3000 V. To avoid singularities, the initial conditions for the fluxes and velocities are not chosen to be exactly zero, but 'sufficiently' small. Empirically, it was seen that the choice of initial conditions is not critical for the overall profile farther away from the sheath edge. Generally, a good agreement can be seen for all modelled quantities. The largest deviations are found close to the sheath edge and can certainly be attributed in part to the simplified boundary conditions. Another reason for deviations of the fluid model close to the plasma-sheath boundary is the neglect of the bulk electrons. Those can penetrate into the sheath for a short distance and shield the positive space charge at the boundary. In principle the agreement can be improved by choosing more precise boundary conditions and including the bulk electrons into Poisson's equation with the corresponding Boltzmann factor. But to follow this approach detailed knowledge is necessary, which for now is only available from the PIC/MCC results (e.g. the electron density and temperature of the bulk electrons). Because the fluid model is intended to rely as little on the particle simulation as possible, it is not adjusted further and the only input from the PIC/MCC is the (approximated) Townsend α coefficient.
As can be seen in figure 7, with the model developed here the ion behaviour and the electric field profile in the cathode sheath can be estimated well without the need to perform a full PIC/MCC simulation (the numerical solution of the fluid model takes less than a second on a standard laptop). Furthermore, the model can be used to gain additional insights. It is clear that the ion composition in the sheath is important. Here, both ion species have comparable densities, but their dynamics are completely different. It can be seen already by inspecting the velocities shown in figure 7 that N + 2 and N + 4 ions have different collisionality. This is emphasised in figure 8, where the various terms in the momentum equations for the individual species are compared for the case investigated in this work. While for the N + 2 ions the collision terms-in particular the charge exchange collisionsdominate, for the N + 4 the largest contribution is by the inertia term. The significance of the inertia term makes the development of analytical models very difficult, such that we were not able to find an analytical solution for the case investigated here. An approximate analytical solution was found for the special case that both ions are highly collisional and that the conversion happens with a constant mean free path. This solution is not applicable here, and therefore only presented in the appendix for the interested reader. Furthermore, it is worth to note that the drift diffusion approximation often used in fluid simulations is in general not valid in the cathode sheath at the pressure investigated in this work.
The assumption of a constant ionisation coefficient, is reasonable for the given conditions. Also, the value of the average α should be mostly set by the neutral gas density, but it might be desirable to check these assumptions. For this purpose, a simple single particle Monte Carlo (MC) simulation is sufficient. A possible approach could be to first assume a constant α = 1.6 × 10 5 m −1 as was observed under the conditions here. Then the MC simulation can be run for the resulting field profiles to check the assumed ionisation coefficient and correct it if necessary for an additional run of the fluid model. This approach is similar as the one used previously in DC discharges [20]. For the frame of this work, a very simple MC simulation is used, which tracks single electrons emitted at the cathode and those produced in the sheath for a given space-dependent electric field. The collision processes are the same as in the PIC/MCC simulation. By counting the ionisation processes in given grid intervals, the Townsend α coefficient can be obtained. Finally, it is observed that the α coefficient does not vary strongly with the applied voltage. This is demonstrated in figure 9 where the sheath voltage in dependence of the discharge current is plotted for α = 1.6 × 10 5 m −1 and for the iterative procedure explained above. It can be seen that the correction due to slightly different Townsend α coefficient is low.

Negative glow and
Faraday dark space. The regions following after the sheath edge are called negative glow and Faraday dark space-analogous to the DC glow discharge. In these regions the electric field drops close to zero, but the secondary electrons still have enough energy for ionisation. Because transport of the ions can be neglected in the Faraday dark space, in the steady state the densities of the charged particles are determined by the balance of production and loss processes for e − , N + 2 and N + 4 : Figure 9. Sheath voltage in dependence of the discharge current. The line represents the (fluid) model calculations for a fixed Townsend α coefficient α = 1.6 × 10 5 m −1 . The symbols denote the sheath voltage that is obtained for an iterative correction. First the fluid model is solved for the fixed ionisation coefficient above. For the resulting field profile the new coefficient is determined with a MC particle simulation. Finally, the fluid model is solved a second time for the corrected α.
This set of equations is easily solved: Here, it should be noted that the ionisation rate cannot simply be written as R iz = ν iz n e , because, here, n e denotes the bulk electron density, while the ionisation happens through the secondary electrons from the cathode. The densities obtained by balancing the production and loss processes are shown in figure 10 for the region between the sheath edge (s ≈ 20 µm) and the bulk (x ≈ 200 µm). As can be seen the agreement is very good in most of the Faraday dark space. Deviations close to the sheath edge can be attributed to the neglect of the flux derivative. The smaller deviations at the transition to the bulk are most likely due to the slight change of the densities in time, which can be seen in figures 11 and 12 in the next section. Finally, the importance of recombination processes should be emphasised. Because the plasma in the negative glow is essentially externally sustained through the secondary electrons produced in the cathode sheath [19], the ionisation rate is mostly uncoupled from the conditions in the negative glow. Because losses through diffusion are negligible on the time scales in this work, other processes are necessary to balance the production through ionisation. Under the conditions of this work, the dominant process is dissociative recombination with the N + 4 ions (the predominant ion species). PIC/MCC simulations without including recombination resulted in a plasma density in the negative glow more than an order of magnitude higher than in figure 10 and still increasing roughly linearly with time at t = 50 ns.

Bulk.
At a pressure of 200 mbar, the bulk covers about 80% of the discharge volume and has homogeneous electron density and electric field profiles (see [15,17,22]), which do not vary drastically after the initial phase of the discharge. These properties make the bulk an interesting target for zero-dimensional plasma models. In figure 13  the electron energy distribution function (EEDF) at the discharge centre during two different times-t = 30 ns and t = 50 ns-is compared to results obtained by the Boltzmann solvers, Bolsig+ [47] and MultiBolt [48] for E/N = 94.5 Td and E/N = 98.4 Td, respectively, which are obtained in the PIC/MCC simulations at these conditions. Generally, the agreement is very good. This suggests that non-local effects visible in the cathode sheath and the NG do not influence the EEDF in the bulk and the bulk electrons are in local equilibrium with the electric field. Therefore, if the bulk field is known, a large portion of the discharge can be simulated efficiently by the use of Boltzmann solvers. But note that results of Bolsig+ [47] (two-term approximation) deviate from the PIC/MCC and MultiBolt [48] (six-term approximation) EEDFs. Pitchford et al [49] discussed this precision loss of the two term approximation for elevated electric fields under similar conditions as are given here. While seemingly small, the difference can have a significant impact on all rates with high energy threshold, such as electronic excitation and ionisation.
The local electron kinetics in the bulk, can now be used to estimate the bulk electric field and electron density. Assuming that the conversion, N + 2 + N 2 + N 2 −→ N + 4 + N 2 , happens fast enough, the densities of the electrons and the cluster ions should be roughly equal, n e ≈ n i4 , as can also be seen in figure 4. Using the spatial homogeneity in the bulk, for a steady state ν iz (E) = n e k rec (E).
Additionally, the current in the bulk is mainly given by the electron conduction current. Therefore, the electron density in equation (14) can be eliminated and one obtains a (non-linear) equation for the electric field: Once the electric field is known, the electron density can be calculated with (14). Both quantities, determined for the steady state, are compared to the values in the PIC/MCC (averaged over the region x = 0.5-0.8 mm) in figures 11 and 12. The overshoot of the electric field and the density directly after the ignition is due to the high current, which is caused by the discharging of the stray capacitance. Clearly, the steady state is not entirely reached in the 50 ns covered in the simulation, but the variation is relatively small. Note that because the field is constant along the bulk, for a known bulk length L b the voltage drop U b is given by U b = EL b . This yields the current-voltage characteristic I(U b ) of the bulk or-by solving the equation (15)-the bulk voltage dependence on the discharge current.

Current-voltage characteristic of the discharge setup.
Because the voltage drop over the Faraday dark space and the anode sheath is small, the total discharge voltage is essentially given by the voltage drops over the cathode sheath and the positive column: where U s is the cathode sheath voltage and U b the bulk voltage. The dependence of U s and U b on the discharge current, I d , is calculated with the models in the previous sections. The applied voltage, U 0 , drops over the series resistor with resistance R s and the discharge Additionally, the current through the series resistor consists of the current through the discharge plus the current through the parallel resistor, I Rp = U d /R p , such that we obtain This finally yields the current-voltage characteristic which can be solved (numerically) for the discharge current. This current voltage characteristic is able to reproduce the DC current (7 A) for the 3 kV applied voltage in the PIC/MCC. Additionally, it is compared to measured current voltage characteristics in the experimental setup for different applied DC voltages in figure 14. The agreement between the model and the measurements is reasonably good, allowing to extrapolate the theoretical results obtained in this work for a single condition to different applied voltages. Furthermore, it can be seen that including the spatial profile of the ionisation coefficientgreen dots in figure 14, see also section 4.3.1-leads to just a slightly better agreement with the experiment.

Conclusions
In this work we presented particle simulation results for a nspulsed discharge operated in nitrogen at 200 mbar. Because of the pulsed nature of the discharge, the plasma density changes drastically in the time frame of the simulation. To avoid a severely increasing number of superparticles and the resulting computational load, a weight adjustment scheme was included. The experimental current and voltage waveforms were well reproduced by the PIC/MCC simulation, which is an indication for the correctness of the discharge model and its implementation in the simulation code. During the ignition, three phases were observed. In the first two, when the ionisation wave travels from the anode to the cathode, the discharge current is given by the displacement current. As the ionisation wave approaches the cathode, it stops. If the ion density at the cathode is not high enough to supply a significant conduction current during this time, the discharge current reaches a plateau. Finally, when the ion density is high enough, the conduction current takes over and the ignition can proceed.
After the ignition, which only takes a few nanoseconds, a steady quasi-DC phase was found to be reached. The spatial profiles of different quantities, like particle densities and field, are essentially constant during this phase and show strong similarities to the traditional DC glow discharge. But in contrast to low pressure DC discharges, here the production of charged particles is not balanced by diffusion losses. Instead, dissociative recombination is the dominant loss process.
In the sheath the Townsend ionisation coefficient, α, is not in local equilibrium with the field due to the high electric field gradients. Nonetheless, for the high fields in this work, an average α was found to reproduce well the growth of the electron current within this region.
The secondary electrons produced in the cathode sheath act as a steady ionisation source in the negative glow, even though the field there is very low. The low field and theconsequently-cold electrons allow for efficient recombination, which finally balances the ionisation via the secondary electrons.
The non-local behaviour in the sheath and the negative glow was found to diminish within the Faraday dark space. Therefore, the bulk, which makes up about 80% of the discharge volume under the given conditions, can be modelled by local Boltzmann solvers. In agreement with [49], it was seen that a solver which uses higher order terms, yields a better agreement with the EEDFs obtained from the particle simulation. The difference was especially pronounced in the high energy part of the distribution, so that the impact is largest for excitation and ionisation processes.
Combined, the simple models for the different discharge regions can be used to calculate the current-voltage characteristic of the setup. The resulting characteristic agrees reasonably well with experimental measurements, so that in this way the results of the simulation can be extrapolated to different applied voltages.
Further, the steady-state fluid models can be used in global chemistry models for ns-pulsed glow discharges, as they provide valuable information on the spatial power dissipation profile along the discharge gap. Incorporating the additional information in chemistry models will produce more reliable results-without increasing the computational load as for example coupling to a PIC/MCC or a fully time resolved fluid model would do.
Finally, we note that the models and procedures presented in this work can readily be applied to different gas mixtures.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors. within the framework of the CRC 1316 'Transient atmospheric pressure plasmas-from plasmas to liquids to solids'. Z D is grateful for being supported in part by the Mercator Fellowship program within the CRC1316.