Coarse-grained analysis of a lattice Boltzmann model for planar streamer fronts

0.49 MB·English

by  Wim Vanroose

Coarse-grained analysis of a lattice Boltzmann model for planar streamer fronts
Wim Vanroose,∗ Giovanni Samaey, and Pieter Van Leemput
Department of Computer Science, Katholieke Universiteit Leuven, Celestijnenlaan 200A, B-3001 Heverlee, Belgium
WestudythetravelingwavesolutionsofalatticeBoltzmannmodelfortheplanarstreamerfronts
that appear in the transport of electrons through a gas in a strong electrical field. To mimic the
physical properties of the impact ionization reaction, we introduce a reaction matrix containing
7 reactionratesthatdependontheelectronvelocities. ViaaChapman–Enskogexpansion,oneisable
0 tofindonlyaroughapproximationforamacroscopicevolutionlawthatdescribesthetravelingwave
0 solution. We propose to compute these solutions with the help of a coarse-grained time-stepper,
2 whichisaneffectiveevolutionlawforthemacroscopicfieldsthatonlyusesappropriatelyinitialized
n simulationsofthelatticeBoltzmannmodelovershorttimeintervals. Thetravelingwavesolutionis
a found as a fixed point of the sequential application of the coarse-grained time-stepper and a shift-
J back operator. The fixed point is then computed with a Newton-Krylov Solver. We compare the
2 resulting solutions with those of the approximate PDE model, and propose a method to find the
minimal physical wave speed.

]
h
p I. INTRODUCTION willbeabletoaccuratelypredictthemicroscopicphysics
- of electron impact on molecular targets such as N and
p 2
m When a gas of neutral atoms or molecules is exposed O2, the most important molecules in the composition of
air. This progress in the understanding of the impact
o toastrongelectricalfield,asmallinitialseedofelectrons
ionization reaction, however, has not been incorporated
c canleadtoanionizationavalanche. Indeed,theseedelec-
. tronsareacceleratedbythefieldandgainenoughenergy inthedescriptionofthemacroscopicbehaviorsuchasthe
s
minimalstreamerfrontofEbertetal. Instead,suchmod-
c to ionize the neutral atoms when they collide. The two
si slow electrons that emerge from this reaction, i.e. the els still make use of a phenomenological approximation
y impactandtheionizedelectron,areagainacceleratedby to the reactions, such as the Townsend approximation.
h the fieldandcause,ontheir turn, anionizationreaction. This article extends the minimal streamermodel andin-
p corporatesmore microscopic information. We model the
Simultaneously the electrical field is locally modified be-
[ system by a Boltzmann equation, which is constructed
cause of the charge creation. This interplay between the
such that the cross sections in the collision integral re-
1 dynamicsoftheelectronsandtheelectricalfieldcanlead
v to a multitude of phenomena studied in plasma physics semble the true microscopic cross sections.
1 such as arcs, glows, sparks and streamers. Tofindthetravelingwavesolutionsofthismoremicro-
3 scopic model, we exploit a separation of time scales be-
In this article we will focus on the initial field driven
0 tweentherelaxationoftheelectrondistributionfunction
ionization that can lead to traveling waves known as
1
to a local equilibrium and the evolution of the macro-
0 streamerfronts. These waveshavepreviously beenstud-
scopic fields (electron density and electrical field). It is
7 ied by Ebert et al. [1] who introduced and analyzed the
0 minimal streamer model, a one-dimensional model for known from kinetic theory that the first process is fast:
/ the propagation of planar streamer fronts. This model once initialized, it takes a molecular gas not more than
s
a few collisions to relax to its equilibrium state.
c consists of two coupled non-linear PDEs: a reaction-
i In kinetic theory, the fast times scales are often elimi-
s convection-diffusion equation for the evolution of the
y electron density and a Poisson-like evolution equation nated from the problem by assuming a local equilibrium
h for the electrical field. The reaction term is based on distribution function which leads to a reaction-diffusion
p model with transport coefficients that depend on the lo-
the Townsend approximation that expresses the growth
v: ofthe number ofelectronsas afunction ofthe localelec- cal electrical field. This reduction method, however, is
i trical field. only successful in the absence of steep gradients in the
X
electron density [18], an assumption not valid for the
Duringthelasttwodecades,however,alotofprogress
r planar streamer fronts that have, typically, very steep
a has been made in the microscopic understanding of im-
increases in the electron density.
pact ionization reactions in atomic and molecular sys-
In this article, we take an alternative route and find
tems. Inthis reactionanimpactelectronionizesthe tar-
the traveling wave solution through a so-called coarse-
get and kicks out an additional electron. There are sev-
grainedtime-stepper (CGTS)thatexploitstheseparation
eral successful theories that can predict the exact prob-
oftimescalestoextracttheeffectivemacroscopicbehav-
ability distribution of the escaping electron [7, 8, 9, 10].
ior. This method was proposed by Kevrekidis et al. [5]
In the next decade, we expect that the theoretical tools
and the numericalaspects of its applicationto find trav-
eling wave solutions of lattice Boltzmann models have
recently been studied [3]. The time-stepper uses a se-
∗Present address: Departement Wiskunde-Informatica, Univer- quence ofcomputationalsteps to evolvethe macroscopic
siteitAntwerpen,Middelheimlaan1,2020Antwerpen,Belgium state. This sequence involves: (1) a lifting step, which
2
creates an appropriate electron distribution function for
a given electron density, (2) a simulation step, where
the lattice Boltzmann model is evolved over a coarse- )
s
grained time step ∆T, and (3) a restriction step, where nit
u
the macroscopic state is extracted from the electron dis-
b
r
tribution function. This method does not derive effec- a
( R
tive equations explicitly, and therefore allows steep gra- n
o
dients to be present. We compare our results with those cti
e
s
obtained by deriving an approximate macroscopic PDE s
s
modelthroughthemoretraditionalChapman-Enskogex- ro
c
pansion. Thepaperthereforeillustratestheapplicability
of the coarse-grainedtime-stepper on a non-trivialprob-
0
lem where the exact macroscopic equations are hard to
0 20 40 60
derive.
Theoutlineofthepaperisasfollows. InsectionII,we Energy (eV)
shortly review the physics of the impact ionization reac-
tion andthe streamerfronts, andrecapitulate the Boltz-
mannequations. InsectionIII,wederivefromtheBoltz- )
s
t
mann equation a lattice Boltzmann model with multiple ni
u
velocities and discuss how the ionization reaction, exter- b
r
nal field and electron diffusion are incorporated in the a
(
model. Section IV derives a macroscopic PDE from the on
ti
model using the Chapman-Enskog expansion and dis- c
e
s
cusses the minimal velocity of the traveling waves. Sec- s
s
o
tion V formulates the coarse-grained time-stepper and r
c
VI how the traveling wave solutions are found. Finally
in section VII, we have some numerical results.
0
0 5 10 15 20 25
Energy (eV)
II. MODEL
Figure 1: A sketch of the typical shape of the impact ion-
A. The physics of the impact ionization reaction
ization cross section ( see the experimental results in [6]).
On top, we show the total cross section as a function of the
The impact ionization reaction is a microscopic reac- impact energy where below a threshold energy of 25eV no
tion where electrons with, typically, an energy around reactions take place. In the proposed model we distinguish
50eV collide with an atom or a molecule and ionize this between slow particles with an energy below this threshold
that do not react and fast particles with their energy above
target. The reaction rates of this process depend sensi-
this threshold. The fast reacting particles experience a cross
tively on a number of parameters. Let us consider the
section of R, as indicated by the dashed line. In the bottom
simplest system: an electron hits a hydrogen atom with
figure, we show the energy differential cross section for the
aboundelectroninitsgroundstate. Whentheincoming
escapingelectron,wherethetotalenergyoftheescapingelec-
electron has an energy larger than the binding energy of
trons is 25eV. Since the two electrons are indistinguishable
the electron in the atom, it can kick out, with a certain there is a symmetry. In our five speed model, we make the
probability,theboundelectronthatwillescape,together approximation that the two electrons can only escape with
with the impact electron, from the atom. The total en- equalenergy sharing.
ergy of the two electrons after the collision is equal to
the energy of the incoming electron minus the original
binding energy of the bound electron. When this cross section is integrated over all angles
Thereactionratesofthisprocessareexpressedbycross (θ ,φ )and(θ ,φ )oftheescapingelectronsandallpos-
1 1 2 2
sections; these are probabilities that a certain event will sible ratios of E /E of the electron energies, we get the
1 2
take place. One such cross section is the triple differ- total cross section. This is the total probability that the
ential cross section, which is the probability to find af- incoming electron will cause an ionization event. This
ter the collision one electron escaping in the direction total cross section depends on the energy of the incom-
(θ ,φ ) with an energy E and a second electron escap- ing electronandis zerowhen the energyof the incoming
1 1 1
ing in the direction (θ ,φ ) with an energy E , where electron is below the binding energy of the bound elec-
2 2 2
the angles are measured with respect to the axis defined tron. Just above this binding energy, there is a steep
by the momentum of the incoming electron. Since the rise in the cross section that is known as a threshold.
twoelectronsrepeleachother,it ismore likelythatthey Just above this threshold the cross section is the largest
escape in opposite directions [13]. and as we further increase the energy the cross section
3
diminishes. This is illustrated in figure 1 (top). reaction-diffusion equation for the evolution of the elec-
Whenthecrosssectionisintegratedoftheanglesonly, tron density with a Boltzmann equation for the one-
butnotovertherelativeenergies,wegetthesocalleden- electron distribution function f(x,v,t), that counts the
ergy differential cross section, whichis the probabilityof number of electrons in the phase-space volume element
causing an ionization event with a given relative energy bounded by position x and x+dx and by speed v and
of the two electrons. In contrast with the electron direc- v+dv. The Boltzmann equation is
tions, there is no pronounced preference for the energy
sharingbetweenthe twoelectrons,seefigure1(bottom). ∂f(x,v,t) ∂f(x,v,t) ∂f(x,v,t)
+v· +E(x,t)· =Ω(x,t),
Itisonlyslightlymorelikelythattwoelectronswillcome ∂t ∂x ∂v
(1)
out with unequal energy.
where E(x,t) is the external electrical field and Ω(x,t)
Recently,severaltheoreticalmethodshavesuccessfully
is the collision operator, an integral operator that inte-
predictedthedirectionsofthe escapingelectrons,theto-
grates the cross sections of the ionization reaction over
tal cross section and the energy differential cross section
the velocity space. This Boltzmann equation is coupled
in the hydrogen atom. We name exterior complex scal-
to an evolution equation for the electrical field. Because
ing [7],timedependent close coupling [8],HRW-SOW [9]
additional electrons are created, the local chargedensity
and convergent close coupling [10].
changes the electrical field through the Poisson law
When the electron hits a molecular system instead of
anatom, the physics is complicatedby the extra degrees ∇·E(x,t)=q(x,t),
of freedom. The cross sections now depend on both the
orientation and internuclear coordinates of the molecule
where q(x,t)=(n −n )e/q is the charge distribution.
+ e 0
atthe momentofelectronimpact,as is seeninprocesses
Here,n representsthenumberofions,n isthenumber
+ e
wheretwo electronsareejected frommolecules afterit is
of electrons and q a unit of charge.
0
hit by a photon [11, 12]. Therefore, there will be some
random terms in the reaction cross section, which will We now connect the change in electricalfield with the
need to be included in a realistic microscopic model. In change in the charge distribution. We have
this paper, we will model the microscopic interactions
∂q(x,t)
using a Boltzmann equation, which is still deterministic; +∇·j(x,t)=0
extensionsthataccuratelyaccountforrandomeffectswill ∂t
be treated in future work. However, we note that the
inwhichj(x,t) is the chargeflux. This leads to anequa-
coarse-grainedtime-stepperapproachthatisusedinthis
tion for the evolution of the electrical field,
work has already been applied successfully to study sys-
tems with stochastic effects [20]. ∂E(x,t)
+j(x,t)=0. (2)
∂t
B. Review of the physics of streamer fronts. Since we assume that the ions are immobile, the flux
j(x,t) is solely determined by the one-particle distribu-
tion function f(x,v,t) of the electrons.
Ebertetal.[1]introducedtheminimalstreamermodel.
It consists of two coupled non-linear PDEs: a reaction- Our extension of the minimal streamer model is now
convection-diffusion equation for the evolution of the the coupled evolution of eq. (1) and (2). Note that the
electron density and an equation that relates the change set of coupled equations is very similar to the Wigner-
intheelectricalfieldtothechargeflux. Theelectronden- Poisson problem [15] used to model electron transport
sityevolvesbecauseofthedriftduetotheelectricalfield, through diodes.
the electron diffusion and the ionization reaction, which
is formulated in the Townsend approximation. The re-
actionrate is then givenby an exponentialthat depends
on the strength of the local electrical field. The evolu-
tion of the electrical field is determined by Poisson’s law
of electrostatics where the field changes because of the III. LATTICE BOLTZMANN DISCRETIZATION
chargecreationby the ionizationreaction. This minimal
streamer model exhibits both negatively and positively
charged fronts. The first moves in the direction of the Togetherwiththeimpactionizationcrosssections,the
electricalfield, while the positively chargedmoves in the coupled equations (1) and (2) are a non-linear integro-
oppositedirectionandcanonlypropagatebecauseofthe differential equation coupled to a scalar partial differen-
electron diffusion and the ionization reactions. Each of tial equation for the electrical field. This equation in its
these fronts appears as a one-parameter family of uni- fulldimensionis hardto solve,bothanalyticallyandnu-
formly translating solutions (since any translate of the merically. Asafirststep,welookatthe one-dimensional
wave is also a solution). streamer fronts of the Boltzmann equation in the lattice
In our extension of the Ebert model, we replace the Boltzmann discretization.
4
A. Discretization of the Boltzmann equation where we have introduced the shorthand f (x,t) for
i
f(x,v ,t) and Ω (x,t) for Ω(x,v ,t) with i∈S.
i i i
In this section, we discretize the one-dimensional
Boltzmann equation (1). The distribution functions
B. The collision term
f(x,v,t) are discretized on a lattice in space, velocity
and time. The grid spacing is ∆x in space and ∆t in
time. The velocity grid of v is chosen such that the dis- The collision term consists out of two parts
i
tance traveled in a single time step, v ∆t, is a multiple
i
∆tΩ =Ωdiff+Ωreaction (6)
of the grid distance ∆x, or in short i i i
∆x wherethefirsttermwillmodeltheelectrondiffusion,the
v =i , with i∈S.
i second term the ionization reactions and the third term
∆t
the influence of the external force. Note that we also
Typically, only a small set of discrete velocities is used. incorporated ∆t in the notation. We will now discuss
A discretization with three grid points on the velocity the two terms individually.
grid has a set S = {−1,0,1} and is called a D1Q3 The first term models the electron diffusion as a BGK
model. A discretization with five grid points has S = relaxation process [21]. In this approximation, it as-
{−2,−1,0,1,2}andisaD1Q5model. Thesizeoftheset sumed that the distribution is attracted to a local equi-
S is denoted by m. We will also denote ci = vi∆t/∆x librium distribution function feq,
i
for the dimensionless velocity.
1
Note that,for easeofnotation,wewillalsousethe set Ωdiff =− (f −feq), (7)
S to index matrices. For example, the result of a linear i τ i i
operator A working on a vector vj∈S will be denoted as with feq(x,t) the equilibrium distribution for electron
A v , where both the indices i and j are in S. i
i∈S ij j diffusion. In the five speed model, we choose
The means that the A−2,−2 matrix element with S =
P
{−2,−1,0,1,2} is the matrix element in the upper left feq =weqρ with weq ={0,1/4,2/4,1/4,0} (8)
i i i
corner of the matrix.
Westartfromthecontinuousequationforthedistribu- and ρ(x,t) the electron density. This choice of equi-
tionfunction in the discrete point(x+v ∆t,v )in phase librium weights conserves the number of electrons, but
i i
space at time t. In the absence of external forces, the doesnotconservemomentumasatraditionalfluidwould
Boltzmann equation in this point reads do. Indeed, the electrons diffuse because they randomly
change their direction during the elastic collisions with
∂f(x+vi∆t,vi,t) ∂f(x+vi∆t,vi,t) the much heavier neutral molecular particles. We fur-
+v
∂t i ∂x ther chose the weights such that there are no fast parti-
=Ω(x+v ∆t,v ,t) (3) clesunderdiffusiveequilibrium. Therelaxationtimeτ is
i i
related to the electron diffusion coefficient
A discrete lattice Boltzmann equation is now obtained
1 D ∆t
by replacing the time derivative with anexplicit forward τ = + . (9)
difference, the introduction of an upwind discretization 2 i∈Sc2iwieq ∆x2
of the convection term and a downwind discretization of
Notethat,intheliteratuPre,therelaxationtimeτ isoften
the collision term Ω(x+v ∆t,v ,t) and replace it with
i i
characterizedby its inverse ω =1/τ.
Ω(x,v ,t) [14],
i The second term in (6) is the reaction term Ωreaction
i
f(x+v ∆t,v ,t+∆t)−f(x+v ∆t,v ,t) that is modelled with a m×m matrix R
i i i i
∆t
Ωreaction =∆t R f , (10)
f(x+v ∆t,v ,t)−f(x,v ,t) i ij j
i i i
+vi v ∆t =Ω(x,vi,t). (4) jX∈S
i
which represents the velocity dependent reaction rates
Note that this discretization of the spatial derivative
and allows us to select between slow and fast particles.
becomes less accurate for the largest speeds in the set
For the five speed model, we choose a reaction matrix
S. Indeed, in the five speed model, for example, the
largest speed is v±2 = ±2∆x/∆t, and the convec- −R 0 0 0 0
tion term will be calculated from the difference between
R 0 0 0 R
f(x+v±2∆t,v±2,t)andf(x,v±2,t),whichis2∆xapart. R= 0 0 0 0 0  (11)
The discretization error is then proportionalto 2∆x.
R 0 0 0 R
 
Equation (4) reduces to  0 0 0 0 −R
 
 
f (x+v ∆t,t+∆t)−f (x,t)=∆tΩ (x,t), (5) that describes how the reactioncrosssections depend on
i i i i
5
the velocities ofthe particles. Eachtime step ∆t, a frac- The Galerkin condition leads to the linear system
tion R of the particles with speed v±2 will collide and
cause a ionization reaction. The reaction rate R is cho- m 0 α 0 β a vt · ∂f
0 0 ∂v
sen to match the height of the cross section, see figure 0 α 0 β 0 a vt · ∂f
1. Since the colliding fast particles transfer their energy  α 0 β 0 γa1=v1t · ∂∂vf  (12)
2 2 ∂v
tfooret,htehbeonuunmdbeelrecotfropna,rttichleeys wwiitlhl lsopoeseedevn±e2rgdy.imiTnhisehrees-  β0 β0 γ0 γ0 0δ aa34 vv34tt ·· ∂∂∂∂vvff 
with a rate−R∆t andwe haveR−2,−2 =R+2,+2 =−R.     
Atthesametime,the numberofslowelectronsincreases where α = v2,β = v4, γ = v6 and δ =
i∈S i i∈S i i∈S i
becauseboththeimpactelectronandtheionizationelec- v8.
tronemergeasslowparticleswithspeedv±1. Becauseof i∈S i P P P
the Coulomb repulsion, the two slow electrons are more PTo calculate the right-hand side of (12) we make a
likelytoemergeinopposite directionsandwechoosethe detour around the continuous representation. We note
rates such that one electron emerges with speed of v−1 that
and the other with v+1. ∂f +∞ ∂f(x,v,t)
This choice of model parameters ensures that the en- vt· = vl dv, (13)
l ∂v ∂v
ergy balance during the ionization reaction is not vio- Z−∞
lated. As discussed in section IIA, the energy of the
where l ∈ {0,1,...,m−1}. Because of our particular
incoming electron is larger that the sum of the energies
choice of basis vectors and the fact that there are no
of the escaping electrons because some of the impact en-
particles with infinite velocities, we have that
ergy coversthe binding energy of the bound electron. In
the above model, a single ionization reaction transforms +∞ ∂f(x,v,t) +∞ ∂vl
one electron with speed v into two electrons with, re- vl dv + f(x,v,t) dv
+2
∂v ∂v
spectively, speed v+1 and speed v−1. The energy of the Z−∞ Z−∞
impact electron is mv2 /2 = 2m∆x/∆t, while the sum =vlf(x,v,t)|+∞ =0 (14)
+2 −∞
of the escaping electrons is merely 2mv2/2 = m∆x/∆t,
1
where m is the mass of the electron. So a portion or, in other words,
m∆x/∆toftheimpactenergycoversthebindingenergy.
For our model, this is half of the initial impact energy; vt· ∂f = −i +∞f(x,v,t)vl−1dv
for more general problems other values are possible. l ∂v −∞
Z
= −l vl−1f (15)
j j
j∈S
X
C. External Force
With the help of
We now derive a discretization of the E∂fi term that
∂v 0 0 0 0 0
modelstheexternalforceintheBoltzmannequation. We
−1 −1 −1 −1 −1
start by expanding N = −2v−2 −2v−1 −2v0 −2v1 −2v2 , (16)
∂f −3v−22 −3v−12 −3v02 −3v12 −3v22
∂v =a0v0+a1v1+...+am−1vm−1, −4v−23 −4v−13 −4v03 −4v13 −4v23
 
where V = {v0,v1,...,vm−1} forms a linear indepen- we can now define
ad0e,nat1,s.e.t.,oafmv−e1ctboyrseninforRcimng. thWeeGafilnerdkitnhceoncdoietffiiocni.ents 11 vv−−21 vv−−2122 vv−−2133 vv−−2144 m0 α0 α0 β0 β0 −1
∂f m−1 V=1 v0 v02 v03 v04  α 0 β 0 γ N,
− aivi ⊥V 1 v1 v12 v13 v14  0 β 0 γ 0
∂v ! 1 v v 2 v 3 v 4  β 0 γ 0 δ 
i=0  2 2 2 2  
X
   (17)
In the current paper, we choose a particular set of vec- and calculate the external force term as
tors in V, namely the polynomials {1,v,v2,...,vm−1}
discretized in the points vi. For the five speed example E(x,t)∂fi(x,t) =E(x,t) V f (x,t), (18)
ij j
the vectors are, besides their powers of ∆x/∆t, ∂v
j∈S
X
1 −2 4 −8 16
where the elements of S denote matrix elements.
1 −1 1 −1 1
V =1, 0 ,0, 0 , 0  Fromeq. (18),itisclearthatwecanincludetheexter-
1  1  1  1   1  nalforceasanadditionalcollisiontermintheright-hand
1  2  4  8  16 side of the lattice Boltzmann equation.
         
         
6
D. Flux scopic variable of interest. The transformation between
distribution functions f and moments ̺ can be written
i l
as a matrix transformation M. In the five speed model,
The evolution of the electrical field E(x,t) is deter-
this matrix is
mined by the net flux j(x,t) of electrons as expressed in
(2). Wediscretize(2)onastaggeredgridwithgridpoints
1 1 1 1 1
halfwaybetweenthegridpointsofthelatticeBoltzmann
−2 −1 0 1 2
model. The flux is defined as the number of particles M = 4 1 0 1 4  (23)
that move between grid points (pass through an inter-
−8 −1 0 1 8
 
face)withinasingletime step. Forthefive speedmodel,  16 1 0 1 16
we have  
 
such that ̺ = M f and f = m−1 M−1 ̺ .
j(x+∆x/2, t) =f1(x+∆x,t)−f−1(x,t) l i∈S li i i l=0 il l
An evolution law for ̺ (x,t) is now easily constructed
+f2(x+∆x,t)−f−2(x,t) (19) by the followingPsequencel: first transPform̺(cid:0)into f(cid:1) (x,t)
l i
+f2(x+2∆x,t)−f−2(x−∆x,t) usingM−1,thenusethelatticeBoltzmannequation(20)
toevolvef (x,t)tof (x,t+∆t)and,then,transformback
i i
to the moments ̺ (x,t+∆t).
l
E. Coupled equations It has been observed phenomenological that the ion-
ization wave can approximately be described by a PDE
The coupled equations (1) and (2) for the evolutionof in the density. This suggests that, in practice, the evo-
theelectrondistributionfunctionsandtheelectricalfield lution of these moments is rapidly attracted to a low
is now, after discretization, dimensional manifold described by the lowest moment
̺ (x,t), which is the density. The higher ordermoments
0
f (x+v ∆t,t+∆t)−f (x,t)= have then become functional of this density and the dy-
i i i
namics of the system can effectively be described by the
1
− τ fi(x,t)−fieq(x,t) + ∆tRijfj(x,t) evolution of this macroscopic moment.
j∈S Ingeneral,however,itis veryhardtofindanalyticex-
(cid:0) (cid:1) X
(E(x−∆x/2,t)+E(x+∆x/2,t)) pressionsforthislowdimensionaldescriptionintheform
− ∆tV f (x,t)
2 ij j ofaPDEwithoutmakingcrudeapproximations. Forthe
jX∈S problem at hand, we illustrate these difficulties in this
(20)
section where we apply the Chapman-Enskog expansion
and derive a macroscopic PDE in terms of electron den-
sity. Onlyafterdroppingseveralcouplingterms,aclosed
∆x ∆x ∆x
E(x+ ,t+∆t)=E(x+ ,t)−∆tj(x+ ,t), (21) PDE is derived.
2 2 2
The model discussed in the previous sections can be
where j(x+∆x/2,t) is calculated from (19). The equa- summarized by the lattice Boltzmann equation
tionsarecoupledbecausetheelectricalfieldappearsasan
external force in the first equation, while the flux drives fi(x+ci∆x,t+∆t)−fi(x,t)=
the evolution of the electrical field in the second equa- 1
− (f (x,t)−feq(x,t))+ A f (x,t), (24)
tion. This coupling makes the evolution of the system τ i i ij j
non-linear. Xj∈S
for ∀i ∈ S. Here, A can be the reaction term R , a
ij ij
force term V , or a combination of both.
ij
IV. A PDE MODEL THROUGH
A second order Taylor expansion of the term f (x+
i
CHAPMAN-ENSKOG EXPANSION
c ∆x,t+∆t) in (24) around f (x,t) leads to
i i
Themodel(20)–(21)evolvestheelectricalfieldE(x,t) ∂f ∂f c2∆x2 ∂f
c ∆x i +∆t i + i i
and the distribution functions fi∈S(x,t) from t to t + i ∂x ∂t 2 ∂x2
∆t simultaneously. Alternatively, the evolution of the ∂2f ∆t2∂2f
i i
distribution functions can be rewritten in terms of the +ci∆x∆t∂x∂t + 2 ∂t2
corresponding (dimensionless) velocity moments defined
1
as =− (f −feq)+ A f , ∀i∈S.(25)
τ i i ij j
j∈S
̺ (x,t)= clf (x,t), (22) X
l i i
i∈S We then expand fi in terms of increasingly higher order
X
contributions as follows
wherel ∈{0,1,...,m−1}. Thezerothmoment̺ (x,t)
l=0
correspondstotheelectrondensityρ(x,t),i.e.themacro- f =f(0)+ǫf(1)+ǫ2f(2)+... (26)
i i i i
7
with ǫ a small tracer parameter. In fluid dynamics, ǫ and are found as the the solution of the linear system
typically refersto the Knudsennumber. The spatialand
time derivatives are scaled respectively as (1−τA )f(0) =feq. (30)
ij j i
j∈S
∂ ∂ ∂ ∂ ∂ ∂ X
= +ǫ +ǫ2 and =ǫ , (27)
∂t ∂t ∂t ∂t ∂x ∂x Since feq = w ρ (8) only depend on the density, we find
0 1 2 1 i i
f(0) =w(0)ρ with the weights w(0) defined by
where we explicitly presume that a zeroth order time i i i
scale t is present in the system. As we will show later
on, thi0s scale corresponds to the observed exponential wi(0) = (1−τA)−1 wjeq. (31)
ij
growth of the electron density. Xj∈S(cid:16) (cid:17)
Becauseofthemultipletimescalest ,t andt ,allthe
0 1 2
terms in the expansion will couple to all the time scales,
whichcomplicatesthederivationofaneffectiveequation.
Since the matrix A can include the ionization reac-
In our search for a reduced second order PDE model, tion that does not conserve the particle number, the
we only keepthe terms up to secondorderin ǫ2. For the sum of the weights w(0) is not necessarily equal
same reason, we also drop the second derivative w.r.t. i∈S i
to one. This means that f(0) 6= ρ. We choose
time from (25). Substitution of (26) and (27) into (25) P i∈S i
leads to to rescale the weights w(0) with a normalization factor
i P
N = w(0) = (1−τA)−1 weq such that
ǫci∆x∂∂fix(0) +ǫ2ci∆x∂∂fix(1) +ǫ2c2i∆2x2∂∂2fxi(20) Pf(i0∈)S=iρ. WPithi,jt∈hSis(cid:16)rescaling w(cid:17)eijfindj the zeroth
i∈S i
∂f(0) ∂f(1) ∂f(2) order term of the Chapman-Enskog expansion
+∆t i +ǫ∆t i +ǫ2∆t i P
∂t0 ∂t0 ∂t0 f(0) =w(0)ρ= 1 (1−τA)−1 weqρ (32)
+ǫ∆t∂∂fti(0) +ǫ2∆t∂∂fti(1) +ǫ2∆t∂∂fti(0) i i N jX∈S(cid:16) (cid:17)ij j
1 1 2
∂2f(0) ∂2f(1)
+ǫc ∆x∆t i +ǫ2c ∆x∆t i
i i
∂x∂t ∂x∂t
0 0 The rescaling, however, forces us to reconsider equa-
+ǫ2ci∆x∆t∂2fi(0) tion (30) because fi(0), as defined above, fails to be a
∂x∂t1 solution. Still, we keep our f(0) of (32) as our zeroth-
i
= −1 f(0)+ǫf(1)+ǫ2f(2)−feq order term and find for the evolution of the system at
τ i i i i time scale t
0
(cid:16) (cid:17)
+ A f(0)+ǫf(1)+ǫ2f(2) (28)
ij j j j ∂f(0) 1 1
Xj∈S (cid:16) (cid:17) ∆t ∂ti =−τ (1−τAij)fj(0)+ τfieq
0 j∈S
We will now group the terms order by order and derive X
expressionsfor f(0), f(1) and f(2) and the corresponding =−1 (1−τA ) 1 (1−τA)−1 feq+ 1feq
evolution equatioins foir ρ at thei different time scales. τ Xj∈S ij kX∈S N (cid:16) (cid:17)jk k τ i
Wewillusethefactthatif∆tand∆x2 areofthesame 1 1
orderofmagnitude—whichisthecaseforourexamples = τ 1− N fieq
—the terms that havefactorsas ∆t∆x areeffectively of (cid:18) (cid:19)
order∆x3 and canbe neglectedcomparedto terms with =αfieq∆t, (33)
∆t or ∆x2.
where the growth factor is
α=(1−1/N)/(τ∆t). (34)
1. Zeroth order contribution
Fromexpansion(28),wecollectthezerothorderterms Summation of (33) over the set S leads to the zeroth
order PDE for the evolution of ρ
(0)
∂f 1
∆t ∂ti0 =−τ (cid:16)fi(0)−fieq(cid:17)+Xj∈SAijfj(0). (29) ∂∂tρ0 =αρ. (35)
(0)
We choose f such that the right-hand side of (24) is This is a growth equation if α is positive, which is the
i
zero; these distributions will not evolve on time scale t case for the ionization reaction.
0
8
2. First order contribution 3. Second order contribution
Finally, we derive the second order evolution. We col-
lect from (28) the second order terms and find
Toderivethe firstorderequation,wecollectthe terms
∂f(1) c2∆x2∂2f(0) ∂f(2)
that are first order in ǫ from (28) c ∆x i + i i +∆t i
i ∂x 2 ∂x2 ∂t
0
∂f(0) ∂f(1) ∂f(1) ∂f(0) ∂2f(1)
ci∆x ∂ix +∆t ∂ti0 +∆t ∂ti1 +∆t ∂ti2 +ci∆t∆x∂x∂it0
∂f(0) ∂2f(0) ∂2f(0) 1
+∆t ∂ti1 +ci∆t∆x∂x∂it0 +ci∆t∆x∂x∂it1 =−τ j∈S(1−τAij)fj(2) (41)
1 X
=− (1−τA )f(1). (36)
τ jX∈S ij j The terms ∂2fi(1)/(∂x∂t0), ∂2fi(0)/(∂x∂t1), and
∂f(1)/∂t (when replacing f(1) by (37)) are of order
We dropthe termci∆t∆x∂2fi(0)/(∂x∂t0)becauseitisof ∆ti∆x, w1hich is smaller thani ∆x2 for our parameter
order ∆t∆x, which is smaller than the other terms that (2)
settings. We also neglect ∆t∂f /∂t because it can be
i 0
are first order in ∆t or ∆x. The second term we neglect
shown, again a postiori, that it is of order ∆t∆x.
is ∆t∂f(1)/∂t because we will show below, a postiori,
i 0 The second order expansion term then becomes
that it is also of order ∆t∆x.
∂f(1) c2∆x2∂2f(0)
c ∆x i + i i
i ∂x 2 ∂x2
We now have ∂f(0) 1
+∆t i =− (1−τA )f(2), (42)
c ∆x∂fi(0) +∆t∂fi(0) = −1 +A f(1) ∂t2 τ Xj∈S ij j
i ∂x ∂t τ ij j
1 jX∈S(cid:18) (cid:19) When replacing f(1) and f(0) by (40) and (32), we get
i i
what leads to a first order term
c ∆x (−1/τ +A)−1 w(0)(c ∆x−C∆t) ∂2ρ
∂f(0) ∂f(0) i ij j j ∂x2
fi(1) =Xj∈S(cid:16)(−1/τ +A)−1(cid:17)ij cj∆x ∂jx +∆t ∂tj1 !. jX∈S(cid:16) +c2i∆(cid:17)x2wi(0) ∂2ρ +∆tw(0) ∂ρ
(37) 2 ∂x2 i ∂t
2
We now see that it is justified to neglect the term
1
∆t∂fi(1)/∂t0 in (36) because it is of order ∆t∆x. Using =−τ (1−τAij)fj(2)(43)
(32) and f(1) = 0 (the latter because f = Xj∈S
i∈S i i∈S i
f(0)+ǫf(1)+ǫ2f(2) = ρ and f(0) = ρ), and the expression for the second order term becomes
i∈S iP i i i∈SPi
Pwe find(cid:16)the following PDE a(cid:17)t time scalePt1 for the evolu- ∂2ρ
tion of the system f(2) = B c B (c ∆x−C∆t)∆xw(0)
i ij j jk k k ∂x2
j∈S j∈S
∂ρ ∂ρ X X
+C =0, (38) ∆x2 ∂2ρ
∂t ∂x + B c2w(0)
1 ij j j 2 ∂x2
j∈S
where the advection coefficient c equals X
∂ρ
+ B w(0)∆t , (44)
i,j∈S (−1/τ +A)−1 cjwj(0) ∆x Xj∈S ij j ∂t2
ij
C = . (39)
P i,j∈S(cid:16) (−1/τ +A)−1(cid:17) wj(0) ∆t where we use Bij = (−1/τ +A)−1 ij. If we define
ij
P (cid:16) (cid:17) (cid:0) (cid:1)
With the help of (38), the first order contribution (37)
can be written alternatively as D =  BijcjBjk(ck∆x−C∆t)∆xwk(0) (45)
i,j,k∈S
X
f(1) = (−1/τ +A)−1 w(0)(c ∆x−C∆t)∂ρ. 
i Xj∈S(cid:16) (cid:17)ij j ∂x + Bijc2jwj(0)∆x2/2/− Bijwj(0)∆t,
(40) k,i∈S i,j∈S
X X
  
9
we obtain the evolution at time scale t (because 0.05
2
f(2) =0)
k 0.04
P ∂ρ ∂2ρ or
=D , (46) ct 0.03
∂t2 ∂2x hfa
t
w 0.02
where D acts as a diffusion coefficient. o
r
g
Wearenowinthepositiontocombinetheevolutionat 0.01
the different timescales t , t and t into a single PDE.
0 1 2
Because the matrix A of the model equation (24) con- 0
tains both the reaction term and the external force, the 0 0.2 0.4 0.6 0.8 1
transportcoefficientsα,C andDwilldependonthelocal Electric field |E|
electrical field E(x,t). For our example the dependence 0
of the transport coefficients is shown in figure 2, we find
that D hardly depends on the strength of the field and nt −0.25
can be set equal to the electron diffusion D used in (9) cie
ffi
to define the relaxation of the lattice Boltzmann model. e
In a similar way,we find that c, the transport coefficient co −0.5
n
o
oftheadvectionterm,isapproximatelyequaltothe−E, ti
c
the local electrical field that causes the drift. Only the dve −0.75
growthfactorα,definedin(34),dependsonthestrength a
of the localelectricalfield. With the help of these obser-
−1
vations we get the coupled PDE
0 0.2 0.4 0.6 0.8 1
∂ρ ∂ρ ∂2ρ Electric field|E|
= α(E(x,t))ρ+E(x,t) +D 1.004
∂t ∂x ∂x2
∂E ∂ρ
∂t = −E(x,t)ρ−D∂t (47) ent 1.003
ci
ffi
fortheevolutionofE(x,t)andρ(x,t). Thesecondequa- e
co 1.002
tion expresses the flux with the help of the transport n
o
coefficients. si
u
These coupled equations are similar to minimal diff 1.001
streamer model of Ebert, Van Saarloos and Caroli [1],
except that the growth rate is now defined by (34) in- 1
stead of the Townsend approximation. In figure 2, we 0 0.2 0.4 0.6 0.8 1
illustrate how the growth coefficient depends on the lo- Electric field |E|
cal electrical field and compare with a Townsend ap-
proximation. We find that a Townsend reaction term Figure 2: Top: The growth factor α(E) (solid) of the
0.111·|E|exp(−1/|E|)approximatelydescribesasimilar PDE model derived from the lattice Boltzmann model. The
growth term as the PDE model derived from the lattice growth depends on the strength of the local electrical field
Boltzmann model. and is similar to the Townsend approximation with 0.111·
|E|exp(−1/|E|)(dashed). We have a reaction rate R = 100
and model parameters given in section VII. Middle: Thead-
vectioncoefficientCisequaltotheexternalfield−E. Bottom:
ThediffusioncoefficientDonlychangeswiththeexternalfield
4. Traveling wave solutions in the fourth significant figure.
The system (47) is non-linear and it is well-known
thatithasaone-parameterfamilyoffrontsolutionsthat
translateuniformlywithaspeedc[2]. Thereisaminimal
speed c∗ that is usually found by looking at the asymp- where the transport coefficients have become constants.
totic region → +∞. In this limit, the electrical field In a co-moving coordinate frame that travels with the
becomes constant and is denoted by E+ — the same same speed c along the x-axis, we define ξ = x − ct.
notation as in [1] — and the equation for the electron The PDE (48) becomes stationary and the solution in
density becomes the asymptotic region fits
∂ρ =α(E+)ρ+E+∂ρ +D∂2ρ, (48) 0=α(E+)ρ(ξ)+(c+E+)∂ρ(ξ) +D∂2ρ(ξ). (49)
∂t ∂x ∂x2 ∂ξ ∂ξ2
10
ThelatterisasecondorderODEthatcanbetransformed Algorithm 1 Constrained runs scheme for LBM
into two coupled first order ODEs by denoting ∂ρ/∂ξ as initialize f[0] =weqρ(x,t) ∀i∈S
v and ρ as u . The system of coupled equations is i i
repeat
f[k+1] =LBM(f[k]) a single LBM step
u 0 1 u
(cid:18)vξξ(cid:19)=(cid:18)−α(DE+) −c+DE+ (cid:19)(cid:18)v(cid:19) (50) ̺ρ[[kk++11]] ==Mρ(xf,[kt)+1] mreaseptitnhteodmeonmsiteynts
whereu andv denotederivativesof,respectively,uand f[k+1] =M−1̺[k+1] map into distributions
ξ ξ
until convergence heuristic
v to ξ. The matrix has two eigenvalues
−(c+E+)± (c+E+)2−4Dα(E+) Table I: Lifting. The constrained runs algorithm computes
λ± = . the distribution functions fi(x,t) corresponding to a given
2D
p densityρ(x,t). Thesuperscriptkindicatestheiterationnum-
There are two cases, if (c+E+)2 < 4Dα(E+) the two ber.
eigenvalues are complex, otherwise, they are real.
The electrondensityinthe asymptotic regionis nowa A. Lifting
linear combination of two exponentials
SincetheelectricalfieldE(x,t)isthesameinboththe
lim ρ(x)=Aeλ+x+Beλ−x. (51) lattice Boltzmann and the macroscopic model, we can
x→+∞
ignore it for the discussion of the lifting and restriction
When the two eigenvalues are complex, the asymptotic operators. In the lifting step, the particle distribution
density is oscillating and can becomes negative. This is functions are initialized starting from the initial density
unphysicalbecausewecannothaveanegativenumberof
particles and it is concluded that the speed c has to be µ:Rn 7→Rn×m :ρ(xj,t)7→fi(xj,t)
above a minimal speed
with i ∈ S, m the number of speeds in S, and j ∈
c>c∗ =c(E+)+2 Dα(E+), (52) {1,2,...,n}denotingthediscretespatialgridpoints. Be-
causelifting is a one-to-manymapping problem,itis the
tokeepbotheigenvaluesreal. Noptethatbotheigenvalues mostcriticalstepinthecoarse-grainedtime-stepper. We
λ+ and λ− coalesce at the critical speed c=c∗. use the constrained runs scheme, an algorithm proposed
in [17] in the context of singularly perturbed systems.
Here, it is wrapped around a single time step ∆t of the
lattice Boltzmann model [4].
TheprocedureisgivenintableI. Startingfromanini-
tialguessρ(x,t) forthe density, weobtaininitialguesses
forthedistributionfunctionsusingtheBGKequilibrium
V. THE COARSE-GRAINED TIME-STEPPER
weights (8). This choice determines the initial guess for
the higher order moments through the transformation
matrix M, see equation (23). (In principle, the initial
In this section, we describe an alternative way to per-
guesses for the higher order moments can be chosen ar-
formthe analysisofthe macroscopicbehaviorofthe sys-
bitrarily; the scheme is designed precisely to convergeto
tem. It is based on the work of Kevrekidis et al. [5] who
thecorrectvalueofthesemomentsforthegivendensity.)
developed a coarse-grainedtime stepper (CGTS), which
We then performthe followingiteration. First, we use
is an effective evolution law for the density. This evolu-
tion law F is not an analytic expression such as a PDE, the lattice Boltzmann model to evolve fi[∈k]S from t to
but the following sequence of computational steps: 1) t+∆t. The result is transformed back into the moment
lifting, 2) simulation and 3) restriction, denoted by the representationbyamatrixmultiplicationwithM,which
operators µ, LBM and M respectively (figure 3). Note gives us ̺[k+1]. Next, the zeroth moment of the vector
that the simulation time ∆T is in general a multiple of ̺[k+1] is reset to the initial value ρ(x,t). Transforming
∆t, the lattice Boltzmann time step. Formally, this is this modified moment vector ̺[k+1] back into distribu-
written as tion functions gives us the next f[k+1]. We repeat this
i∈S
iterationuntilthehigherordermomentshaveconverged.
U(x,t+∆T) = F(U(x,t),∆T) Theconvergencebehavioroftheconstrainedrunsalgo-
= M(LBM(µ(U(x,t)),∆T)), (53) rithm from table I is analyzed in [4] for one-dimensional
reaction-diffusion lattice Boltzmann models with S =
where we have introduced U(x,t) = (ρ(x,t),E(x,t)) as {−1,0,1} (D1Q3 stencil) and a density dependent re-
a shorthand notation. The time-stepper F evolves the actionterm. Forsuchsystems,the algorithmisuncondi-
macroscopic density ρ(x,t)= f (x,t) and the elec- tionallystableandconvergesuptothefirstordertermsin
i∈S i
trical field E(x,t) from time t to t+∆T. theChapman-Enskogexpansionofthedistributionfunc-
P