A comparison of 3D particle, fluid and hybrid simulations for negative streamers

The Monte Carlo particle model is of PIC-MCC type. It describes the motion and collisions of free electrons in a streamer discharge in air without photoionization. 'Particle in cell' (PIC) means that the electric charge of electrons and ions is mapped to an electric charge density on a numerical grid; this charge density changes the electric potential (from which the electric field is calculated) according to the Poisson equation. 'Monte Carlo collision' (MCC) means that collisions of the free electrons with the neutral background molecules occur randomly, so neutral molecules do not need to be simulated. Ions are treated as immobile. The Monte Carlo procedure and the differential cross-sections for the particle model are described in detail in section 2.1 of [25].

The particles represent single electrons during the initial stages of the simulation, as indicated in the text. However, when the particle number becomes too large for the computer memory, super-particles are introduced that represent several real particles. While particles carry their generic physical distributions and fluctuations of density and energy, the fluctuations of super-particles are unphysically increased and can generate artifacts when fluctuation effects or rare events become important. To avoid this problem, we have developed the spatially hybrid model. The different models are compared in sections 3 and 4.

The classical fluid model for streamers has a long tradition in streamer modeling, much longer than the more microscopic particle model. Originally it was a phenomenological model based on the essential physical mechanisms and conservation laws. It can be traced back at least to the 1930s. The classical fluid model approximates the electron dynamics by a reaction–drift–diffusion equation for the electron density, and the reaction and transport coefficients are assumed to depend on the local electric field, in the so-called 'local field approximation'. The model is completed with the reaction equation for the ions and with the Poisson equation for the electric potential.

In order to approximate the particle dynamics well, the coefficients in the fluid model should be derived from the particle model. This can be performed either through averaging over the Boltzmann equation for the electron distribution in configuration space, or through evaluating swarm simulations in a Monte Carlo particle model; in a swarm simulation the evolution of an approximately Gaussian electron distribution in a constant electric field is traced by a particle model. For this paper the transport coefficients and the reaction rates for the fluid model are derived by fitting the fluid coefficients to the swarm dynamics in a Monte Carlo particle model, as evaluated in [25]. (Note that [25] contains some corrections to the reaction and transport rates described in [10, 33].) The coefficients in [25] were derived up to a field of 250 kV cm⁻¹ in air at standard temperature and pressure, and for stronger electric fields fit formulae from the same paper are used.

Furthermore, in a reactive plasma one needs to distinguish between bulk and flux coefficients [35].While bulk coefficients characterize the dynamics of a swarm as a whole including its reactions, flux coefficients characterize the dynamics of individual electrons within a swarm. Robson et al [35] express the general opinion that bulk data should not be used in low-temperature plasma simulations; see also figure 2 in [35].

However, essential physics is missing in the classical or minimal model as we found in [33], and as will be discussed much more extensively and systematically in a forthcoming paper [36]. For example, in strong electric fields, the ionization rate cannot be computed accurately from just the local electric field and the local electron density. Therefore in [33] we extended the generation term in the fluid model with a gradient expansion term while in [36] a higher-order model is derived in a systematic manner by averaging over more moments of the Boltzmann equation. An important conclusion is that accurate results cannot be expected from the classical fluid model, independently of whether flux or bulk coefficients are used, as the functional form of the equations is not sufficiently general.

For comparing the classical fluid model with the other models, we have chosen to use bulk coefficients. From studies of planar (1D) fronts, we know that bulk coefficients in the classical fluid model will lead to approximately correct front velocities. This is the case because the leading edge of an ionization front, which pulls the front along, propagates under swarm-like conditions; for details, we refer to [10] or for a full mathematical analysis of pulled front dynamics to [11]. A particle swarm should be parametrized with bulk coefficients in the classical model (recall that the bulk coefficients are constructed to give agreement for particle swarms). Therefore, the streamer ionization front should also be modeled with bulk coefficients in the classical model. The front velocity is then close to that of a particle model in the same electric field, although the ionization density behind the front is too low [10]. If flux coefficients would be used, the ionization density behind the front would be a little higher. (Both sets of coefficients produce the same amount of ionization per unit time in a given electric field. The flux mobility is typically lower, therefore a larger amount of electrons per unit length is produced, and the electron density is higher.) But the front velocity would be significantly too low with flux coefficients, and therefore we use bulk coefficients. We stress again that accurate results cannot be expected from the classical fluid model with either set of coefficients.

The ionization term in the classical fluid model for streamers is calculated in the local field and local density approximation, but the comparison with particle models shows that the ionization densities in the streamer interior are too low behind a planar front in a fixed electric field [10, 33, 37]. This is because the mean electron energy varies even within a swarm in a constant electric field: at the front edge of the swarm, the electrons have higher energies and are more likely to ionize the neutrals, while the electrons at the back end of the swarm are slower on average and less likely to ionize. By including the first term of a gradient expansion in the electron density in the impact ionization rate, both particle swarms and planar ionization fronts are approximated well [33]. The extra term was derived in a model with flux coefficients; therefore we need to use flux coefficients for consistency with the model derivation.

The hybrid model connects the particle model with the extended fluid model through a moving model interface with a buffer zone, as described extensively in [25] and briefly recalled in the introduction. When the flux of electrons across the interface between particle and fluid model is calculated, the same definition of coefficients should be used on both sides in order to be physically consistent; otherwise the physical inconsistency becomes visible in the form of a discontinuity of the electron density at the model interface [33]. For further details and the numerical implementation, we refer to [10, 25, 33]. The important feature of the spatially hybrid model is that it follows the particle dynamics only in the dynamically relevant region, and that it therefore can continue to track the single-electron fluctuations much longer than the (super-)particle model.

We simulate the evolution of a negative streamer in air without photoionization, at standard temperature and pressure. The streamer propagates in a background field of −100 kV cm⁻¹, or about 400 Td; this field is well above the breakdown value. The simulation volume is 1.17 mm long in the z-direction parallel to the electric field and extends up to ±0.29 mm outwards from the axis in the x- and y-direction where homogeneous Neumann boundary conditions are applied to the electric potential.

The initial distribution of electrons and ions is generated in the following manner: first 500 electrons and ions are placed at a distance of 0.115 mm from the cathode on the z-axis and followed by the particle model for 60 ps. At this time, there are about 2500 electrons and ions with spatial distributions close to a Gaussian; the electrons are at this time all in the interval of 0.120–0.136 mm from the cathode. These electron and ion distributions are then used as an initial condition for the simulations in all four models; for the fluid models, the swarm is mapped to densities on the numerical grid.

The models were already described in the previous section, and references to more detailed discussions were given there as well. Electric field and electron and ion densities are calculated on a uniform grid of 256 × 256 × 512 points with Δx = Δy = Δz = 2.3 µm, using the numerical schemes described in section 2.2 of [25]. The time step is Δt = 0.3 ps. The Poisson equation for the electric field is solved in all models at each time step with the same fast elliptic solver FISHPACK [38].

In both particle and hybrid simulation, the particle model with single electron is used in the early stages, until the number of electrons reaches 2 × 10⁷; this occurs at about 0.46 ns. At this time, the particle model switches to super-particles, while the hybrid model switches to the full hybrid scheme: the fluid model is applied inside the streamer channel where the electric field is less than 0.95 E_b or where the electron density is larger than 0.7 n_e,max, and the particle model in the remaining part of space; here E_b stands for the background field and n_e,max is the maximal electron density in the complete simulation volume. When the hybrid model is activated, electrons inside the streamer channel are removed from the particle list, and the particle model continues to trace all single electrons at the streamer head. In the particle model, on the other hand, super-particles are introduced at the time 0.46 ns by removing at random every second electron, and by doubling the weight of the remaining particles.

Figure 1 shows the evolution of the electron density in the streamer in seven stages, from time 0.72 up to 0.9 ns, with time increments of 0.03 ns, and figure 2 shows the electric charge densities and the electric fields at 0.72 and at 0.9 ns. The rows in figures 1 and 2 present from top to bottom: the classical fluid model, the extended fluid model, the particle model and the hybrid model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows the electron density and the electric field on the z-axis for the four different models at 0.72 and 0.9 ns.

Zoom In Zoom Out Reset image size

Fluctuation and destabilization effects can be seen more clearly in figures 4–6 that zoom into the propagating streamer heads. The figures show electron density, negative space charge density and electric field for the same seven time steps as in figure 1 for the extended fluid model, the particle model and the hybrid model. The classical fluid model is not included since figures 1–3 demonstrate clearly that it does not approximate the particle dynamics well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At the earlier stage of 0.72 ns, figures 1 and 2 show that a streamer has emerged and grown to about the same length in all models, and that it has approximately the same radius and field enhancement at the tip. The streamer in the classical fluid model (upper row) has stayed a little behind, and electron and charge density are lower, though the field enhancement is still similar. The propagation differences can be seen more clearly in figure 3, which shows the electron density and electric field on the z-axis at 0.72 and 0.9 ns. At 0.72 ns, the profiles of the extended fluid, particle and hybrid model are about the same, but the classical fluid model has a lower field enhancement, lower electron density and shorter propagation length. That the streamer in the classical fluid model grows more slowly both in space and in electron density, while the other models have comparable results, is also reflected in the total number of electrons: it is (5.3 ± 0.2) × 10⁸ in the particle, hybrid and extended fluid model (more precisely 5.1, 5.6 and 5.3 × 10⁸), while it is only 2.9 × 10⁸ in the classical fluid simulation. We will discuss the dynamics of the classical fluid approximation in more detail in section 4.1.

At this stage, the particle model uses 1.6 × 10⁷ super-particles with a weight of 32 real electrons, while the hybrid model follows 2 × 10⁷ real electrons and leaves the rest to the fluid region. The super-particles in the particle model already create visible fluctuations of the space charge density, as discussed earlier in [32]. The fluctuations of the local field create numerical heating, and this effect increases as time evolves. Such numerical artifacts can be somewhat suppressed if the super-particles are formed adaptively using particle coalescence techniques [14, 39–41]. But as the number of particles increases, the increased fluctuations will affect the simulations, especially for negative streamers where perturbations in the electron density can grow as they move outwards. The hybrid model does not suffer from such artifacts.

During the further evolution up to time 0.9 ns, the streamer ionization front destabilizes in three of the four models, but in characteristically different manners.

The streamer in the classical fluid model is destabilizing into off-axis branches, which are quite symmetric (as we expect in the deterministic fluid model, and as we have seen previously in the simulations of Montijn et al [5]). The actual branching can also be seen in the plot for time 0.9 ns in figure 3: the electron density on the z-axis in the classical fluid model starts to decrease for z > 0.69 mm. As the ionization density is determined by the electric field at the front at the moment when it passed that particular position, the electric field at the front increases until the position 0.69 mm is reached, and decreases thereafter on the axis as the lateral protrusions grow and screen the electric field on the axis.

The streamer in the extended fluid model propagates in a stable manner until the last time step 0.9 ns. The off-axis branching of the classical fluid model is suppressed by the higher ionization rates on the axis in the extended fluid model, which are due to the gradient correction in the ionization term. However, it cannot be excluded a priori that the streamer later destabilizes along a different mode.

The streamers in the particle and in the hybrid model both clearly show density fluctuation effects, both at the front and in the interior charge density. The fluctuations in the particle model are unphysical due to super-particle artifacts while the fluctuations in the hybrid model stay physical. The fluctuations at the front destabilize the streamer into several branches in both models, but the streamer tips close to the axis keep the strongest field enhancement and screen new branches up to the end of the simulation at 0.9 ns. It should be mentioned here that the hybrid model operates with real particles up to 0.69 ns, but has introduced super-particles of weight 8 at 0.9 ns. The weight of the super-particles in the particle model at 0.9 ns is 16 times higher, namely 128.

The extended fluid model and the particle model agree very well in propagation velocity, field enhancement and ionization density, up to the large super-particle fluctuations in the particle model. The hybrid model has a similar ionization density and electric field profile, but is ahead of the other models. At 0.9 ns, there are 3.9, 4.6 and 3.5 × 10⁹ electrons in the particle, hybrid and extended fluid model, and only 1.3 × 10⁹ in the classical fluid model.

To obtain the results presented here, the computing time for the fluid simulations was about 1 week, for the hybrid simulation about 1.5 weeks, and for the particle simulation about 2 weeks (all on an Intel Q6600 2.4 GHz quadcore processor). The cost of solving the Poisson equation at every time step dominated the total computational cost, therefore computing times are similar for the different models. All the simulations ran sequentially on a single core.

Our results show that the streamers in the classical fluid model develop lower velocities, field enhancement and ionization densities than in the other models; the model clearly approximates the microscopic dynamics quite badly, as also found previously in [10]. This is the case even though the transport and reaction coefficients were derived from swarm simulations in the particle model for consistency, and though bulk coefficients were used. As the fields do not exceed 250 kV cm⁻¹ in the simulations with the classical model, the fluid coefficients were only used in the parameter range in which they were actually derived in [25]. By construction, the classical fluid model with bulk coefficients models electron swarms in a constant electric field well, and earlier numerical studies as well analytical arguments have shown that also the velocity of a planar front in a fixed electric field is well approximated [10, 33]. However, the ionization density behind a planar front in a fixed field is too low in the classical fluid model with bulk coefficients when the maximal field exceeds 50 kV cm⁻¹ [10]; this is always the case in the present calculations.

In the 3D simulations, the deviation from the other models is larger than in 1D [10]. We argue that this is because in our 1D front simulations, the electric field ahead of the front is fixed, while in 3D the field falls off with distance and varies in time. As the same field in the front creates a lower ionization density in the classical fluid model, also the conductivity in the streamer channel and the consecutive electric screening are lower than in the other models. Therefore the field enhancement is less in 3D and leaves an even lower ionization density behind as the ionization level depends on the field at the front. The lower field enhancement also explains why the streamer is slower than in the other models.

Choosing flux rather than bulk coefficients had not resolved the discrepancy with other models either, as already discussed in section 2.2. The model with flux coefficients does not reproduce the swarm results in a constant electric field. Furthermore, with flux coefficients the electron mobility had been considerably lower (see figure 3 in [25]); therefore the front had been even slower than with bulk coefficients.

The streamers within the classical fluid model destabilize and branch at about the same time as in hybrid and particle model, but in a more symmetric manner, as the destabilization is not supported by electron density fluctuations. This deterministic branching in a fluid model is well approximated by moving boundary models as studied in [6] and reviewed in [3].

We note that the front destabilization in the present fully three-dimensional simulations of the classical fluid model for negative streamers without photoionization occurs in a very similar manner as in previous high accuracy calculations under the constraint of cylindrical symmetry [5]. Earlier simulations [4] demonstrated the mechanism but suffered from lower numerical resolution. In [5], streamers in the same background field were studied, but with a more ionized initial condition and attached to a planar electrode. In those simulations, the branching instability occurred after 1.09 ns, while here it occurs after 0.81 ns. The somewhat longer evolution time until branching in [5] could be due to the different initial and boundary conditions, or to the less accurate transport and reaction coefficients in [5] or to the symmetry constraint; this is a subject of future research. In any case, the results support the argument given in [42] that the branching time under the symmetry constraint is an approximation and upper bound of the branching time in the fully 3D calculation.

The extended fluid model was constructed to cure the deficiencies of the classical model. It was shown already in [10, 33] that with the extension in the reaction term and with flux coefficients, it approximates the growth and propagation of particle swarms and of planar streamer fronts well, including the ionization density behind an ionization front. It should be noted though that the reaction and transport coefficients are used here for up to 400 kV cm⁻¹ while they were derived only for up to 250 kV cm⁻¹ in [25] and extrapolated to higher field values. The hybrid model uses the fluid coefficients only in the range in which they were derived.

We already discussed above that destabilization into off-axis branches in the extended fluid model is less likely than in the classical fluid model, but we have currently no explanation why branching does not occur at all; possibly this is a mere coincidence and branching does occur at some time after the end of the present simulations at 0.9 ns. The front destabilization in particle and hybrid model occurs at a similar time as in the classical fluid model, but in a different manner: the fastest propagating branch stays close to the axis and screens the other branches that therefore keep staying behind.

The figures show that the extended fluid model approximates the particle and hybrid model well up to the moment of destabilization; after this moment there are characteristic differences due to the density fluctuation effects caused by the discreteness of the electrons. These fluctuations have an unphysical distribution when the particle model needs to use super-particles, and therefore the hybrid simulations should be closer to the true dynamics. The stronger destabilization in the particle model seems to be compensated by the increased noise in the space charge distribution (see figure 5), so that in the end the front moves with essentially the same velocity as in the extended fluid model.

The front positions (see figure 7(a)) and the density and field profiles on the axis (see figure 3) agree very well between the extended fluid model and the particle model even at the latest stages, when the particle model shows strong fluctuation effects; the streamer in the hybrid model is a little faster at the latest stages. This could be due to two different reasons.

• The single-particle fluctuations resolved in the hybrid model cover rare electron runaway effects better. This could create more new avalanches ahead of the front, so that eventually the front jumps forward when the electron density within the avalanches has increased sufficiently. In this case the front would be moving faster than expected from reaction, drift and diffusion in the local electric field.
• Alternatively, the single particle fluctuations create more branching and thinner streamers with more field enhancement. The front would then propagate faster because the local field is higher, and not because electrons run away.

Zoom In Zoom Out Reset image size

The maximal electric field as a function of time is plotted in figure 7(b). Comparison with the front position in figure 7(a) already points to the second statement: where the hybrid model is ahead of the other models, the maximal field enhancement is higher as well. The question is further analyzed in figure 8. For planar fronts in a slowly varying electric field E and with a sufficiently rapidly decay of the electron density ahead, the front velocity in the classical fluid model is given by [43]

$$\begin{equation} v^*=\mu_{\rm e}|E|+2\sqrt{D_{\rm e} \mu_{\rm e}|E|\alpha}, \end{equation} \tag{ 1 }$$

where μ_e is the electron mobility, D_e is the electron diffusion constant and α is the effective ionization coefficient. We now use this equation for all models, not only for the classical fluid model. We insert the maximal electric fields E(t) on the z-axis of the respective models into this equation and evaluate it with our flux coefficients for μ_e(E), D_e(E) and α(E). As the front velocity and the maximal field in the hybrid and particle model fluctuate heavily during the late stages, we do not compare velocities, but the resulting front displacements of the models in figure 8. Here the lines indicate the simulation results and the extended symbols the front displacement as predicted by equation (1).

Zoom In Zoom Out Reset image size

Up to time 0.6 ns and a maximal field of about 200 kV cm⁻¹, the models agree very well with the analytical approximation of (1). During the further evolution, the classical fluid model shows almost no deviation from (1). (We should note that this is somewhat accidental, as the classical fluid model uses bulk coefficients and the approximation flux coefficients.) For extended fluid model, particle model and hybrid model, the deviations follow a similar trend at later times: the simulation models are always a little ahead of the analytical approximations, and in both the hybrid model is ahead of the particle and extended fluid model.

One can conclude (i) that the hybrid model is ahead of the others because the field enhancement is higher, and (ii) that equation (1) is a reasonable approximation of the front velocity in all models, but somewhat too low at higher fields. This might be related to the fact that the transport and reaction coefficients were extrapolated from 250 kV cm⁻¹ to higher fields.

Some high energy electrons can 'run away' from the front in the hybrid and in the particle model. The first electrons with energies above a typical runaway threshold of 200 eV appear at time 0.585 ns in the hybrid model, when the maximal electric field reaches 190 kV cm⁻¹, in agreement with the results in [25]. In the particle model the first electrons pass the 200 eV threshold somewhat later, at time 0.65 ns. At the end of the simulation the energies of some electrons exceed 1 keV both in the hybrid and in the particle model. The role of runaway electrons for the front speed will be investigated further in future work, but it does not seem to have a significant influence on the results presented here, according to our velocity analysis above.

Finally, it should be noted that we used the maximal electric field on the axis in our analysis, while after front destabilization at 0.75 ns, the true maximum fluctuates over some near-axis positions.

We have tested four 3D models for negative streamers in air without photoionization in overvolted gaps, and we have found a clear advantage for the hybrid model. It offers a fast and accurate method to model streamers in cases when rare events are significant, like electron runaway or ionization avalanches created by single electrons. Simulations with the hybrid model run faster than with the particle model, with computing times approaching those of the fluid models, and they do not suffer from super-particle artifacts. The extended fluid model is a good approximation up to the moment of front destabilization, but lacks realistic fluctuations. In the classical fluid model propagation speeds and ionization densities and field enhancement are too low.

We have studied short negative streamers in high fields, because this allows us to run physically meaningful 3D simulations with all models, without the need to introduce (adaptive) grid refinement; at this moment, adaptive grid refinement is only available in our fluid model [7], but introducing it across the model boundary in the hybrid model still poses a major challenge. This forced us to exclude photoionization from the model for reasons discussed below in the outlook, and in that sense our present results directly apply to discharges in gases like high-purity nitrogen [26, 27] or the atmosphere of Venus [47] where photoionization is very weak.

Our study lays a basis for future quantitative investigations of electron runaway from streamers and streamer branching, based on the methodology of analysis and numerical models presented here.

As the presented investigations do not include photoionization, the question rises as to how they change when the non-local photoionization mechanism is included for realistic models in atmospheric air. Basically this depends on whether the background field is below or above the breakdown field, i.e. whether the gap is undervolted or overvolted.

Now our choice of system parameters was constrained by the fact that we had to resolve both the space charge structure of the ionization front and the complete system with the developed streamer in 3D without grid refinement. Therefore we chose a high background field to initiate the streamer sufficiently rapidly. This field is well above the breakdown value where impact ionization balances electron attachment, i.e. the gap is overvolted. A consequence is that the results will change substantially if photoionization is included: wherever the non-local photoionization mechanism creates a new electron ion pair, a new local ionization avalanche will appear. The emerging structure was demonstrated with the 3D hybrid model in figure 1 of [44] and in figure 14 of [25], and with Luque's density fluctuation model in figure 2 of [8]. The figures show how eventually the whole space is filled with new avalanches at stochastic locations, suppressing the field enhancement at the streamer tip. Depending on initial conditions, eventually the whole region above the breakdown field will be filled with plasma. This prediction is completely in agreement with experiments where the high-field region around a needle electrode fills up with an ionization cloud that breaks up into streamers only beyond some critical radius [48, 49]. (Of course, this behavior also depends on the initial electron density distribution while streamer development is rather insensitive to this distribution in undervolted gaps [45].)

The onset of branching of positive streamers in air with photoionization in an undervolted gap was recently investigated in [8] with a model accounting specifically for density fluctuations. There it was found that the density fluctuations that are due to the discrete nature of electrons accelerate streamer branching. Similarly, the present results show that the fluctuations of electron densities and energies destabilize the ionization front earlier than without fluctuations, but this does not need to create permanent branching in negative streamers. These observations open up many questions: are negative streamers more self-stabilizing than positive ones as various experiments also seem to indicate? What differences are caused by over- or undervolted gaps or by the presence or absence of photoionization, next to the polarity differences? What role is played by runaway electrons? Our study lays a methodological basis for future studies of these questions.