JOURNAL OF APPLIED PHYSICS 101, 123305 共2007兲

Deviations from the local field approximation in negative streamer heads

Chao Lia兲
Centre for Mathematics and Computer Science (CWI), P.O. Box 94079, 1090 GB Amsterdam, The
Netherlands
W. J. M. Brok and Ute Ebert
Centre for Mathematics and Computer Science (CWI), P.O. Box 94079, 1090 GB Amsterdam, The
Netherlands and Department of Applied Physics, Eindhoven University of Technology, P.O. Box 513, 5600
MB Eindhoven, The Netherlands
J. J. A. M. van der Mullen
Department of Applied Physics, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven,
The Netherlands
共Received 15 February 2007; accepted 4 May 2007; published online 21 June 2007兲
Negative streamer ionization fronts in nitrogen under normal conditions are investigated both in a
particle model and in a fluid model in local field approximation. The parameter functions for the
fluid model are derived from swarm experiments in the particle model. The front structure on the
inner scale is investigated in a one-dimensional setting, allowing reasonable run time and memory
consumption and high numerical accuracy without introducing superparticles. If the reduced electric
field immediately before the front is 艋50 kV/ 共cm bar兲, solutions of fluid and particle model agree
very well. If the field increases up to 200 kV/ 共cm bar兲, the solutions of particle and fluid model
deviate, in particular, the ionization level behind the front becomes up to 60% higher in the particle
model while the velocity is rather insensitive. Particle and fluid model deviate because electrons
with high energies do not yet fully run away from the front, but are somewhat ahead. This leads to
increasing ionization rates in the particle model at the very tip of the front. The energy overshoot of
electrons in the leading edge of the front actually agrees quantitatively with the energy overshoot in
the leading edge of an electron swarm or avalanche in the same electric field. © 2007 American
Institute of Physics. 关DOI: 10.1063/1.2748673兴

I. INTRODUCTION Another recent result supporting the fluid approximation

for streamers was that even streamer branching24–26 can be
Streamers1,2 are growing filaments of weakly ionized understood in terms of an inherent instability of the fully
nonstationary plasma produced by a sharp ionization front deterministic fluid equations.18,26–29 These studies have
that propagates into non-ionized matter. Streamers are used shown that a streamer in nitrogen can reach a state in which
in industrial applications such as lighting3,4 or gas and water the width of the space charge layer that creates the field
purification,5,6 and they occur in natural processes as well enhancement at the streamer tip is much smaller than the
such as lightning7–9 or transient luminous events in the upper streamer diameters; the streamer then can branch spontane-
atmosphere.10 Therefore accurate modeling and simulation of ously due to a Laplacian interfacial instability.
streamers is of high interest. However, despite success and progress of fluid approxi-
Most streamer models 共see, e.g., Refs. 11–20兲 are so-
mations and simulations for streamers, there are three major
called fluid models for the densities of different particle spe-
reasons to reinvestigate the local field approximation.
cies in the discharge. These models build on the assumption
of local equilibrium: transport and reaction coefficients in the 共1兲 Not all streamers are alike. Experiments as well as simu-
continuity equations are functions of local parameters only. lations show that rapidly applied high electric voltages
If this parameter is the local electric field, we refer to this can create streamers that are more than an order of mag-
assumption as the local field approximation. This assumption nitude faster and wider than streamers at lower
is commonly considered to be valid as long as equilibration voltages.30 Whether earlier findings on streamers in
length or time scales are much smaller than the spatial or lower potentials apply to those fast and wide streamers
temporal gradients in the electric field. For the strongly vary- as well has to be investigated.
ing electric fields within a streamer ionization front, the va- 共2兲 The detection of x rays emanating from lightning
lidity of the local field approximation was investigated in strokes31–34 indicates that electrons can gain very high
Refs. 21–23 The general “sentiment” in these studies is that energies within early stages of the lightning event.
the approximation suffices for practical purposes and that Therefore runaway electrons within streamer and leader
more detailed methods tracking the behavior of individual processes should be investigated. Runaway electrons by
particles lead to just minor corrections. definition violate a local approximation.
共3兲 Streamer branching is an inherent instability of a fully
a兲
Electronic mail: [email protected] deterministic fluid model. However, fluctuations of par-

0021-8979/2007/101共12兲/123305/14/$23.00 101, 123305-1 © 2007 American Institute of Physics

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-2 Li et al. J. Appl. Phys. 101, 123305 共2007兲

this outer scale problem concerns only the electric field and
can be dealt with through an inner-outer matching
procedure.38–40 Planar fronts allow us to investigate indi-
vidual particle kinetics and fluctuations within the front and
its specific strong spatiotemporal gradients in a systematic
way and within reasonable computing time.
In this paper, we concentrate on negative streamer fronts
in pure nitrogen under normal conditions. We thoroughly dis-
cuss the case where the reduced electric field at the streamer
tip is 100 kV/ 共cm bar兲, and we summarize results for fields
ranging from 50 to 200 kV/ 共cm bar兲. The paper is organized
as follows. In Sec. II, first our Monte Carlo particle code and
its numerical implementation are described. Then the deriva-
tion of the fluid model is recalled, and the numerical imple-
mentation of the fluid model is summarized. Then swarm or
avalanche experiments in a fixed field are performed in the
FIG. 1. 共Color online兲 The relation between the full streamer problem and
the planar fronts described in this paper. The left picture shows the narrow particle model, the approach of electrons to a steady state
space charge layer surrounding the negative streamer head;18,28,29 the width velocity distribution is investigated, and the parameter func-
of the layer is much smaller than the streamer diameter which creates the tions for the fluid model are generated. This sets the stage for
characteristic field enhancement ahead and field suppression behind the a quantitative comparison of front solutions in particle and
front. The right picture shows a zoom into the inner structure of the space
charge layer with an essentially planar ionization front as treated in this fluid model in Sec. III. Here first the setup of planar front
paper. In the transversal direction, periodic boundary conditions are applied. simulations is described, then the results of the planar front
The simulation can start initially with charge Q evenly distributed in a thin simulations within fluid and particle model and analytical
layer of transversal area A. Q is so large that it screens the field below the
results are presented and compared. The emphasis lies on
charge layer.
front profile, front velocity, and ionization level behind the
front. It will be shown that differences can be attributed to
ticle densities might trigger this instability earlier than the electron kinetics in the leading edge of the front where
they would occur in the fully deterministic fluid model. the electric field does not vary, and that the electron energy
In particular, in the leading edge of an ionization front, distribution there agrees quantitatively with that in the lead-
particle densities are very low and the fluid approxima- ing edge of an ionization avalanche or swarm. Section IV
tion eventually breaks down. As the front velocity of this contains our conclusions on the validity of the fluid approxi-
so-called pulled front17,35 is determined precisely in the mation. An Appendix contains analytical approximations on
leading edge region, single particle dynamics and fluc- the ionization level behind an ionization front.
tuations should be accounted for.

These three observations motivate our present reinvesti- II. SETUP OF PARTICLE MODEL AND FLUID MODEL
gation of the local field approximation for streamers. The IN LOCAL FIELD APPROXIMATION
starting point is a Monte Carlo model for the motion of In this section, we summarize the features of particle and
single free electrons in nitrogen. We note that complete fluid models, their numerical implementation and mutual re-
streamers have been simulated with Monte Carlo particle lation as a basis for the quantitative comparison of ionization
models before;36 however, a drawback of such models is that fronts in particle and fluid model in Sec. III.
the computation time grows with the number of particles and Our starting point is a model that contains all micro-
eventually exceeds the CPU space of any computer. This scopic physical mechanisms that are thought to be relevant
difficulty is counteracted by using superparticles carrying for the propagation of a negative impact ionization front in
charge and mass of many physical particles, but superpar- pure nitrogen. It models the generation and motion of single
ticles in turn create unphysical fluctuations and stochastic free electrons and positive ions in the neutral background
heating. gas. While propagating freely, the electrons follow a deter-
In the present paper, we compare the results of a Monte ministic trajectory according to Newton’s law. The collision
Carlo particle model and a fluid model. We circumvent the of electrons with neutral molecules is treated as a stochastic
problems caused either by a too large particle number or by Monte Carlo process. Because the mobility of the positive
the introduction of superparticles by investigating a small, ions is two orders of magnitude smaller than that of the elec-
essentially one-dimensional section of the ionization front as trons, ions are treated as immobile within the short time
illustrated in Fig. 1. We suppress effects of lateral boundaries scales investigated in this paper. Neutral molecules are as-
by periodic boundary conditions. As the electric field essen- sumed to have a uniform density with a Maxwellian velocity
tially does not deviate from the planar geometry within the distribution. The electron-neutral collisions, including all rel-
region where the particle densities vary rapidly, a planar ion- evant elastic, excitation, and ionization collisions, are treated
ization front37 is a very good approximation of this inner with the Monte Carlo method. Electron-electron or electron-
structure. Of course, a planar front will not incorporate the ion processes as well as density changes of the neutral gas
electric field enhancement caused by a curved front,11,37 but are neglected as the degree of ionization stays below 10−5

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-3 Li et al. J. Appl. Phys. 101, 123305 共2007兲

even at atmospheric pressure.17,26 This well-known model In an elastic collision, the longitudinal scattering angle ␹
will be summarized in Sec. II A. The space charges can and the azimuthal scattering angle ␸ relative to the direction
change the local electric field; this is accounted for by solv- of the incident electron, are given in Ref. 43. In an inelastic
ing the Poisson equation. The particle model gives a very collision, the energy loss of incident electrons has to be taken
detailed and complete description at the expense of signifi- into account, but the scattering angle is calculated in the
cant computational costs where we stress that one particle is same way as for an elastic collision.
one electron and superparticles are not used. In an ionizing collision, energy conservation dictates
If densities are high enough and fields vary slowly in
␧1 + ␧2 = ␧c − ␧ion , 共1兲
space and time, the average behavior of the electrons can be
modeled by a fluid approximation for electron and ion den- where ␧c, ␧1, and ␧2 are the energy of the incident, the scat-
sities whose parameters depend on the local electric field tered, and the ejected electrons, respectively, and ␧ion is the
only. The derivation of the fluid approximation can be for- ionization threshold energy. We use Opal’s empirical fit47 for
malized by taking the zeroth and the first moment of the the distribution of the ejected electron energy,
Boltzmann equation. However, for the practical purpose of
determining mobilities, ionization rates, and diffusion coef-
ficients as a function of the electric field, we directly perform
冋
␧2 = B tan p1 arctan
2B
册
␧c − ␧ion
, 共2兲

swarm experiments with the particle model in a constant where B ⬇ 13 eV in the energy range of interest and p1 is a
electric field. This procedure together with the averaging random number equally distributed between 0 and 1. For the
processes involved are described in Sec. II C. Here also the scattering angles, Boeuf and Marode45 assumed that 共i兲 the
relaxation of an electron swarm to steady state motion and incident, ejected, and scattered electron velocities are copla-
the velocity distribution of steady state motion in a given nar, and 共ii兲 the scattered and ejected electron velocities are
homogeneous field are discussed. perpendicular. These assumptions lead to
␧1 ␧2
A. The Monte Carlo particle model cos2 ␹1 = , cos2 ␹2 = , 共3兲
␧c − ␧ion ␧c − ␧ion
1. Physical processes
where ␹1,2 are the respective scattering angles. The set of
In the particle scheme, at each instant of time t, there is Eqs. 共1兲–共3兲 determines the distribution of energies and scat-
a total number of Ne共t兲 electrons and N p共t兲 ions. The single tering angles of the scattered and the ejected electron in an
electrons are numbered by i = 1 , . . . , Ne共t兲 at time t; they are ionizing collision.
characterized by a position xi共t兲 and a velocity vi共t兲, each
within a continuous three-dimensional vector space. Between
2. Numerical implementation
collisions, electrons are accelerated and advanced according
to the equation of motion. For the positive ions, only their The particle code moves electrons within the applied
position x pj , j = 1 , . . . , N p共t兲 is taken into account while their plus the self-induced field and includes their collisions.
mobility is so low that their velocity can be neglected. The Therefore the numerical calculation consists of three parts:
electric field is determined from the Poisson equation to- the Newtonian electron motion within the field, the field gen-
gether with appropriate boundary conditions. erated by the charged particles, and collisions. At each time
The collisions account for the impact of free electrons on step of length ⌬t, the field is calculated from the charge
neutral nitrogen molecules. As the neutrals are abundant, densities on a lattice with grid size ⌬ᐉ. Then the electrons
their density determines the probability of an electron-neutral move in continuous phase space according to the field, pos-
collision. The collision can be elastic, inelastic, or ionizing. sibly interrupted by Monte Carlo collision processes. Elec-
In inelastic collisions, electron energy is partially converted trons can experience more than one collision during one time
into molecular excitation; in ionizing collisions, electron en- step ⌬t.
ergy is consumed to split the neutral into a positive ion and a In more detail, position and velocity of the electrons are
second free electron. The probability distribution of the dif- determined in continuous phase space from their Newtonian
ferent collision processes depends on the electron energy at equation of motion according to the electric field at their
the moment of impact; we use the cross section data from the initial position within the time interval. The commonly used
41
SIGLO database. As the collisions are random within a prob- integration is the leap-frog method,48 in which the electron
ability distribution, the actual occurrence of a specific colli- positions and velocities are offset in time by ⌬t / 2.
sion within a sample is determined by a Monte Carlo pro- For the electron-neutral collisions, time, type and scat-
cess. tering angles are determined in a Monte Carlo process by
Once an elastic or inelastic collision process is chosen, sequences of random numbers. For N2, the maximal average
the energy loss of the electron and therefore the absolute collision frequency ␯max is about 9.7⫻ 1012 / s, therefore the
value of its velocity after the collision is fixed. However, minimal average collision time Tmin = 1 / ␯max is about 0.1 ps.
model results will depend on the angular distribution of the By introducing the so-called null collisions, in which no col-
emitted electrons, which again follows a probability distribu- lisions occur, Tmin can be chosen as average collision time
tion. Different scattering methods have been discussed in the independently of the electron energy ␧. The probability P共t兲
literature.42–46 Here we will only focus on the scattering that an electron will travel a time t without collision 共includ-
method used in the present paper. ing null collisions兲 is

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-4 Li et al. J. Appl. Phys. 101, 123305 共2007兲

P共t兲 = e−␯maxt . 共4兲 je = − ␮Ene − D ⵜ ne , 共7兲

Therefore the next collision time ⌬tcollision of an electron is where ne is the electron density, je = une is the flux density,
drawn in a Monte Carlo process from the distribution and u = 具v典 is the mean velocity of electrons. S is the source
of electrons due to collisions and impact ionization, ␮ rep-
1 resents the mobility and D is the diffusion matrix.
⌬tcollision = Tmin ln , 共5兲 The coefficients S, ␮ and D appearing in Eqs. 共6兲 and
p2
共7兲 are to be obtained from elsewhere. One common
where p2 is again a random number. When a collision occurs, approach56 is to solve the Boltzmann equation for a homo-
the energy of the incident electron is calculated, and the dis- geneous and constant electric field E within a background
tribution of the collision processes is determined according gas of constant density. In a uniform electric field, the elec-
to the collision frequencies; then a random number deter- trons gain energy from the field and loose it in inelastic col-
mines the collision type 共null collision, elastic, excitation or lisions, reaching some steady state transport conditions.54,55
ionization collision兲. At the collision, the electron velocities We will derive our coefficients directly from swarm experi-
are changed according to the processes discussed in Sec. ments in the particle model in the next section. Furthermore,
II A 1. Then the next collision time for the particle is deter- the electron source term can be written as
mined. This approach is described in more detail in Refs. 45 S = 兩ne␮共E兲E兩␣共E兲 共8兲
and 49.
At each time step of length ⌬t, the electric field is cal- when attachment and recombination can be neglected. Using
culated on the grid with mesh ⌬ᐉ. First, the number of el- these coefficients in a given gas and density as a function of
ementary charges n p − ne within a grid cell is counted; it di- the electric field is called the local field approximation.
rectly determines the charge density within the cell. Then the Of course, this fluid model has to be extended by conti-
change of electric field components normal to the cell faces nuity equations for other relevant excited or ionized species.
are determined from the densities within the cells through the For a nonattaching gas with neglected ion mobility, the con-
Poisson equation. This simple interpolation on cells of ap- tinuity equation for the density n p of positive ions has to be
propriate size causes no artifacts as we are dealing with par- included,
ticles carrying just one elementary charge e, not with super- ⳵n p
particles. The condition on the cell size is that 共i兲 it is large = S. 共9兲
⳵t
enough that a elementary charge in the cell center does not
create substantial fields on the cell boundary, and 共ii兲 it is Alternatively, the fluid model 共6兲–共9兲 can also be motivated
small enough that no strong density gradients occur between by physical considerations and conservation laws.17,18,57 The
neighboring cells. Here it should be noted that density gra- continuity equations coupled to the Poisson equation for the
dients due to particle number fluctuations are strongly sup- electric field,
pressed when we deal with real particles, not superparticles.
e共n p − ne兲
Therefore more involved interpolation methods such as par- ⵜ·E= . 共10兲
ticle in cell48 共PIC兲 are not required. ⑀0
The choice of the spatial and temporal mesh determines In the present paper, the highest possible consistency
the computational accuracy as well as the computational between particle and fluid model is achieved by determining
costs. We have tested different meshes in planar fronts as the transport coefficients and ionization rate ␮共E兲, matrix
described in Sec. III A. The results, most prominently the D共E兲, and ␣共E兲 for the fluid model from swarm experiments
ionization density behind the front, converge for a sufficient in the particle model; this will be done in Sec. II C.
discretization. However, a balance has to be found between Solutions of the particle model and of the fluid model in
computational accuracy and computational costs. We choose local field approximation will differ when the electric field or
the time step as ⌬t = 0.3 ps, which is of the same order as the the electron density vary rapidly in space or time as the elec-
minimal average collision time Tmin, and the cell size as trons then will not fully “equilibrate” to the local electric
⌬ᐉ = 2.3 ␮m, which is the basic length scale according to field.21 We will investigate these deviations further below.
dimensional analysis in Ref. 17 On this mesh, the electron
density behind planar fronts has an error of less than 0.2%.
2. Numerical implementation
The fluid equations are solved on a uniform grid where
B. The fluid model the electron densities ne and ion densities n p are calculated at
1. Basic equations the centers of the grid cells. The densities can be viewed as
averages over the cell like in the particle model. The field
Fluid models in general can be derived from the Boltz- strength is also calculated at the cell centers. The electric
mann equation.50–53 They approximate the motion of charged field components are taken on the cell faces, where they
particles by continuity equations determine the mass fluxes.
The equations for the particle densities are discretized in
⳵ne space with the finite volume method, based on mass balances
+ ⵜ · je = S, 共6兲
⳵t for all cells. The particle densities are updated in time using

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-5 Li et al. J. Appl. Phys. 101, 123305 共2007兲

the third order upwind-biased advection scheme combined

with a two-stage Runge-Kutta method. For the details and
the tests of the algorithm, we refer to Ref. 18
Analytical studies17,35 and numerical investigations18,29
show that the ionization front is a pulled front; therefore, a
very fine numerical grid is required in the leading edge re-
gion of the front, and standard refinement techniques refining
in the interior front region fail. Like for fronts in the particle
model, we also have tested different numerical meshes for
fronts in the fluid model; these fronts are treated in Sec.
III A. On a too coarse grid, the front moves too fast and is
too smooth due to numerical diffusion of the electron den-
sity. To achieve the same numerical accuracy below 0.2% as
for the particle model, the fluid model requires an approxi-
mately four times finer mesh, namely, ⌬ᐉ = 0.575 ␮m and
⌬t = 0.075 ps. This mesh will be used below. FIG. 2. Average electron energy as a function of time for three different
electron swarms in a field of 100 kV/ cm. Starting with electron swarms
C. Kinetics and transport of electron swarms in moving into the drift direction with an identical kinetic energy of 50 eV
共solid兲, 5 eV 共dotted兲, and 0.5 eV 共dashed兲, all electron swarms approach a
constant fields mean electron energy characteristic for the applied field within 2 – 3 ps.
Swarm experiments deal with electron swarms moving
and multiplying in a constant electric field without changing number of electrons flies in the direction of the field and a
it. If the field is high enough such that the electron number decreasing number against it, and the number of electrons
grows measurably, such a swarm is also called an avalanche. with high kinetic energy increases.
Swarms or avalanches in homogeneous fields are important
experimental and theoretical tools to investigate the electron
dynamics. 2. Gaussian swarm profiles and transport coefficients
Swarm experiments are used as well to determine mo-
1. Particle swarm kinetics: Approach to steady state bilities, reaction rates, and diffusion constants
Particle swarm experiments can be used to determine experimentally.58 An electron swarm drifts, broadens, and
transport coefficients and also to study the particle kinetics. grows under the influence of a constant electric field. The
We will first study the second issue, namely, the relaxation of same experiment is performed here for this purpose, but now
electrons to a steady state velocity distribution and the ve- numerically with the particle model.
locity distribution itself. Indeed the electron swarm will rap- We here recall the essentials: A single electron will ge-
idly equilibrate to the applied field. In such a balanced state, nerically evolve into a swarm that has a Gaussian profile in
the electrons on average gain as much energy from the elec- space. In terms of the fluid model 共6兲–共9兲, this Gaussian dis-
tric field as they loose in inelastic and ionizing collisions; tribution is given by
this is how they reach an energy and velocity distribution
specific for the electric field. The time that the electrons need
to get in balance with the local electric field is an important
indication for the validity of the local field approximation.
We therefore test it here within a particle swarm experiment.
Figure 2 shows how different electron swarms converge
to the same mean energy within a field of 100 kV/ cm. The
experiment starts with a group of electrons of identical ve-
locity directed in the electron drift direction; their kinetic
energy is 50, 5, and 0.5 eV. When the swarms start to drift,
their average energies converge to the same constant value
within at most 3 ps in all three cases. Here the average en-
ergy is used as the simplest indication of their energy and
velocity distribution function.
After reaching steady state motion, swarms demonstrate
a steady state distribution of electron velocities and energies
in a given electric field. Figure 3 shows the distribution of
the longitudinal electron velocity vz in an electric field of 50,
100, 150, and 200 kV/ cm. The figure shows that with in-
FIG. 3. 共Color online兲 Distribution of the electron velocity vz in the longi-
creasing field, the electron velocity distribution deviates
tudinal direction in particle swarm experiments in fields of 50, 100, 150, and
more and more from the Maxwellian profile and therefore 200 kV/ cm. The higher the field, the more the distribution deviates from
from symmetry about velocity vz = 0; rather an increasing Maxwellian and from symmetry about velocity vz = 0.

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-6 Li et al. J. Appl. Phys. 101, 123305 共2007兲

ne共x,y,z,t兲 ⬀ e ␮E␣共E兲t e
−共x2+y 2兲/关4␲DT共t−t0兲兴

4␲DT共t − t0兲
再
function共E兲 = unit · exp a + b ln
E
kV/cm

⫻
e−共z − z0 − ␮Et兲
冑4␲DL共t − t0兲
2/关4␲D 共t−t 兲兴
L 0
, 共11兲 +
c共kV/cm兲
E
± 冋
d共kV/cm兲
E
册冎
2
, 共16兲

where “⫾” is “⫹” for ␣ and ␧, and “⫺” for ␮, DT, and DL.
共The discharge specific context of this solution can be found
Their units and parameters a, b, c, and d are listed in Table I.
in Refs. 57 and 59.兲 Here the center of the package is at
These coefficients will be used in the fluid model 共6兲–共9兲
共x , y , z兲 = 共0 , 0 , z0兲 at time t = 0, and the field is in the z direc-
to reach optimal agreement between particle and fluid model.
tion. The longitudinal and transversal components of the dif-
fusion matrix are denoted as DL and DT.
The transport and reaction coefficients can be deter-
mined from this equation by
III. SIMULATIONS OF PLANAR FRONTS
具z共t2兲典 − 具z共t1兲典
␮共E兲兩E兩 = , 共12兲 A. Concepts and setup of planar ionization fronts
t2 − t1
1. The role of planar fronts in the inner analysis of the
streamer structure
1 ln Ne共t2兲 − ln Ne共t1兲
␣共E兲 = , 共13兲
␮共E兲兩E兩 t2 − t1 Fluid model simulations of streamers within the past
20 years 共see, e.g., Refs. 11, 13, 18, and 20兲 have shown that
具x2共t2兲 + y 2共t2兲典 − 具x2共t1兲 + y 2共t1兲典 the streamer head is surrounded and preceded by an ioniza-
DT共E兲 = , 共14兲 tion front that propagates into the nonionized gas. Within the
4共t2 − t1兲
ionization front, ionization grows until space charge effects
具关z共t2兲 − 具z共t2兲典兴2典 − 具关z共t1兲 − 具z共t1兲典兴2典 set in. The formed space charge layer is much thinner than
DL共E兲 = , 共15兲 the radius of the streamer, it leads to a screening of the elec-
2共t2 − t1兲
tric field in the interior of the streamer head and to a field
where Ne共t兲 is the total number of electrons at time t, and enhancement ahead of it. Therefore the field dependent ion-
具¯典 denotes the average over all particles. ization reaction coefficient ␣共E兲 is enhanced ahead of the
space charge layer and suppressed behind it. The space
charge layer around the streamer is shown in Fig. 1; for a
3. Fluid parameters determined from particle swarms further discussion of the three-dimensional structure and
growth of streamers, we refer to the literature, see, e.g., Ref.
We have determined ␮共E兲, DT共E兲, DL共E兲, and ␣共E兲, and
26.
also the average electron energy ⑀共E兲 in particle swarm ex-
It is clear that the full configuration of the electric field
periments for 42 different background electric fields ranging
can only be analyzed within a two-or three-dimensional set-
from 2 to 205 kV/ cm.
ting. On the other hand, within the inner structure of ioniza-
To obtain the transport coefficients and mean values with
tion front and space charge layer, the electric field does not
satisfactory statistics, one needs a sufficient number of elec-
deviate much from a planar configuration. To analyze the
trons that have experienced an adequate number of colli-
processes within the ionization front in detail, it is therefore
sions. The experiments start from a number of electrons at
advisable to study the inner structure of a planar front. This
the same position 共i.e., located in a single point, which is, a
will be done here. The results can be put in further physical
Gaussian with zero width兲, and end with a swarm of elec-
context through a separate analysis on the inner and the outer
trons with a Gaussian distribution as described in Eq. 共11兲.
scales of the structure as commonly done in hydrodynamic
Because the ionization rate depends strongly on the electric
boundary layer analysis, reaction-diffusion systems, etc.38–40
field strength, the number of initial electrons and the simu-
lation time is chosen according to the fields. For example,
the simulation starts with 106 electrons at 2 kV/ cm and lasts
for 1500 ps, but for 205 kV/ cm, the simulation starts with 2. Construction of planar fronts in the particle
102 electrons and ends with 4 ⫻ 106 electrons after 30 ps. model
As there is some initial transient during which the elec-
trons equilibrate to the field and approach a Gaussian density The construction of a planar front is straightforward in
profile, the transport and reaction coefficients are evaluated the fluid model: gradients ⵜ are simply evaluated in one
according to Eqs. 共12兲–共15兲 at appropriate times t1,2. We spatial direction. We choose this direction to be the z direc-
choose t2 as the end of a swarm experiment, and t1 = t2 / 2 in tion. In the particle model, on the other hand, electrons move
the middle of an experiment. In view of the relaxation times in all three spatial dimensions. Therefore a three-dimensional
below 3 ps evaluated above, this choice of t1 is on the very setting has to be retained. An essentially one-dimensional
safe side. setting is achieved by considering only a small transversal
The numerical results for ␮, DT, DL, ␣ and ⑀ for differ- area A of the front and by imposing periodic boundary con-
ent electric fields E = 兩E兩 are presented in Fig. 4 together with ditions at the lateral boundaries. Furthermore, the electric
empirical fit formulas. These formulas are in the form of field is calculated only in the forward direction z through

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-7 Li et al. J. Appl. Phys. 101, 123305 共2007兲

FIG. 4. Electron mobility, diffusion

rates, ionization rate, and average
electron energy in nitrogen. Plotted are
the reduced coefficients ␮ p, ␣ / p, DT,
DL, and ␧ as a function of reduced
field E / p at room temperature. Our
units are related to other commonly
used units such as 1 kV cm−1 bar−1
= 1.316 V cm−1 torr−1 = 0.424 Td
共1 Td= 10−17 V cm2兲 at T = 300 K.

Ez共z,t兲 = Ez共z0,t兲 + 冕 ⬘冕
z

z0
dz
A
dxdy e共n p − ne兲共x,y,z⬘,t兲
A ⑀0
,
The density fluctuations projected onto the forward di-
rection depend on the transversal area A over which the av-
erages are taken. When A increases, the total number of elec-
共17兲 trons in the simulation increases proportionally to A, while
where Ez is the electric field in the z direction, and z0 can be the relative density fluctuations decrease like 1 / 冑A. There-
any arbitrary position. This means that fluctuations of the fore some intermediate value of the area A has to be chosen:
transversal field due to density fluctuations in the transversal On the one hand, there should be a sufficient number of
direction are not included. In fact, the numerical implemen- electrons to reach a satisfactory statistics, but on the other
tation of Eq. 共17兲 is performed on a grid in the forward hand, there should not be so many electrons that the com-
direction only as discussed in Sec. II A 2. puter runs out of memory within the time interval of interest.
In the simulation, we use a small A for high electric field and
TABLE I. The units and parameters in the empirical fit formula Eq. 共16兲 for a large A for low electric field. For example, we choose the
reduced coefficients and average energy. transversal averaging area as A = 6⌬ᐉ ⫻ 6⌬ᐉ for
−100 kV/ cm, but for −50 kV/ cm, the transversal averaging
Function Unit a b c d
area is A = 20⌬ᐉ ⫻ 20⌬ᐉ, here ⌬ᐉ = 2.3 ␮m.
␮ m2 / 共V s兲 −4.02 0.21 5.44 2.42 In the z direction, the system length is 500⌬ᐉ which
␣ 1/m 12.5 0.16 −200 19.2 allows the front to propagate freely for all runs reported in
DT m2 / s −4.71 0.64 4.80 1.84
this paper. The electric field in the nonionized region at large
DL m2 / s −6.75 1.18 7.89 2.49
␧ eV −1.37 0.78 −4.44 3.46
z is specified by E = E+ẑ, where E+ ⬍ 0 and ẑ is the unit
vector in the z direction.

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-8 Li et al. J. Appl. Phys. 101, 123305 共2007兲

FIG. 5. 共Color online兲 Spatial profile

of a particle front within a field of
−100 kV/ cm at time t1 = 450 ps 共left兲
and t2 = 900 ps 共right兲. 共a兲 Electron
density 共solid兲 and ion density
共dashed兲, 共b兲 electric field 共dashed兲
and net charge 共solid兲, 共c兲 density of
electrons with energy above 0, 15.6,
20, 30, and 50 eV, and 共d兲 zoom into
the front region with high electron en-
ergies, same quantities as in 共c兲.

In the simulations, two different initial conditions are = n−p . 共Upper indices ⫾ indicate quantities before ⫹ and be-
used. In the first, Ne electrons are evenly distributed in a thin hind ⫺ the ionization front.兲 The qualitative features of the
layer of area A with an extension of 19.5⌬ᐉ ⬍ z ⬍ 20.5⌬ᐉ in front are the same as those in the fluid model.17
the field direction. Choosing Electrons with energies above the ionization threshold of
15.6 eV are shown in the lower two panels in Fig. 5; they
Ne/A = 兩E+兩⑀0/e, 共18兲
exist essentially only in the high field region. Electrons with
the field behind the layer is screened according to Eq. 共17兲. energies above 30 eV are so rare that they are hardly seen
Another choice is to begin with a few seed electrons which even on the scale of panel 共d兲. Electrons with energy above
will create an ionization avalanche and form a charge layer 50 eV exist, but cannot be distinguished within this plot. The
later. profiles of high energy electrons also move with the whole
front without change of shape, up to fluctuations.
B. Planar fronts in the particle model Following the track of single electrons of high energy,
1. Qualitative discussion of typical results we found that they gain and loose energy in few collisions
within a few picoseconds and do not run away. In that sense
We first present results in a field of E+ = −100 kV/ cm. they are in a fast dynamic equilibrium with the electrons at
The initial condition is a thin electron layer with total elec- lower energies. This observation agrees with the fast relax-
tron number Ne as in Eq. 共18兲, screening the electric field ation of 50 eV electrons travelling in the forward direction
behind the layer. Figure 5 presents the evolution at times t1 whose energy relaxation is shown in Fig. 2.
= 450 ps 共left兲 and t2 = 900 ps 共right兲. Panel 共a兲 shows the
density distribution of the electrons 共solid兲 and the ions
共dashed兲. Panel 共b兲 shows the net negative charge distribu-
tion 共solid兲 and the electric field 共dashed兲. Panel 共c兲 shows
the total charge density of electrons and the charge density of 2. Quantitative results in different fields
electrons with an energy higher than 0, 15.6, 20, 30, and
50 eV, where 15.6 eV is the ionization energy. Panel 共d兲 The ionization front in a given field E+ is characterized
zooms into panel 共c兲, both in space and in electron densities. by a velocity v, a degree of ionization n−e = n−p behind the
The figure shows an ionization front propagating to the front, and an electron energy distribution in the high field
right. Up to fluctuations, the spatial profiles are essentially region. We now present these quantities in detail.
unchanged; therefore, the front velocity v is essentially con- We define the front position as the position of the maxi-
stant as well. The front carries an overshoot of electrons, mal electron density. Table II summarizes our numerical re-
generating a thin space charge layer that screens the electric sults on the front velocity v as well as on the saturated elec-
field behind the front. In this screened interior region, the tron density n−e behind the front as a function of the electric
electron and ion density reach an equal constant density n−e field E+ immediately ahead of the front.

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-9 Li et al. J. Appl. Phys. 101, 123305 共2007兲

TABLE II. Numerical results on planar fronts in the particle model: front TABLE IV. The degree of ionization ne− behind the front in the fluid model
velocity v and ionization level ne− behind the front as a function of the as a function of field E+: numerical result ne− and analytical upper bound
−
electric field E+ ahead of the front. ne,bound as derived in the Appendix.

E+ v ne− E+ ne− −
ne,bound
共kV/cm兲 共m/s兲 共1 / cm3兲 共kV/cm兲 共1 / cm3兲 共1 / cm3兲

50 共2.773± 0.007兲 ⫻ 105 共5.923± 0.031兲 ⫻ 1011 50 5.43⫻ 1011 5.80⫻ 1011
75 共4.845± 0.023兲 ⫻ 105 共4.372± 0.011兲 ⫻ 1012 75 3.83⫻ 1012 3.98⫻ 1012
100 共7.258± 0.062兲 ⫻ 105 共1.422± 0.003兲 ⫻ 1013 100 1.17⫻ 1013 1.22⫻ 1013
125 共1.012± 0.010兲 ⫻ 106 共3.233± 0.007兲 ⫻ 1013 125 2.49⫻ 1013 2.61⫻ 1013
150 共1.365± 0.008兲 ⫻ 106 共6.014± 0.006兲 ⫻ 1013 150 4.35⫻ 1013 4.58⫻ 1013
175 共1.745± 0.027兲 ⫻ 106 共9.875± 0.020兲 ⫻ 1013 175 6.70⫻ 1013 7.09⫻ 1013
200 共2.262± 0.063兲 ⫻ 106 共1.486± 0.004兲 ⫻ 1014 200 9.50⫻ 1013 10.1⫻ 1013

C. Planar fronts in the fluid model D. Comparison of planar fronts in particle and fluid
model
For planar fronts in the fluid model, there are not only
1. Detailed investigation in a field of −100 kV/ cm
numerical, but also analytical results; both agree within the
numerical accuracy. First, the velocity of the front in a field Now the stage is set to compare planar fronts in the
E+ is given by particle model to those in the fluid model. The comparison is
first done in detail for a planar front propagating into a field
v* = ␮共E+兲兩E+兩 + 2冑DL共E+兲␮共E+兲兩E+兩␣共E+兲, 共19兲 of −100 kV/ cm. Both fluid and particle simulations are car-
ried out in the same setup starting from the same initial con-
according to Refs. 17 and 37 with a slight generalization ditions, i.e., from an electrically neutral group of 100 elec-
along the general lines discussed in Ref. 35. Note that for the trons and ions evenly distributed within the thin layer
initial conditions treated in this paper, the velocity v* is ap- 19.5⌬ᐉ ⬍ z ⬍ 20.5⌬ᐉ and within the transversal area A
proached from below after an algebraically slow relaxation.35 = 6⌬ᐉ ⫻ 6⌬ᐉ.
The electron and ion densities decay exponentially like Figure 6 shows the temporal evolution of the planar
−z/ᐉ*
e in the leading edge of the front where the electric field front. We compare the spatial profile of the electron density
is approximately E+.17,35 The decay length is 共solid line兲 and ion density 共dotted兲 in a particle simulation

冑
with the electron density 共dashed兲 and ion density 共dot
DL共E+兲 dashed兲 in a fluid simulation. Two features are clearly vis-
ᐉ* = . 共20兲
␮共E+兲兩E+兩␣共E+兲 ible: First, the particle and the fluid front move with approxi-
mately the same velocity and the same density profile in the
These analytical results are summarized in Table III. leading edge of the front where the electric field does not
For the degree of ionization behind the front n−e , there is vary yet. Second, the maximal electron density in the front
no closed analytical solution. The ionization behind the front and the saturation level of the ionization behind the front in
can be derived numerically. Furthermore, in the Appendix, the particle model are about 20% higher than in the fluid
we derive an analytical upper bound for this quantity, model. These results of visual inspection agree with those of
namely, Tables I and III for a field of −100 kV/ cm.

n−e
ne,bound
−
=
⑀0
e
冕
0
兩E+兩
␣共e兲de. 共21兲

Numerical result and analytical bound are summarized in

Table IV.

TABLE III. Exact analytical results for planar fronts in the fluid model,
evaluated with the transport coefficients from Figs. 4共a兲–4共d兲: front velocity
v* and electron density decay length ᐉ* as a function of the electric field E+.

E+ v* ᐉ*
共kV/cm兲 共m/s兲 共m兲

50 2.68⫻ 105 7.83⫻ 10−6

75 4.70⫻ 105 3.91⫻ 10−6
100 7.14⫻ 105 2.72⫻ 10−6
125 9.88⫻ 105 2.15⫻ 10−6 FIG. 6. Temporal evolution of the electron and ion densities in a planar front
150 1.29⫻ 106 1.84⫻ 10−6 in a field of E+ = – 100 kV/ cm. Shown are the spatial profiles of electron and
175 1.63⫻ 106 1.66⫻ 10−6 ion densities derived with the particle or the fluid model at time steps t
200 1.98⫻ 106 1.54⫻ 10−6 = 0.09, 0.36, 0.63, and 0.9 ns 共solid lines, ne,part; dashed, ne,fluid, dotted,
n p,part; and dot dashed, n p,fluid兲.

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-10 Li et al. J. Appl. Phys. 101, 123305 共2007兲

FIG. 9. Relative difference of ionization level behind the front 共22兲 共dotted,
FIG. 7. Zoom into the particle ionization front of Fig. 6 in a field of above兲, front velocity 共23兲 共dashed, below兲, and mean electron energy in the
−100 kV/ cm at time 0.63 ns: shown are electron density distribution ne leading edge 共24兲 共solid, middle兲 between particle and fluid model as a
共solid line兲, local average electron energy ␧part 共dotted line兲, local average function of the electric field E+ ahead of the front.
electron energy ␧LFA according to the local field approximation 共dashed
line兲, and electric field strength E 共dot-dashed line兲.
for the particular spatial region of the front where the electric
As we have excluded other reasons of the discrepancy field has decreased by 10% to a value of −90 kV/ cm. More
such as numerical discretization errors or inconsistent trans- precisely, electrons are collected from the first cell where
port and reaction coefficients, deviations must be due to the 兩E兩 ⬎ 90 kV/ cm, searching with increasing z. The average
approximations in the fluid model, and a closer inspection field in these cells is about E = −91.47 kV/ cm. To reach a
shows that we should focus on the electron energies. satisfactory statistics, the electrons are collected from ten
Figure 7 zooms into the ionization front shown in Fig. 6 different instants of time with a temporal distance of 30 ps
at time t = 0.63 ns. Here we show the electron density 共solid between the two consecutive sampling times to ensure statis-
line兲 and electric field 共dot-dashed line兲 in the particle model. tical independence. While these data are plotted as a solid
Furthermore the local mean energy of the electrons in the line, for comparison the electron energy distributions in the
particle model is indicated with a dotted line. Finally, the swarm experiments in a constant field of −90 and
mean electron energy according to the local field approxima- −100 kV/ cm are plotted as a dashed or dotted line, respec-
tion ⑀共E兲 is derived from the local field E and Eq. 共16兲 for tively. Analyzing the electron energy distribution at low en-
⑀共E兲; it is indicated with a dashed line. It can be seen that the ergies ␧ ⬍ 10 eV, the ionization front at a local field around
average electron energy nicely follows the local field ap- −90 kV/ cm and the swarm experiment in a constant field of
proximation in the interior of the ionized region while it is −90 kV/ cm are quite similar while the distribution for a
considerably higher in the region where the electric field is swarm in a field of −100 kV/ cm is clearly lower. On the
large and the electron density decreases rapidly. other hand, for electron energies above 20 eV, the energy
This deviation is analyzed in more detail in Fig. 8 where distribution in the ionization front at a local field of
the full electron energy distribution is shown. This is done −90 kV/ cm actually lies closer to the distribution of the
swarm at −100 kV/ cm than to that at −90 kV/ cm. This ob-
servation not only confirms that the average electron energy
in the leading edge of the ionization front is higher than in
the local field approximation, as is consistent with Fig. 7, but
it indicates that the higher mean energies correspond to a
higher population of electron energy states above the ioniza-
tion threshold. We therefore expect that also the ionization
rates are higher in the particle model than in the local field
approximation.

2. Results for other fields

FIG. 8. 共Color online兲 The electron energy distribution function is measured Having analyzed the front propagating into a field of
as P共␧兲 / 冑␧, where P共␧兲 is the probability of electron energy ␧. Here we E+ = −100 kV/ cm in detail, we now summarize equivalent
show the electron energy distribution at the front region of Fig. 7 where E
⬇ −90 kV/ cm 共solid line兲 and in the swarm experiments in constant fields of results for fields ranging from −50 to − 200 kV/ cm in Fig. 9.
−90 kV/ cm 共dashed line兲 and −100 kV/ cm 共dotted line兲. The figure shows the relative difference

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-11 Li et al. J. Appl. Phys. 101, 123305 共2007兲

− − is an electron energy overshoot in the leading edge of the

ne,part − ne,fluid
− 共22兲 ionization front. So the most prominent deviation from the
ne,fluid
fluid model visible in Fig. 7 occurs in the region were the
of the saturated electron density behind the front between the electric field is almost constant; furthermore, if it would be a
particle and the fluid model, the relative difference consequence of the local field approximation, the deviation
would have the opposite sign.
vpart − vfluid
共23兲
vfluid
of the front velocity and the relative difference 2. Density gradient and relation between front and
swarm experiments
具␧part典 − 具␧LFA典
共24兲
具␧LFA典 We conclude that the electron energy overshoot in the
leading edge of the ionization front has to be due to the
of the mean electron energy between the particle model and electron density gradient rather than to the field gradient. The
the local field approximation at a point in the ionization front effect of density gradients can be tested in swarm experi-
where electron densities are low and the electric field is just ments 共cf. Sec. II C兲 in a constant electric field as well. It
slightly screened to the value E = 0.98E+. All differences in- turns out that such a test goes beyond qualitative results and
crease with increasing field. actually allows for a quantitative comparison of fronts and
swarms as will be explained below.
E. Interpretation of results The key to the quantitative comparison is based on two
1. Electric field gradient, electron motion in the front, closely related facts: 共i兲 The density profiles in the leading
and electron energy overshoot edge of a swarm or avalanche in a constant field and in the
leading edge of an ionization front penetrating the same field
A few authors have dealt with particle discharge models
are both given by
with strong gradients both in the electric field and in the
*
particle density. They analyzed deviations between particle ne ⬃ e−z/ᐉ , 共25兲
models and the fluid model in local approximation by means
of the Boltzmann equation. They analyzed the case of sta- where ᐉ is given by Eq. 共20兲 and in Table II. 共ii兲 Also the
*

tionary gradients of the electric field54,55 and the case of posi- velocities of swarm and front have the same value v* from
tive streamer ionization fronts.21 They suggested corrections Eq. 共19兲 and Table II.
to the fluid model both due to field gradients and to density This relation between velocities and decay rates of
gradients. swarms and fronts holds generally for any so-called pulled
The situation in negative streamer ionization fronts is front whose velocity is determined in the linearly unstable
somewhat different: a negative ionization front in nitrogen region ahead of the front, as is discussed in a general setting
moves approximately with the electron drift velocity deter- in Ref. 35, specifically in Sec. 2.5.1. The statement can be
mined by the electric field in the leading edge of the ioniza- verified on the explicit form of the Gaussian profile 共11兲. In
tion front. This means that the electrons in this leading edge the coordinate
region on average move within a stationary electric field. ␰ = z − z 0 − v *t 共26兲
More precisely, the front velocity v*共E+兲 from Eq. 共19兲 is *
slightly larger than the electron drift velocity ␮共E+兲兩E+兩 in moving with velocity v , the Gaussian swarm distribution
the electric field ahead of the front. This is because the front 共11兲 can be written as
is not carried by electron drift only, but also by diffusion and −共z − z0 − ␮Et兲2/共4DLt兲 −␰2/共4DLt兲
␮兩E兩␣t e −␰/ᐉ* e
ne ⬀ e =e 共27兲
creation of new free electrons. In a coordinate system mov- 冑4␲DLt 冑4␲DLt .
ing with the ionization front, the existing electrons therefore
on average move backwards. The faster this motion is, the We now test the above theoretical predictions on the
further they are behind. particle model for a swarm and a planar front in a field of
Figure 8 shows that the high energy tail of the electron −100 kV/ cm. The results are shown in Fig. 10. It should be
energy distribution in the ionization front at the position remarked that our electron swarm has a Gaussian density
where the field is −91.47 kV/ cm is closer to the swarm ex- distribution both in the longitudinal and in the transversal
periment at −100 kV/ cm than to that at −90 kV/ cm. One direction. To focus on the profile in the longitudinal direc-
could interpret these results by assuming that the electrons tion, the density in the swarm in Fig. 10共a兲 is taken as the
have gained their energy distribution in a field of close to density on the longitudinal axis. Figure 10共a兲 shows the elec-
−100 kV/ cm and then are transported backwards to where tron density profile of the swarm decaying at both edges and
the field is only −90 kV/ cm. the profile of the planar front that grows and saturates to a
However, Fig. 7 shows that this interpretation cannot be constant level. As the densities are plotted on a logarithmic
true. If the electron energy distribution would first equilibrate scale, an exponential decay like in Eq. 共25兲 amounts to a
to a field of −100 kV/ cm and then partially be transported straight line with slope −1 / ᐉ* in the plot. Despite density
backwards, then the mean electron energy 具⑀part典 would ev- fluctuations and slow transients in the buildup of the
erywhere be below the mean electron energy 具⑀LFA典 at profile,35 Fig. 10共a兲 indeed shows that swarm and front have
−100 kV/ cm in local field approximation. But clearly there a similar decay profile on the right hand side.

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-12 Li et al. J. Appl. Phys. 101, 123305 共2007兲

Now the electron energies in the leading edge substan-

tially exceed the local field approximation and indicate the
presence of a larger fraction of electrons with energies above
the ionization threshold. This leads also to higher ionization
rates than estimated by the local field approximation ␣共E兲.
As the ionization level behind the front n−e = n−p is well ap-
proximated by the effective ionization rate ␣ integrated over
the fields within the front 共21兲 共recall also the argument in
the Appendix兲, it is clear that the ionization level behind the
front is higher in the particle model than in the fluid model.

IV. CONCLUSION

Negative streamer ionization fronts in pure nitrogen

FIG. 10. Comparison of an electron swarm and a planar ionization front, were investigated in planar approximation, both following
both in the particle model and in a field of −100 kV/ cm. 共a兲 Electron density the kinetics of single electrons in a particle model and in a
on a logarithmic scale for swarm 共solid兲 and front 共dotted兲 as a function of
z. 共b兲 Mean electron energy in swarm 共solid兲 and front 共dotted兲 and in the fluid model in local field approximation. As parameter func-
front assuming the local field approximation 共dashed兲 as a function of z. tions for the fluid model were derived from swarm experi-
ments in the particle model and numerical errors are under
This sets the stage to compare now the electron energies strong control, a discrepancy between results of particle and
of swarm and front in Fig. 10共b兲. The dotted and the dashed fluid models must be attributed to the approximations made
lines reproduce lines from Fig. 7 for the planar particle front; in the fluid model.
they show the actual mean electron energy 具␧part典 共dotted兲 For electric fields immediately ahead of the front of
50 kV/ cm or lower—the statements hold for nitrogen under
and the mean electron energy according to the local field
normal condition, they are easily extended to other densities
approximation 具␧LFA典 共dashed兲. For the electron swarm, the
by introducing the reduced electric field E / N—particle and
mean electron energy along the swarm axis is indicated as a
fluid model essentially agree in front speed, profile, and ion-
solid line in Fig. 10共b兲. It is remarkable that this energy
ization level behind the front. When the field increases, the
within the swarm is not constant in space though the electric
velocity does not vary much between both models, but the
field has a constant value. Rather the energy increases almost
ionization level behind the front is substantially larger, as is
linearly along the axis. While the average energy of the
demonstrated for a field of 100 kV/ cm in Fig. 6. In Fig. 9 the
swarm is around 8.5 eV, producing the respective value for
differences between particle and fluid models are summa-
the local field approximation ␧共E兲 in Fig. 4共e兲, the average
rized for the full field region investigated.
energy at the leading edge of the swarm reaches values
Can the discrepancy between the models be attributed to
above 10 eV.
a particular physical mechanism and to a particular part of
Now the leading edge regions of swarm and front need a
the front region? In Fig. 7 we find the largest discrepancy
closer inspection. The leading edge is defined as the region
between fluid and particle models in the leading edge of the
where the electric field is 共almost兲 constant and where the
front where the electron density is very low and where the
electron densities both approach the profile 共25兲. It is in this
electric field is hardly screened. Here the electrons have an
region where the electrons attain the highest mean energies,
average energy considerably higher than estimated by the
and Fig. 10共b兲 shows that the mean energy profiles of swarm
local electric field. This effect can only be explained by an
and front in this region are virtually identical up to fluctua-
effect of the electron density gradient, not of the electric field
tions. We conclude that the high electron energies within the
gradient. Indeed, an electron swarm or avalanche in a con-
leading edge of the ionization front that are shown in Figs. 7
stant electric field has the same velocity and the same density
and 8, are due to the electron density gradient, and can be
gradients in its leading edge and shows the same electron
studied in the leading edge of a swarm experiment as well.
energy overshoot in its front part, as is illustrated in Fig. 10.
In essence, in high fields the gradient length of front and
avalanche ᐉ* becomes of the order of the electron energy
3. The ionization density behind the front
relaxation length; therefore, the fast and energetic electrons
The discussion above shows that within the electron den- can get ahead of the slower ones.
sity gradient, the fast electrons are on average a bit ahead of These higher electron energies in the leading edge of the
the slower electrons as they have a higher average velocity in front lead to higher ionization rates than the ionization rate
the forward direction. This is understandable since on the ␣共E兲 in local field approximation—indeed, ␣共E兲 is the total
one hand, each single fast electron looses its energy over a effective ionization rate of a complete electron swarm 共as
length of ⬃1.5 ␮m 共which is the electron drift velocity derived in Sec. II C兲 while the electrons in the leading edge
␮共E兲E times the relaxation time ⬃2 ps from Fig. 2兲, but on of the swarm have higher ionization rates.
the other hand, the gradient length ᐉ* ⬇ 2.7 ␮m of the den- Now the ionization level behind the front can be ap-
sity profile is of the same order. Clearly, this effect becomes proximated by integrating the effective ionization rate ␣共E兲
more pronounced, the higher the field. over the electric field strength E, independently of the pre-

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-13 Li et al. J. Appl. Phys. 101, 123305 共2007兲

cise spatial structure of the front. This result is derived in the S j = 兩je兩␣共E兲, 共A2兲
Appendix. This implies that higher electron energies and
higher effective ionization rates in the leading edge of the while keeping the equations unchanged otherwise. This
front directly lead to higher ionization levels behind the change of the source term means that more ionization is cre-
front, even though only very few electrons are involved in ated at each part of the front, and therefore the ionization
this dynamics in the low density region. level behind the front will be higher. This is because in a
The effect of few electrons of high energy is much more negative front, drift current and diffusion current point into
severe for the ionization rates ␣ than for the average drift the same direction, and therefore
velocity ␮E; therefore, the most pronounced effect is seen in = 共E兲 · ⵜne − ␮共E兲Ene兩  ␮共E兲兩E兩ne .
兩je兩 = 兩− D 共A3兲
the ionization levels behind the front and less in the front
velocities. The dominant part of the current is the drift or Ohmic part,
We finally remark that next to explicit predictions, our therefore the results with the changed source term 共A2兲
work has delivered two useful insights. 共i兲 The physical dis- should be close to the original problem.
crepancies between particle and fluid model lie in the leading According to Maxwell’s equations, the divergence of the
edge of the front, though the effect is not so much seen in the total current vanishes,
velocity, but much more in the ionization level behind the ⵜ · 共⑀0⳵tE + j兲 = 0, 共A4兲
front. This gives a clue for a numerical strategy combining
efficient features of fluid and particle models. 共ii兲 The essen- where j is the electric current; in our case it is j = −eje with je
tial features in the leading edge of the front are equally from Eq. 共7兲. If the front is planar and if the electric field E+
present in the leading edge of an electron swarm or ava- in the nonionized region does not change in time, then it
lanche in a constant field, where it can be studied much follows immediately that
easier.
⑀0⳵tE = eje . 共A5兲
Now the growth of the ion density 共6兲 is bounded by the
ACKNOWLEDGMENTS source term 共A2兲 and the identity 共A4兲 is inserted,
⑀0
The authors thank W. Hundsdorfer for helpful discus- ⳵tn p = S = ␮共E兲兩E兩ne␣共E兲
S j = je␣共E兲 = ⳵tE␣共E兲.
sions about the numerical solution of the fluid model. They e
acknowledge support by the Dutch national program BSIK, 共A6兲
in the ICT project BRICKS, theme MSV1.
The first and last expressions in this inequality can be inte-
grated in time with the result

APPENDIX: APPROXIMATING THE IONIZATION

LEVEL BEHIND THE FRONT
n p共x,t兲 − n p共x,0兲

⑀0
e
冕 E共x,t兲

E共x,0兲
␣共e兲de. 共A7兲
In previous papers,17 we have established that for a pla-
nar, uniformly propagating front with field independent elec- Now the time interval 关0 , t兴 is taken in such a way that it
tron mobility ␮共E兲 = const and for vanishing electron diffu- contains the time range in which the ionization front passes
sion D共E兲 = 0, there is a unique relation between the field over the point x of observation. n p共x , 0兲 is then the ion den-
ahead of the front E+, the ionization rate function ␣共E兲, and sity before and n p共x , t兲 the ion density behind the front. As a
the degree of ionization n−e = n−p behind the front, namely, result, we find that the expression 共A1兲 indeed is an upper
bound and good approximation for the ionization level be-
n−e = n−p =
⑀0
e
冕
0
兩E+兩
␣共e兲de 共A1兲
hind the front.
We finally note if there is electron attachment and posi-
tive and negative ions n p,n are formed, the statement stays
in dimensional units. true for the total ion charge density n = n p − nn if the total
Such a simple identity does not hold for the full fluid source and sink term for the ion density n can be written in
model 共6兲–共10兲 with electron diffusion, but we can establish the form 共A2兲.
the expression 共A1兲 as an upper bound for the free electron 1
density n−e behind the front; this upper bound is actually a L. B. Loeb and J. M. Meek, The Mechanism of the Electric Spark 共Clar-
endon, Oxford 1941兲.
very good approximation to the real value as Table IV shows. 2
H. Raether, Electron Avalanches and Breakdown in Gases 共Butterworths,
The upper bound is constructed as follows: London, 1964兲.
In the ionization source term S 共8兲, we follow the usual
3
B. Lay, R. S. Moss, S. Rauf, and M. J. Kushner, Plasma Sources Sci.
Technol. 12, 8 共2003兲.
procedure to take only the drift term of the electron current 4
A. Bhoj and M. J. Kushner, J. Phys. D 37, 2510 共2004兲.
into account. We note at this place that this assumption re- 5
Electrical Discharges for Environmental Purposes: Fundamentals and Ap-
quires some reconsideration, and one could argue that the plications, edited by E. M. van Veldhuizen 共Nova Science, New York,
complete current should be taken into account and that a 2000兲.
6
L. R. Grabowski, E. M. van Veldhuizen, A. J. M. Pemen, and W. R.
decomposition into a drift and a diffusion part is artificial. Rutgers, Plasma Chem. Plasma Process. 26, 3 共2005兲.
This argument will have to be elaborated at some other 7
E. M. Bazelyan and Yu. P. Raizer, Lightning Physics and Lightning Pro-
place. We here just change the source term into tection 共Institute of Physics, Bristol, 2000兲.

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp
123305-14 Li et al. J. Appl. Phys. 101, 123305 共2007兲

8
J. R. Dwyer, Geophys. Res. Lett. 30, 2055 共2003兲. 35
U. Ebert and W. van Saarloos Physica D 146, 1 共2000兲.
9
E. R. Williams, Plasma Sources Sci. Technol. 15, S91 共2006兲. 36
O. Chanrion and T. Neubert, “2d axisymmetrical particle modelling of the
10
D. D. Sentman, E. M. Wescott, D. L. Osborne, D. L. Hampton, and M. J. production of thermal runaways electron by sprite streamers,” Geophysi-
Heavner, Geophys. Res. Lett. 22, 1205 共1995兲. cal Research Abstracts, Vol. 9, p. 08389 共2007兲.
11
S. K. Dhali and P. F. Williams, J. Appl. Phys. 62, 4696 共1987兲. 37
A. N. Lagarkov and I. M. Rutkevich, Ionization Waves in Electrical
12
J. M. Guo and C. H. J. Wu, IEEE Trans. Plasma Sci. 21, 684 共1993兲. Breakdown of Gases 共Springer, Berlin, 1994兲.
13 38
P. A. Vitello, B. M. Penetrante, and J. N. Bardsley, Phys. Rev. E 49, 5574 C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for
共1994兲. Scientists and Engineers 共McGraw-Hill, New York, 1978兲.
14
A. A. Kulikovsky, J. Phys. D 27, 2556 共1994兲. 39
P. C. Fife, Dynamics of Internal Layers and Diffusive Interfaces 共SIAM,
15
A. A. Kulikovsky, J. Phys. D 28, 2483 共1995兲. Philadelphia, 1978兲.
16
N. Yu. Babaeva and G. V. Na
dis, J. Phys. D 29, 2423 共1996兲. 40
U. Ebert and W. van Saarloos, Phys. Rep. 337, 139 共2000兲.
17
U. Ebert, W. van Saarloos, and C. Caroli, Phys. Rev. E 55, 1530 共1997兲. 41
W. L. Morgan, J. P. Boeuf, and L. C. Pitchford, SIGLO cross sections
18
C. Montijn, W. Hundsdorfer, and U. Ebert, J. Comput. Phys. 219, 801 共http://www.siglo-kinema. com.兲.
共2006兲. 42
E. E. Kunhardt and Y. Tzeng, Phys. Rev. A 34, 2158 共1986兲.
19 43
P. Ségur, A. Bourdon, E. Marode, D. Beeieres, and J. H. Paillol, Plasma A. Okhrimovskyy, A. Bogaerts, and R. Gijbels, Phys. Rev. E 65, 037402
Sources Sci. Technol. 15 648 共2006兲. 共2002兲.
20 44
A. Luque, U. Ebert, C. Montijn, and W. Hundsdorfer, Appl. Phys. Lett. 90, A. V. Phelps and L. C. Pitchford, Phys. Rev. A 31, 2932 共1985兲.
081501 共2007兲. 45
J. P. Boeuf and E. Marode, J. Phys. D 15, 2169 共1982兲.
21
G. V. Na
dis, Tech. Phys. Lett. 23, 493 共1997兲. 46
M. Surendra, D. B. Graves, and G. M. Jellum, Phys. Rev. A 41, 1112
22
E. E. Kunhardt, J. Wu, and B. Penetrante, Phys. Rev. A 37, 1654 共1986兲. 共1990兲.
23
N. Liu and V. P. Pasko, J. Geophys. Res. 109, A04301 共2004兲. 47
C. B. Opal, W. K. Peterson, and E. C. Beaty, J. Chem. Phys. 55, 4100
24
E. M. van Veldhuizen and W. R. Rutgers, J. Phys. D 35, 2169 共2002兲. 共1971兲.
25 48
T. M. P. Briels, E. M. van Veldhuizen, and U. Ebert, IEEE Trans. Plasma C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simula-
Sci. 33, 264 共2005兲. tion 共Adam Hilger, Bristol, 1991兲.
26 49
U. Ebert, C. Montijn, T. M. P. Briels, W. Hundsdorfer, B. Meulenbroek, A. W. J. M. Brok, Ph.D. thesis, Eindhoven University of Technology, The
Rocco, and E. M. van Veldhuizen, Plasma Sources Sci. Technol. 15, S118 Netherlands 2005.
共2006兲. 50
Molecular Theory of Gases and Liquids, edited by J. O. Hirschfelder, C. F.
27
M. Arrayás, U. Ebert, and W. Hundsdorfer, Phys. Rev. Lett. 88, 174502 Curtiss, and R. B. Bird 共Wiley, New York, 1964兲.
共2002兲. 51
S. Chapman and T. G. Cowling, The Mathematical Theory of Non-
28
A. Rocco, U. Ebert, and W. Hundsdorfer, Phys. Rev. E 66, 035102共R兲 Uniform Gases 共Cambridge University Press, Cambridge, 1974兲.
共2002兲. 52
I. P. Shkarofsky, T. W. Johnston, and M. P. Bachynski, The Particle Ki-
29
C. Montijn, U. Ebert, and W. Hundsdorfer, Phys. Rev. E 73, 065401 netics of Plasmas 共Addison-Wesley, Reading, 1966兲.
共2006兲. 53
E. Gogolides and H. H. Sawin, J. Appl. Phys. 72, 3971 共1992兲.
30 54
T. M. P. Briels, J. Kos, E. M. van Veldhuizen, and U. Ebert, J. Phys. D 39, N. Sato and H. Tagashira, J. Phys. D 18, 2451 共1985兲.
5201 共2006兲. 55
N. L. Aleksandrov and I. V. Kochetov, J. Phys. D 29, 1476 共1996兲.
31
J. R. Dwyer, Geophys. Res. Lett. 31, L12102 共2004兲. 56
G. J. M. Hagelaar and L. C. Pitchford, Plasma Sources Sci. Technol. 14,
32
D. M. Smith, L. I. Lopez, R. P. Lin, and C. P. Barrington-Leigh Science 722 共2005兲.
307, 1085 共2005兲. 57
C. Montijn and U. Ebert, J. Phys. D 39, 2979 共2006兲.
33
J. R. Dwyer, Geophys. Res. Lett. 32, L20808 共2005兲. 58
L. G. H. Huxley and R. W. Crompton, The Diffusion and Drift of Electrons
34
G. D. Moss, V. P. Pasko, N. Liu, and G. Veronis, J. Geophys. Res. 111, in Gases 共Wiley, New York, 1974兲.
A02307 共2006兲. 59
Yu. P. Raizer, Gas Discharge Physics 共Springer, Berlin, 1991兲.

Downloaded 26 Jun 2007 to 131.155.108.56. Redistribution subject to AIP license or copyright, see http://jap.aip.org/jap/copyright.jsp