3D PIC-MCC simulation of corona discharge in needle-plate electrode with external circuit

This paper investigates the external circuit effect in corona discharge by means of numerical simulation, while the effect is usually disregarded in previous plasma simulations. A general numerical solution method of serial external circuit, which can compute surface charges on any shape of electrodes adaptively and solve electrode potential and circuit current simultaneously, is developed to support the investigation, due to lack of theoretical solution in irregular electrode systems. Brief description and validation of the scheme are presented in section . In general, the accuracy of numerical solution improves with mesh refinement in orthogonal mesh, due to stair step approximation. Simulations are performed in our self-developed 3D PIC-MCC (particle-in-cell, Monte Carlo collision) code with photoionization and external circuit. The simulation model including recently developed photoionization model is described briefly in section . Then the simulation results in a needle-plate electrode in air, such as spatiotemporal evolution of corona discharge and waveforms of electrode voltage and circuit current, are presented. Physical mechanisms of formation of branched corona discharge and self-consistent adjustment of electrode voltage and circuit current are discussed. The waveforms of voltage and current in this simulation are generally consistent with experimental result. The external circuit has negative feedback effect on the discharge, and can adjusts electrode voltage self-consistently. The application of external circuit improves the applicability of particle simulation method to the simulation of discharge phenomenon in actual devices.


Introduction
Low-temperature plasma discharges are of interest both for fundamental research and for technical applications [1] such as plasma biomedicine [2,3], environmental treatment [4,5], flow control [6,7], material processing [8,9], ignition and assisted combustion [10,11]. Numerical modeling and simulation, which have provided considerable insight into the nature of the discharges [12][13][14][15][16], have already become important and essential methods in research of plasma discharges, due to the limited measurement methods in experiments.
While experimental plasma discharges are always triggered and maintained via an electric circuit (comprising of resistance, capacitance, and inductance) [17], modeling and simulations of such discharges usually disregard external circuit and assume a certain voltage pulse shape as input [18]. The addition of an external circuit usually acts as current-limiting device in discharge models. Without a current-limiting device, the cathode current may converge to a wrong volt-ampere regime, which makes it impossible to validate the calculated cathode current with experimental values and may leads to an incorrect prediction of the plasma characteristics [19].
Until now, external circuit has already been coupled with numerical simulation [20,21]. In 1989, Lawson [20] presented a comprehensive review of the considerations involved in a bounded one-dimensional electrostatic plasma simulation code. The analysis also included external circuit elements. However, the scheme presented by Lawson was first-order accurate and did not solve the internal field and the external circuit equations simultaneously [22]. In 1993, Verboncoeur [23] improved upon Lawson's method. The external circuits were solved with the internal fields simultaneously, that makes the circuit plus plasma secondorder accurate and self-consistent. In 1997, Vahedi [22] extended the one-dimensional coupled circuit model to two dimensions. And later, Thiagarajan [24], Bultinck [19], Derzsi [18], Farouk [17] and Lin [25] had performed 2D/3D simulations of plasma discharges with external circuit. However, the external circuit models in above works can be hardly applied to simulations in irregular electrode geometries directly. Actually, most electrode geometries in discharge experiments, such as needle-plate or needle-needle, are irregular, in which it is hard to get electric field strength with analytical theory.
In this paper, we intend to investigate the external circuit effect in corona discharge in air in a needle-plate electrode geometry using numerical simulation. A general numerical solution method of serial external circuit, which can compute surface charges on any shape of electrodes adaptively and solve electrode potential and circuit current simultaneously, is developed to support the investigation. Brief Description and verification of the numerical solution method are presented in section 2. Simulations are performed with a 3D PIC-MCC code, which is developed based on the 2.5D PIC code [26] with c ++ language. Since the needle electrode has a positive voltage, positive corona discharges are considered in the simulation. Thus, photoionization, the most important nonlocal source of electrons in nitrogen/oxygen mixtures [16], is taken into account. Important components of the simulation model are described briefly in section 3. And then simulation results including spatiotemporal evolution of corona discharge and waveforms of electrode voltage and circuit current are presented and discussed.

Numerical solution of a series external circuit
In the 3D PIC-MCC code, we place a cubic mesh over the computational domain. The particle charge densities are weighted to the grid points; the electric fields are assumed to distribute along the cell edges. A cell is the smallest simulation unit and the dielectric property in it is uniform. Curved boundaries are usually approached by stair step approximation, as shown in figure 1. Cells are labelled as either inside a conductor or not.

Numerical calculation of surface charge
Rather than analytical calculation, surface charge on the boundary in this scheme is considered as the sum of surface charges weighted at each mesh point that is on the boundary surface.
Before the calculation, we shall find out all the grid points that are on the conductor surface. Then based on the work of Vahedi [22], as shown in figure 2, Gauss' law at grid  is the dielectric constant at grid point (i+1/2, j, k); Dx, Dy and Dz are grid sizes in x, y and z directions respectively; Q ijk is the total charge at grid point (i, j, k), and is simply the sum of all particle charges that weighted to grid point (i, j, k) in PIC codes.
If the grid point (i, j, k) is on a conductor surface, as shown in figure 3, Q ijk can be written as where r ijk is the free charge density; DV ijk is the volume associated with the grid point (i, j, k) where free charge exists; s ijk is the surface charge density; DA ijk is the area of a physical boundary at (i, j, k) where surface charge density can accumulate, Q sijk is the surface charge of the conductor at (i, j, k).   x y and 0 at each grid point on the conductor surface respectively. If the grid point is on the boundary along the x direction, then S y and S z are halved, and so do along the other two directions. As for DV, it will add D D D x y z 8 for each surrounded cell that is free space.
In the simulation, electric fields can be written as

Coupling to external circuit
The external circuit considered in the paper is a series RLC circuit, connected to a discharge gap, as shown in figure 4. The general circuit equation is obtained by applying Kirchhoff's law around the circuit where L, R, C represents the inductance, resistance, capacitance respectively, V(t) is a constant or time-varying voltage, f t ( ) is the potential difference between electrodes, Q c is Figure 3. The Gauss' theorem on a box at (i, j, k), the grid point is on a conductor surface, the shaded region is a perfect conductor, the blank region is vacuum Figure 4. A series external circuit connected to a discharge gap.
capacitor charge. The central difference form of equation (10) is The surface charge Q s can be obtained from the Kirchhoff's current loop law where, I is the external circuit current, Q abs is the total charge absorbed by the electrode. The discrete finite differenced form of equation (12) can be expressed as Applying equations (6)- (9) to the equation (13), we obtain

Verification and validation
Because of parameters such as capacitance, electric field strength, surface area and surface charge have mature theoretical solutions in a parallel-plate electrode system, simulations are performed to demonstrate the validity of the scheme in a parallel-plate electrode system, as shown in figure 5. The length and width of the plate electrodes are both 2 mm and the thickness is 0.02 mm. The gap distance between electrodes is 4 mm.  The surface charge calculated in simulation is actually affected by mesh generation, as shown in appendix A. With proper mesh generation, the calculated result in simulation is which agrees well with theoretical solution, the error between them is less than 1e-5.

Electrode voltage and circuit current.
It is assumed that the resistance in the external circuit is 100 Ω, the inductance is 0.506 μH, the capacitance is 5 pF, the amplitude of the  applied voltage is set to be −2000 V and the frequency is 5.0×10 8 Hz.
The capacitance of the parallel-plate electrode is The total equivalent impedance of the circuit is The impedance angle is Then, the circuit current can be calculated as Meanwhile, the waveform of voltage and current in the simulation are plotted in figures 7 and 8 respectively. The electrode peak voltage is about 2087 V, slightly larger than the applied peak voltage. The current phase is about 0.5π ahead of the voltage phase, and the peak value is about 0.0578 A. The simulation results agree well with the theoretical results, and the error is about 5‰.
2.3.4. Coaxial electrode. The calculation precision of surface charge in an irregular electrode system also needs to be confirmed, although it is satisfactory in the parallel-plate electrode system. The model used in the test is a coaxial electrode placed in vacuum. The inner radius R 1 is 0.5 mm, the outer radius R 2 is 1.0 mm, and the length L is 0.1 mm. The voltage V, 1 V, is applied on the inner electrode. The surface charge Q s on the inner electrode satisfies PIC Simulations are performed to compute the surface charge on the inner electrode. The grid points on the inner electrode surface identified in the simulation are plotted in figure 9. The numerical errors in different cases are plotted in figure 10.
As shown in figure 10, the error between the computed and theoretical value of the surface charge decreases with mesh refinement. Generally, the numerical calculation precision of surface charge in irregular electrode systems is affected by mesh generation. The denser the mesh is, the more   accurate the expression of curved boundary is, the more accurate the electric field on the electrode surface is, and the higher the accuracy of the numerical solution of the external circuit is. However, the corresponding simulation burden increases, and the efficiency decreases.  figure 11, the needle electrode is a cone, the base diameter is 0.216 mm, the top diameter is 0.016 mm, the distance between the needle tip and the plate electrode is 0.4 mm. Both the length and width of the plate electrode are 1 mm.
The plate electrode is grounded and the voltage is applied on the needle electrode. The applied voltage is a trapezoidal wave, the rising edge is 0.5 ns, and the peak value is 2000 V.
The grid points on the needle surface identified in the simulation are plotted in figure 12. The distribution of electric fields at xz cross section in the computation domain, when the needle potential is 1 V is shown in figure 13. 3.1.2. Numerical error and simulation cost. With the consideration of simulation cost of streamer discharge, we have chosen to use a relative coarse mesh to keep the computational efficiency. The mesh size in the simulation is 125×125×126, the cell size is 8×8×4 μm, and the simulation time step is 5.0 e −15 s. Due to lack of analytical solution, the numerical error of surface charge on the needle electrode is fuzzy. Then it is estimated as 8%-9%, according to figure 10. Numerical uncertainty of corona discharge simulation is discussed in appendix B.
Because of massive particles in streamer discharge simulation, a simplified 4-to-2 particle merging method was applied [27] and the 3D PIC-MCC code was parallelized using MPI (message passing interface). The simulation    presented in sections 3.2 and 3.3 was performed on 90 threads (Intel Xeon Platinum 8180) using at most 9×10 7 simulation particles. It almost spent 100 h to complete the simulation.

Collisions and photoionization.
The background gas is air, a 4-to-1 mixture of N 2 and O 2 , and the pressure is 1 atm. The pre-ionization density of air is set as 10 13 m −3 . Plasma considered in the simulation is low temperature and weakly ionized, thus electron-electron, electron-ion, and ion-ion collisions in discharges are ignored, collisions considered and corresponding impact cross sections are listed in table 1. Photoionization is considered as the dominant source of free electrons in air in this simulation [27]. The photoionization model used in the simulation, proposed by us earlier [34], has considered the generation and de-excitation of nitrogen molecules in high excited states, which is closer to the physical practice than the numerical model proposed by Zhelezhyak [31]. The comparison of physical processes considered in the two photoionization models is presented in figure 14. Detailed description and verification of the new photoionization model can be found in [34].

Formation of corona discharge
The initial stage of positive corona discharge in the electrode gap has been simulated. Electron and ion densities at various times are plotted in figure 15, corresponding electric fields at central cross section are plotted in figure 16.
The needle voltage rises when simulation starts. Seed electrons move toward the needle electrode and collide with neutral gas molecules in the non-uniform electric field between the needle-plate electrodes. The strongest spatial electric fields are distributed around the needle tip, as shown in figure 13, and decrease rapidly with the increase of the distance from the needle tip. Therefore, the movements and collisions of electrons in the region around the needle tip are relatively more intense than that far away from the needle tip, so is the ionization rate. Thus, electron avalanches can only be formed and develop in the region around the needle tip. When electron avalanches reach the needle electrode, electrons are absorbed, leaving positive ions distributed around the needle electrode, as shown in figure 15. Thus, a corona layer is gradually formed around the needle electrode. The ionization degree and electric field strength outside the corona layer are relatively low, as shown in figure 16.
Corona discharge is a basic form of gas discharge at initial stage in non-uniform electric fields. Its formation and propagation are related to the electric field, in which space charge effect also plays an important role. Positive ions in the corona weaken the electric field and restrain the local discharge around the needle electrode; but intensify the local electric field and accelerate the formation and propagation of electron avalanches in the peripheral area of the corona layer, as shown in figures 15 and 16. These electron avalanches will eventually blend into the corona layer, making the corona layer expanding outward. Photoionization provides source electrons to sustain this process. During expansion, electron  avalanches may convert into positive streamers and go on propagating toward plate electrode in the region where local electric field intensity is high enough. Then the streamer destabilizes into branched streamer channels eventually, forming a brush discharge, as shown in figure 15. In our previous research [27], we have pointed out that inception, propagation and branching of a streamer are the results of merging secondary avalanches, that move toward the primary avalanches. When the plasma density at the head of secondary avalanche is close to that at the streamer head, branches may emerge after the combination. At 3.5 ns, the streamer channels have not reach plate electrode yet, discharge between electrodes presents intense local breakdown characteristic and nonlinearity.

Waveforms of voltage and current
The current collected by needle electrode in the simulation is plotted in figure 17. the waveforms of voltage and current are plotted in figure 18, The applied voltage rises linearly from zero when the simulation starts, and reaches the peak value (2000 V) at 0.5 ns, then it maintains at the peak value. The needle electrode voltage rises from zero relatively slowly, and reaches the peak value (about 1814 V, a bit less than the peak applied voltage) at about 1 ns. Then the needle electrode voltage decreases slowly and tends to reach a stable value (about 1440 V). The circuit current increases with the increase of voltage at 0-0.5 ns, and is mainly caused by the charge of electrode capacitance, which belongs to displacement current; while the current collected by needle electrode is almost zero. At 0.5-1 ns, the applied voltage reaches the peak value, the growth rate of electrode voltage drops, and the collected current remains zero, so the circuit current drops in this   period. After 1 ns, the collected current increases, the corona layer emerges around the needle electrode and expands outward, the electrode voltage tends to be stable. So, the circuit current can be considered as causing by gas discharge, and tends to be stable at about 5.5 mA. Briels [35] had given the waveforms of voltage and current in positive streamer discharges at atmospheric pressure between the needle-plate electrodes, as shown in figure 19. At 0-11 ns, the voltage increases rapidly, the current first increases and then decreases, and can be called as capacitive current. At 11-85 ns, the voltage is slightly oscillating around 80 V, while the current is oscillating around 40 A. The waveforms of voltage and current in this simulation are basically consistent with that in Briels's experiment: the voltage increases rapidly, and the current first increases and then decreases at 0-1 ns; the voltage and current waveforms tend to be stable at 1-3 ns.
However, if the external circuit is not considered, the voltage waveform in simulation is consistent with the applied voltage waveform, and the current waveform is close to that collected by the needle electrode, and the volt-ampere characteristic in simulation goes wrong. The application of external circuit makes voltage and current waveforms in simulation can be compared with that in experiment, and improves the applicability of particle simulation method to the simulation of discharge phenomenon in actual devices.
In this paper, we deal with electrostatic problems by solving Poisson equation with space charge effect. Magnetic field is unsolved and ignored in this simulation. So is the intrinsic inductance of plasma, which is necessary to characterize in a very precise way the overall inductance of the load to reproduce secondary oscillations of the type shown in figure 19.

Conclusion
The gas discharge in a non-uniform electric field usually presents strong local breakdown characteristic, and the external circuit plays an important role in it. In this paper, the numerical solution of external circuit for 3D PIC simulations is proposed and verified. With the numerical model, electrode voltage and circuit current can be solved self-consistently in a regular or irregular electrode system. The computed results agree well with the theoretical results in the parallel-plate electrode system. The accuracy of the numerical solution depends on the mesh generation in an irregular electrode system: generally, the denser the mesh generation, the higher the accuracy.
3D PIC simulation of positive corona discharge between needle-plate electrodes at atmospheric pressure has been performed with a serial external circuit. The spatiotemporal evolution of corona discharge is obtained. Plasma discharge between electrodes presents intense local breakdown characteristic and nonlinearity. The waveforms of voltage and current at initial stage of corona discharge, consistent with experimental results, is obtained in the simulation. The application of external circuit makes voltage and current waveforms in simulation can be compared with that in experiment, and improves the applicability of particle simulation method to the simulation of discharge phenomenon in actual devices.

Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant No. U1537210).

Appendix A. Effect of mesh generation on surface charge calculation
In the 3D PIC-MCC code, simple conformal approximation are taken into account during mesh generation. We have defined a conformal parameter Frac, which equals 0.5 in this simulation, to do the approximation. If the volume proportion of a material in a cell is larger than Frac, we consider the cell is full filled with the material, else the cell is free space.  Because of the simple mesh conformal approximation, mesh generation may cause computational error in parallelplate geometry, as shown in figure A1. In case (a), the transaction enlarge the gap distance, then the electric field strength is reduced accordingly. In case (b), operations are opposite. With proper mesh generation, the actual thickness of a material equals the one used in the simulation, then the electric field can be solved accurately. According to equation (7), the calculation precision of surface charge improves.

Appendix B. Numerical uncertainty of corona discharge simulation
Several numerical constraints must be met in plasma discharge simulation [36]: (1) Grid spacing must be on the order or smaller than the electron Debye length, l D  (2) Time step must be a fraction of the plasma period, D < t w 0.2 pe (w e = e n m pe e 2 0 ); when above two constraints are satisfied, the Courant-Friedrichs-Lewy condition is fulfilled for particles traveling at the thermal velocity: The number of charged particles per cell, N P , must be large enough (at least several tens) to avoid numerical heating.
In our simulation, the maximum electron density is about 3×10 21 m −3 , electron energies distributed along z direction at 3.5 ns is plotted in figure B1. The maximum value of electron energy is about 100 eV at the moment, and the average value is about 1.42 eV. Thus, the grid spacing and time step in our simulation should meet the constraints m D  x 0.16 m and D <´t 6.3 10 s. 14 The time step used in our simulations is 5.0×10 -15 s, which satisfies the time step constraint. However, as mentioned by Garrigues [36], it is practically impossible to satisfy the severe grid spacing constraint with current computer technology in 3D PIC-MCC simulations. Influence of grid spacing on the numerical uncertainty of plasma discharge simulation is discussed below, while the number of particles per cell is set as 48 according to works of Garrigues [36] and Teunissen [16]. In fact, the particle merging method excludes   over-weighted particles, and the weight threshold is set as 5.0×10 6 . So, the particle number in a cell can be larger than 48, as the discharge proceeds. Additional simulations of corona discharge with different cell size (10×10×5 and 12.5×12.5×6.25 μm) are performed to confirm the numerical uncertainty. Electron and ion densities in these two cases are plotted in figures B2 and B3. It can be found out that the discharge channels become coarser with the increase of grid spacing, by comparing figure 15 with figures B2 and B3. The corona layer around needle electrode almost disappears, when the cell size is 12.5×12.5×6.25 μm, which indicates a relatively large numerical error in the simulation.
To confirm the numerical heating caused by grid spacing, we list the average electron energies at 3.5 ns and the target values of needle voltages and circuit currents in table B1. The effect of numerical heating do increases with the increase of grid spacing, which makes corona discharge in the simulation more intense in coarser mesh, according to table A1. Therefore, relatively refined mesh shall be applied in 3D PIC-MCC simulations of plasma discharges, although it is hardly to satisfy the grid spacing constraint.