Numerical simulation of fast atmospheric pressure discharge in gas diode with plane-grid cathode system

In this paper we perform an accurate time-dependent finite-element numerical gas discharge simulation in two-dimensional semi-periodic computational geometry of the gas diode with a plane-grid cathode system. The diode configuration we investigate is similar to previously studied experimentally. Discharge simulation is performed in the framework of two-moment macroscopic (hydrodynamic) discharge plasma model accounting photoionization and autoelectronic emission from nonuniform electrodes surfaces. The results of numerical calculations can be used for further estimations of a runaway electron flows characteristics.


Introduction
The intensive development of modern experimental techniques for the investigating of fast processes in various gas discharges makes it possible to discover a large number of practical applications for a wide class of fundamental physical phenomena. One of the most significant fundamental phenomenon in gas discharge physics is the electron runaway in high-pressure gases [1]. It is represented by the generation of high-energy electron flows that proceed into the regime of continuous acceleration by overcoming the decelerating force caused by various collisions with the molecules or atoms of gaseous medium. The runaway electrons appearance in the Earth atmosphere was predicted first in 1925 [2], and much later, the theoretical prediction was successfully proven in laboratory studies [3]. In nanosecond gas discharges, the runaway electron beams are high-current electron flows having a subnanosecond duration with wide energy spectra [4]. This allows using fast gas discharges as a prospective compact sources of ultra-short powerful X-ray pulses. According to this a high-voltage nanosecond discharge must occur in a specific configuration of gas diode providing high output characteristics of the output runaway electron beam.
In numerous experiments different designs of gas diodes were studied [1]. The three-electrode diode consisting plane cathode, grid cathode and plane anode was found to be the most promissory w.r.t. feasible applications. In the experimental paper [6], it was shown that the runaway electron flow through the anode foil in this configuration has record current values. Moreover, it was shown in [6] that the runaway electron beams were also registered in the backward direction, i.e. beyond the cathode. This in itself is of fundamental interest, since the backward motion of plasma electrons to the cathode, and furthermore into the quasi-equipotential region between plane cathode and cathode grid, is unlikely.
The main objective of this theoretical study is to investigate the instant characteristics of the discharge plasma in the diode with plane-grid cathode system in order to provide input parameters for the hybrid kinetic model interface for the runaway electron beam parameters estimations. In this work, the two dimensional semi-periodic macroscopic fluid model was developed investigating the discharge dynamics in the diode with plane-grid cathode system filled with atmospheric pressure nitrogen with small (less than ~1%) oxygen admixture. Further evaluation of runaway electron beam parameters can be provided by using the hybrid kinetic approach described in [7].

Model development
The schematic of the gas diode geometry with a planegrid cathode system being investigated is depicted in Figure 1 [6]. It consists of two plane electrodescathode and anode, and cathode (grounded) grid made of parallel thin wires of the radius r arranged in parallel to plane electrodes. The anode is connected to the highvoltage nanosecond power source through the 100  ballast load, so the maximum voltage applied to the diode is 10-15 times greater than the static breakdown value. For example, in experimental paper [6] typical diode parameters were l = 3-12 mm, d = 4 mm, r = 0.1 mm, L = 4-8 mm (grid period). The pulsed voltage source amplitude was set to U0 = 140 kV, and the pulse duration at FWHM was about 1 ns while the voltage pulse rise time was equal to 0.3 ns. For the simplicity reasons real three-dimensional diode connected to axisymmetric transmission line in our model is considered to be a two-dimensional semiperiodic in x-axis structure ( Figure 2). Such a simplification is entirely permissible, since the radial size of the diode is much larger than the mesh period of the cathode grid. Reducing the dimensions of the problem greatly simplifies the calculations, since due to the periodicity it is possible to impose a periodic boundary conditions along x-axis at the boundaries 1-4 and 2-3 (in Figure 2), thereby considering only the 1-2-3-4 area in the simulation. This proposed simplification is reinforced by specific experimental measurement conditions in which the runaway electron beam aperture is usually narrowed by using the adjustable diaphragm between the foil and the current collector [6].

Plasma dynamics model
In order to compute discharge plasma formation and obtain instant discharge characteristics, we use the twodimensional drift-diffusion moment macroscopic (liquid) model [8] including a number of electrons and ions continuity equations coupled to Poisson's equation accounting the electrostatic field self-consistently. Since we consider nanosecond and sub-nanosecond atmospheric pressure discharges, only impact ionization reactions have to be implemented in the model to take into the account the production of single charged ions. This reduces the number of equations, so the plasma dynamics in such case is governed by the following nonlinear partial differential equations system where ni and ne are the ion and electron species densities, Ri is the impact ionization reaction source term, Rph is the photoionization reaction production rate,  and D are mobilities and diffusion coefficients, respectively, E  is the electric field, is the electron energy density,  and R are the energy loss during inelastic collisions and the corresponding reaction rate, respectively, t is the time variable, 0 is the vacuum permittivity, q is the electron charge. The generalized boundary conditions for the particle fluxes and the electron energy density on solid electrodes walls are where n  is the normal vector to the wall surface, Te and Ti are the electron and ion temperatures, respectively, me and mi are the electron and ion masses,  = 0.1 is the secondary electron emission coefficient by ion impact, kb is the Boltzmann constant. Electron and ion fluxes are given by the following expressions Since the fast (nanosecond) discharge is investigated, we use a simplified plasma chemistry mechanism accounting for only impact ionization of neutral species (atoms or molecules depending on the gas medium). For reactions involving electron collisions with heavy neutral species, the reaction rate coefficient   e k T in terms of electron temperature 2 3 e T   was obtained using the two-term Boltzmann equation numerical solver BOLSIG+ [9]. The electron mobility dependence on the electron temperature was also calculated using the BOLSIG+. We took IST-LISBON cross-sections database as the source of cross-section data for the BOLSIG+ numerical solver [10]. In this case, the electron impact reaction source term is given by the expression Ri = k(Te)Nnne, where Nn is neutral gas atoms total number of density.
Along with the role of the electron impact reactions, the extremely important contribution of the gas photoionization during the discharge propagation is also well known from the early investigations [11]. In this paper, we simulate the discharge in atmospheric pressure nitrogen that contains less than 1 % of oxygen admixture. Such an addition of oxygen does not change drift-diffusion characteristics as well as the ionization coefficient noticeably, but the changes of the photoelectron production rate are becoming significant as the characteristic absorption length of ionizing radiation increases considerably. Therefore, the role of photoionization processes in the discharge development have to be taken into the account accurately. The model we use is based on the assumption that the major contribution to the photoionization rate gives the radiation in the spectral range 980-1025 Å, where the absorption of radiation by nitrogen can be omitted. The wavelength of 1025 Å is the photoionization threshold of molecular oxygen. Below the 980 Å value, the radiation is strongly absorbed by nitrogen and provides a minor contribution to the photoelectrons production. Taking into account the fine structure of the oxygen absorption spectrum, in the integral over the wavelength range 980-1025 Å in the general expression for the rate of photoionization was calculated and adopted for further computations [12].
The rate Rph can be calculated by either using of the integral or the differential formulations. In the former case, when the integral model is utilized it involves operations with huge matrices and volume integrations at every single time step that is computationally expensive [12]. To avoid the storage of the extensive number of matrix pattern during the time-dependent study, the differential approach have been implemented [13]. It is based on the numerical solution of radiation transport problem and utilizes number of Helmholtz equations approximating the absorption coefficients of the gas. Within this model, the photoionization rate is expressed by three-term contribution to Rph=S1+S2+S3 obtained by solving three following linear Helmholtz equations 2 2 where j=1,2,3 denotes the equation number, j and Ajare constant sets defined for three-term differential approach according to [13], pO2 is the partial pressure of molecular oxygen in air (pO2 = 150 Torr at atmospheric pressure), p is the gas pressure (p = 760 Torr at normal conditions), pq = 30 Torr is the quenching pressure. The set of equations (1)-(4) is complemented by zero Dirichlet boundary conditions Sj = 0 at electrodes and periodic boundary conditions for 1-4 and 2-3 borders at Figure 2.
It is important that we do not use nonphysical and therefore unreliable "stabilization" methods for the numerical solution in our model, e.g. like setting the high level of initial gap pre-ionization (near to 10 10 cm -3 ), implementing an external intensive volume ionization source or using of artificial diffusion terms for electrons.
The implementing of such auxiliary conditions leads to the numerical distortion of the discharge channel shape, its propagation velocity and decreasing the peak of an electric field value. It is evident that a sufficiently high level of the initial pre-ionization causes the discharge parameters to change strongly without certain physical reasons. We base on the fact that the typical initial ionization level in pure nitrogen does not exceeds 10 3 cm −3 distributed uniformly in the interelectrode space.
The calculations also take into account the contribution of field emission from the electrodes surfaces due to the increase the electrostatic field strength near the possible emission centers formed by the natural roughness of the metallic surfaces. As the flux density of field-emission electrons we use the Fowler-Nordheim current density formula from the original paper [14], which contributes to the condition of the boundary electron flux where n  is the normal vector to the electrode boundary, 0 = 2.12510 30 1/(cm 2 c) and B = 6.5410 8 V/cm are constants,  is the dimensionless surface roughness factor, E is the electric field norm. To vary the cathode autoelectron emission properties, we change the value of the surface roughness factor  in the range of 4-20.
The numerical solution of the equations system (1)-(4) is performed using the finite element method in the time domain. The linear finite elements are used in discretization of the equations of the mathematical model. We also use smoothing of boundary fluxes for the Helmholtz and Poisson equations solution. In order to solve the transport equations (1) coupled to the system of stationary Helmholtz and Poisson equations the segregated time-dependent solver was implemented. The relative accuracy of the time-dependent solver was set not more than 0.0001. We use the triangular-shaped finite elements allowing grid refinement in the rectangular region close to the wire cathode electrode. The maximum mesh size was set to 0.02 mm, while the minimum size was 0.0075 mm in the entire model. Nearly 52000 elements were employed and the number of degrees of freedom (DOF) was about 186 thousands. Since we have been simulated the electric discharge in a periodic structure restricted to the area of 1-2-3-4 ( Figure 2), so in order to verify the correctness of the imposed boundary conditions, we also performed our simulation for the computational geometry containing two and four 1-2-3-4 computational cells in a row. Also the reduction of the computational grid size to half and quarter of that of the base case, yielded output distribution profiles identical to the base case, implying spatio-temporal convergence of the problem.

Analysis and modelling
Since we base on the experimental papers [5,6], in our computations the variation of the distance between plane cathode and cathode grid does not affect the discharge dynamics. The diode switching characteristics (discharge current and potential drop) stay unchanged under the variation of l = 3-12 mm due to the discharge current mainly flows between grid and anode. Below we study two operation modes of the diode with plane-grid cathode system: the case when the anode pulse polarity is positive, and the inverse case of the negative polarity. For both the cases of anode polarity these characteristics are similar, the only difference is that with positive anode polarity the discharge passes into the switching stage faster due to the electron velocity coincides in sign with the direction of the ionization wave motion. Figure 3 shows the current and voltage time-profiles for the case of positive anode polarity. The voltage pulse across the gap in the presence of ballast load only in the electric circuit always has a subnanosecond duration (see Figure 3). Following the discharge current profiles, the diode switching time point corresponding to the electric breakdown is at approximately t = 0.35 ns when the ionization wave from the grid cathode reaches anode and the voltage at the gap significantly falls down.  The discharge development for the positive anode polarity case is demonstrated in Figure 4. On the surface plots, one can see how the electron density profiles progressing towards the anode plane electrode (at the top) until the commutation stage begins. This case represents a subnanosecond negative streamer dynamics in strong electric field.
Due to the amplification of the electric field near the wires surfaces, the impact ionization locally increases the plasma density up to the values of 10 14 cm -3 and even greater near the initial time point during about ~100 ps. After that, a sharp ionization wave propagates into a nonionized gas between grid and anode, leaving a nonequilibrium plasma behind it. In this case, when the anode polarity is positive, the ionization wave front moves towards to the direction of electron drift representing classical negative streamer dynamics. It starts when the physical conditions approach so-called Lozansky-Firsov limit of "ideal conductivity" [15]. As it could be also seen from Figure 4 the velocity of ionization wave increases in proportion to its length in full accordance to Lozansky-Firsov theory [15] as well. Due to the periodicity of the plane-grid cathode system, the ionization wave channel (negative streamer) significantly expands in the transverse direction. At the end of the voltage pulse, a dense plasma completely fills the space between the cathode grid and the anode. This plasma does not penetrate under the grid, but due to the amplification of the electric field between the cathode grid and the cathode, a localized ionized region is also generated near the surface of the wire facing plane cathode, where the plasma density reaches 10 13 cm -3 . Plasma layer below the grid does not move towards the plane cathode, and its density monotonically decreases to it down to 10 5 cm -3 . Figure 5 shows the dynamics of the longitudinal electric field distribution on the symmetry axis of the 1-2-3-4 computational cell in the gap above cathode grid (top) and below it (bottom) before the diode commutation phase starts. It is seen that between grid and anode, the electric field resembles a structure of traveling wave with the variable amplitude. Since the ionization wave (streamer) moves towards anode, the electric field screening occurs inside the streamer body according to the Lozansky-Firsov criterion [15]. In the space between grid and plane cathode the electric field distribution also has strongly non-uniform structure. Despite the fact that electric field screening effect is also exists near the wires, the maximum of the electric field is not in motion during the voltage pulse.
If the anode voltage is negative, we obtain similar voltage and current time-profiles as we have in the previous case. In terms of quantity, the voltage across the gap at the diode switching time point is slightly higher. This reflects a somewhat slow ionization wave headway, supported by photoionization processes instead of drift. The Figure 6 shows current and voltage timeprofiles for this case. As it follows from the discharge current profile, the time point of the switching phase onset coincides with the previous calculation, while the peak voltage across the gap is rather higher (80 kV versus 70 kV in the previous case). It affects the ionization wave dynamics (see Figure 7), where the streamer moves towards anode with the same variable velocity as in [15], but starting later. It takes place due to the need of "plasma spot" formation on the wire surface in the absence of electron drift to the anode plane in a reversed anode polarity.
Here the photoionization is only significant factor leading to initiation and further development of the ionization wave towards anode. If we eliminate the calculation of the photoionization source (4), the "plasma spot" at grid also forms, but the discharge is not developing, so the preionization grows further exceeding the reasonable physical limit equal to the density of neutral gas molecules.  To compare the electric field dynamics on the symmetry axis of the problem, Figure 8 shows the longitudinal electric field distribution as it was shown in Figure 5. It is clearly seen that, despite of the different physical nature of positive and negative streamers, the general regularity of the ionization wave dynamics for both polarities persists. The electric field also represents a traveling wave moving from cathode grid to anode. The region under the grid is partially filled with plasma, locally disrupting the general electrical screening provided by the equipotential grid electrode placement.

Conclusions
As a result of the time-dependent high-pressure gas discharge simulation for the diode with plane grid cathode system [5], the detailed dynamics of streamer discharge in strong electric fields was studied for positive and negative polarities of the voltage source pulse. In particular, it was shown that ionization wave front arises at the surface of the grid facing the anode, and it propagates toward the anode, regardless of its voltage sign. It was also shown that, in the absence of external volumetric ionization sources or the specific initial preionized state, the discharge plasma does not penetrate into the quasi-potential space between grid and plane cathodes. In this region plasma is generated simultaneously with the process of streamer discharge formation. This process can contribute to the emergence of a runaway electron flow directed to the cathode. It is assumed the subject of further research within the framework of the hybrid kinetic approach [7].