# Computing the Paschen curve for argon with speed-limited particle-in-cell
simulation

Joseph G. Theis<EMAIL_ADDRESS>Center for Integrated Plasma
Studies, University of Colorado, Boulder, Colorado 80309, USA  Gregory R.
Werner  Center for Integrated Plasma Studies, University of Colorado, Boulder,
Colorado 80309, USA  Thomas G. Jenkins  Tech-X Corporation, 5621 Arapahoe
Avenue Suite A, Boulder, Colorado 80303, USA  John R. Cary  Center for
Integrated Plasma Studies, University of Colorado, Boulder, Colorado 80309,
USA  Tech-X Corporation, 5621 Arapahoe Avenue Suite A, Boulder, Colorado
80303, USA

###### Abstract

Upon inclusion of collisions, the speed-limited particle-in-cell (SLPIC)
simulation method successfully computed the Paschen curve for argon. The
simulations modelled an electron cascade across an argon-filled capacitor,
including electron-neutral ionization, electron-neutral elastic collisions,
electron-neutral excitation, and ion-induced secondary-electron emission. In
electrical breakdown, the timescale difference between ion and electron motion
makes traditional particle-in-cell (PIC) methods computationally slow. To
decrease this timescale difference and speed up computation, we used SLPIC, a
time-domain algorithm that limits the speed of the fastest electrons in the
simulation. The SLPIC algorithm facilitates a straightforward, fully-kinetic
treatment of dynamics, secondary emission, and collisions. SLPIC was as
accurate as PIC, but ran up to 200 times faster. SLPIC accurately computed the
Paschen curve for argon over three orders of magnitude in pressure.

##  I Introduction

Townsend discharge is the conduction of electric current in a gas due to a
cascade of ionizing collisions.Townsend (1900) In typical discharge
experiments, the gas is enclosed between the electrodes of a capacitor, and
there is an initial source of free electrons. The electrons accelerate in an
applied electric field and ionize the neutral atoms, creating more electrons,
which ionize more atoms. The ions accelerate towards and impact the cathode,
where they can produce secondary electrons. When the combined effects of
ionizing collisions and ion-induced cathode emission yield more electrons at
the cathode than initially seeded, a self-sustaining discharge is triggered.
The voltage difference at which this occurs is the breakdown voltage.

The breakdown voltage $V_{b}$ of a gas-filled capacitor is (to a good
approximation) a function of $pd$, where $p$ is the pressure of the gas and
$d$ is the distance between the electrodes. A Paschen curve is the
relationship between $V_{b}$ and $pd$ for a given gas.Paschen (1889) Paschen
curves have a minimum $V_{b}$ at some intermediate $pd$ and higher $V_{b}$ at
the extremes (this shape is shown later in Fig. 2). This characteristic shape
is due to the competition between ionizing collisions, which multiply the
number of free electrons, and non-ionizing collisions, which prevent electrons
from gaining enough energy to ionize atoms. At low $pd$, the ionization mean
free path is too long for a sufficient number of ionizations to occur during
crossing. At high $pd$, the non-ionizing mean free path is too short for
electrons to reach the ionization energy. Paschen curves will be developed
further in §II.1.

Electrical discharge and breakdown have been extensively studied both for a
multitude of applications, including lighting,Sugiura (1993) plasma
processing,Lieberman and Lichtenberg (2005) and lasers,Belasri, Boeuf, and
Pitchford (1993) as well as to avoid detrimental effects, including damage to
electrical devices,Gao  _et al._ (2003) explosions in chemical processing
facilities,Gillman and Le May (2007) and damage to aerospace systems.Ashwell
_et al._ (1966) Following Townsend’s description of the mechanism behind gas
discharge,Townsend (1900) various models and simulation methods have been
employed to study the phenomena.Phelps and Petrovic (1999) Analytical models
have proven effective for computing breakdown and steady-state current
amplification in simple geometries;Phelps and Petrovic (1999) however,
understanding the particle dynamics of discharge and handling complex
geometries requires simulation.

Particle-in-cell (PIC) and fluid methods have been the primary simulation
techniques applied to model gas discharges.Boeuf and Pitchford (1991);
Lieberman and Lichtenberg (2005) Fluid approaches evolve the mean properties
of a particle species and can therefore handle situations involving many
particles and mean free paths that are small compared with the system size.
However, the fluid approaches assume Maxwellian particle velocity
distributions, while discharges often have non-Maxwellian distributions.Boeuf
and Pitchford (1991)

PIC methods evolve the distribution function $f(x,v,t)$ in time via the method
of characteristics and, therefore, can simulate arbitrary distribution
functions.Birdsall (1991) Collisions can be included in PIC using the PIC
Monte-Carlo Collisions framework (PIC-MCC).Birdsall (1991); Bird (1994); Nanbu
(2000) PIC-MCC will be discussed further in §III. PIC-MCC requires that the
smallest mean free path be resolved by the grid cell size $\Delta x$. Further,
in PIC, the timestep $\Delta t$ is limited according to $\Delta
t\lesssim\Delta x/v_{\rm max}$, where $v_{max}$ is the maximum particle speed.
In gas discharge simulations, the large mass difference between the ions and
the electrons results in ion timescales much greater than the timestep imposed
by the maximum electron speed. This makes PIC-MCC simulations computationally
slow. Because PIC-MCC is slow, hybrid methods combining PIC and fluid
approaches (e.g., with fluid electrons but PIC ions) have also been applied to
discharge,Boeuf and Pitchford (1991) but they fail to capture the full
velocity distribution as PIC does.

In this work, we apply the speed-limited particle-in-cell (SLPIC) simulation
methodWerner  _et al._ (2018) to electrical breakdown in argon. SLPIC is a
time-domain algorithm that limits the speed of the fastest particles in the
simulation to enable larger time steps and, therefore, faster computing times.
The SLPIC algorithm facilitates a straightforward, fully-kinetic treatment of
effects such as secondary emission and collisions. A detailed introduction to
SLPIC is provided in §II.2. This work is the first demonstration of the
integration of SLPIC particles with PIC-MCC, showing that speed-limiting
techniques can be used in collisional low-temperature plasma discharges. We
chose to simulate electrical breakdown, and, specifically, Paschen’s law,
because multiple collision mechanisms must be correctly simulated to compute
the breakdown voltage accurately. We were also able to quantify the speed-up
in runtime of the SLPIC simulations relative to the regular PIC simulations.

In this paper, we begin with a review of Paschen curves and SLPIC. We then
introduce the SLPIC Monte Carlo Collisions framework in §III. The simulation
design and parameters will be provided along with our methodology for
determining breakdown in §IV. In §V, we compare our simulation results (a
Paschen curve) with experimental data and quantify the computational speed-up
of the SLPIC algorithm relative to normal PIC.

##  II Background

###  II.1 Paschen Curve

Electrical breakdown in a gas-filled capacitor is caused by the ionization of
the neutral gas by electrons and the secondary emission of electrons at the
cathode due to ion impact. The average number of ions produced per electron
per unit length travelled is termedTownsend (1900) Townsend’s first ionization
coefficient $\alpha$. The average number of secondary electrons emitted per
ion impact is Townsend’s second ionization coefficient $\gamma_{se}$. The
coefficient $\alpha$ depends on the voltage $V$ across the gap and the
collision cross-sections. The coefficient $\gamma_{se}$ depends on the ion
species and the electrode material and treatment. When $\gamma_{se}(e^{\alpha
d}-1)>1$ breakdown occurs.

Paschen curves relate the breakdown voltage $V_{b}$ of a gas-filled capacitor
to the product of the gas pressure $p$ and the distance $d$ between the
electrodes.Paschen (1889) Paschen curves are not monotonic due to a
competition between electrons gaining energy in the electric field and losing
energy in collisions. Paschen curves have a characteristic minimum at an
intermediate $pd$. At low $pd$, most electrons stream through the capacitor
without ionizing enough atoms. At high $pd$, few electrons reach the threshold
ionization energy due to frequent non-ionizing collisions. Paschen’s law is a
function relating $V_{b}$ to $pd$:

$V_{b}=\frac{Bpd}{\ln(Apd)-\ln[\ln(1+1/\gamma_{se})]}.$ (1)

Paschen’s law assumes that, except for the initial seed electrons, ion impact
at the cathode and electron-neutral ionizations are the only sources of
electrons, and that $\alpha$ can be expressed as:

$\frac{\alpha}{n_{g}}=kTA\exp\left(\frac{-kTBn_{g}}{E}\right)$ (2)

where $kT$ is the thermal energy, $n_{g}$ is the neutral gas density, and
$E=V/d$ is the constant electric field. $A$ and $B$ are related to the
collision cross-sections and are determined by fitting Eq. (2) (which gives
the field-intensified ionization cross-section $\alpha/n_{g}$ as a function of
the reduced field $E/n_{g}$) to the experimental results for a given gas. In
Paschen’s law, $A$ and $B$ are assumed to be constant for a given gas species.
Paschen’s law will be plotted in Fig. 2, with the values of $A$, $B$, and
$\gamma_{se}$ taken from Lieberman and Lichtenberg (2005).

Different gases yield different Paschen curves due to different cross-
sections, ionization energies, and secondary emission yields. In this work, we
simulated a neutral argon gas. Argon has been studied extensively in the
context of electrical breakdown. We will compare our simulation results to
experimental results compiled by Phelps and Petrovic (1999).

###  II.2 SLPIC

Before describing how collisions can be implemented within the SLPIC method in
§III, we review the derivation of SLPIC, following Werner  _et al._ (2018).
PIC numerically evolves a particle distribution $f(\mathbf{x},\mathbf{v},t)$
according to the Vlasov equation:

$\displaystyle\partial_{t}f(\mathbf{x},\mathbf{v},t)$ $\displaystyle=$
$\displaystyle-\mathbf{v}\cdot\partial_{\bf
x}f(\mathbf{x},\mathbf{v},t)-\mathbf{a}\cdot\partial_{\bf
v}f(\mathbf{x},\mathbf{v},t),$ (3)

where the acceleration $\mathbf{a}=\mathbf{a}(\mathbf{x},\mathbf{v},t)$ is
typically determined by electromagnetic fields. SLPIC numerically evolves
$f(\mathbf{x},\mathbf{v},t)$ according to an _approximate_ Vlasov equation. We
multiply Eq. (3) by a yet-to-be-chosen function $\beta(v)$ (where
$v=\|\mathbf{v}\|$), which will turn out to limit particle speeds,

$\displaystyle\beta\partial_{t}f(\mathbf{x},\mathbf{v},t)$ $\displaystyle=$
$\displaystyle-\beta\mathbf{v}\cdot\partial_{\bf
x}f(\mathbf{x},\mathbf{v},t)-\beta\mathbf{a}\cdot\partial_{\bf
v}f(\mathbf{x},\mathbf{v},t),$ (4)

and then make the critical “speed-limiting” approximation:

$\displaystyle\beta\partial_{t}f$ $\displaystyle\approx$
$\displaystyle\partial_{t}f.$ (5)

This approximation is valid if either (1)
$\partial_{t}f(\mathbf{x},\mathbf{v},t)\approx 0$ or (2) $\beta(v)\approx
1$.111This approximation is subtly different from that in Werner  _et al._
(2018), but only when $\beta$ has an explicit time-dependence—a regime that
has not yet been investigated and is irrelevant for this paper. Thus we arrive
at the “speed-limited” approximation to the Vlasov equation

$\displaystyle\partial_{t}f(\mathbf{x},\mathbf{v},t)$ $\displaystyle=$
$\displaystyle-\beta\mathbf{v}\cdot\partial_{\bf
x}f(\mathbf{x},\mathbf{v},t)-\beta\mathbf{a}\cdot\partial_{\bf
v}f(\mathbf{x},\mathbf{v},t).$ (6)

In the special case of a steady state (i.e., $\partial_{t}f\equiv 0$), the
speed-limiting approximation is exact and SLPIC, though different from PIC, is
as accurate as PIC. For the special choice of $\beta(v)\equiv 1$, SLPIC is
exactly the same as PIC.

Just as PIC simulation evolves the Vlasov equation via the method of
characteristics, SLPIC evolves Eq. (6). We posit a solution that is a sum over
macroparticles $p$ that follow trajectories
$[\mathbf{x}_{p}(t),\mathbf{v}_{p}(t)]$ with weights $w_{p}(t)$ (a macro-
electron with weight $w_{p}$ represents $w_{p}$ electrons; i.e., it has mass
$w_{p}m_{e}$ and charge $-w_{p}e$):

$\displaystyle f(\mathbf{x},\mathbf{v},t)$ $\displaystyle=$
$\displaystyle\sum_{p}w_{p}(t)S[\mathbf{x}-\mathbf{x}_{p}(t)]\delta^{3}[\mathbf{v}-\mathbf{v}_{p}(t)]$
(7)

where $S$ is the particle shape functionBirdsall (1991) and $\delta$ the Dirac
delta function. Unlike PIC, SLPIC requires the particle weight to vary in
time. Plugging Eq. (7) into Eq. (6), we find that it yields a solution if:

$\displaystyle\dot{\mathbf{x}}_{p}$ $\displaystyle=$
$\displaystyle\beta(v_{p})\mathbf{v}_{p}$ (8)
$\displaystyle\dot{\mathbf{v}}_{p}$ $\displaystyle=$
$\displaystyle\beta(v_{p})\mathbf{a}(\mathbf{x}_{p},\mathbf{v}_{p},t)$ (9)
$\displaystyle\dot{w}_{p}$ $\displaystyle=$ $\displaystyle w_{p}[\partial_{\bf
x}\cdot(\beta\mathbf{v})+\partial_{\bf v}\cdot(\beta\mathbf{a})]=\beta
w_{p}[\partial_{\bf x}\cdot\mathbf{v}+\partial_{\bf
v}\cdot\mathbf{a}]+w_{p}\mathbf{a}\cdot\partial_{\bf
v}\beta=w_{p}\dot{\beta}/\beta$ (10)

where we have assumed that $\partial_{\bf x}\cdot\mathbf{v}+\partial_{\bf
v}\cdot\mathbf{a}=0$—i.e., the true motion must be phase-space-volume
preserving (Liouville’s theorem proves this for any Hamiltonian motion). We
have also substituted
$\dot{\beta}\equiv(d/dt)\beta(v_{p}(t))=\beta\mathbf{a}\cdot\partial_{\bf
v}\beta$.

When $\beta(v_{p})=1$, these are the usual (PIC) equations of motion. When
$0<\beta(v_{p})<1$, a SLPIC macroparticle, representing physical particles
with true velocity $\mathbf{v}_{p}$, travels with pseudo-velocity $\dot{\bf
x}_{p}=\beta(v_{p})\mathbf{v}_{p}$ and pseudo-acceleration $\dot{\bf
v}_{p}=\beta(v_{p})\mathbf{a}(\mathbf{x}_{p}(t),\mathbf{v}_{p}(t),t)$. The
SLPIC macroparticle follows the same path through phase space as the physical
particles it represents, but slower; e.g., in a short time $\Delta t$, the
macroparticle moves $({\bf v}_{p},{\bf a})\beta\Delta t$, while a physical
particle moves $({\bf v}_{p},{\bf a})\Delta t$. For example, a SLPIC electron
in a uniform magnetic field will have the correct gyroradius but a
gyrofrequency that is a factor $\beta(v_{p})$ too small.

Choosing $\beta(v_{p})<1$ invokes the SLPIC approximation [Eq. (5)], reducing
accuracy; however, it limits the simulated speed of macroparticles from
$v_{p}$ to $\dot{x}=\beta(v_{p})v_{p}$, which allows larger timesteps, hence
faster simulation. In most cases, the PIC timestep must be smaller than
$\Delta t\lesssim\Delta x/v_{\rm max}$, where $\Delta x$ is the grid cell size
and $v_{\rm max}$ is the maximum particle speed. With SLPIC, a larger timestep
is allowed: $\Delta t\lesssim\Delta x/[\beta(v_{\rm max})v_{\rm max}]$.
Typically, we choose a speed limit $v_{0}\ll v_{\rm max}$, and then define
$\beta(v_{p})\equiv v_{0}/v_{p}$ to ensure that for all particles, $\|\dot{\bf
x}\|\leq v_{0}$, so that we can use $\Delta t\approx\Delta x/v_{0}\gg\Delta
x/v_{\rm max}$. For particles with $v_{p}<v_{0}$, however, we define
$\beta(v_{p})\equiv 1$ so that they experience “true” motion. In this way, we
speed-limit only the particles (with $v_{p}>v_{0}$) that need to be speed-
limited to allow a large timestep, $\Delta t=v_{0}/\Delta x$. It turns out
that speed-limiting also reduces the plasma frequency so that it is resolved
by the large $\Delta t$ and poses no instability.Werner  _et al._ (2018);
Jenkins, Werner, and Cary (2021) For simplicity, we will use this specific
$\beta(v_{p})$ for the rest of the paper: $\beta(v_{p})=v_{0}/v_{p}$ for
$v_{p}\geq v_{0}$, and $\beta(v_{p})=1$ for $v_{p}<v_{0}$, or, equivalently:

$\displaystyle\beta(v_{p})$ $\displaystyle\equiv$
$\displaystyle\Theta(v_{p}-v_{0})v_{0}/v_{p}+\Theta(v_{0}-v_{p})$ (11)

where $\Theta$ is the Heaviside step function (and $v_{p}=\|\mathbf{v}\|$).

As a macroparticle accelerates so that $v_{p}>v_{0}$, its simulated speed
$\dot{x}_{p}=\beta v_{p}$ remains at $v_{0}$, and its weight $w_{p}(t)$
decreases. The decrease in $w_{p}(t)$ reflects the fact that the physical
particles represented by the macroparticle at time $t$ outrun it, and the
macroparticle represents a different set of fewer physical particles at time
$t+\Delta t$. In this way, SLPIC can correctly simulate a steady-state beam
within which particles accelerate such that their density decreases (like cars
exiting a steady-state traffic jam). In SLPIC, a beam macroparticle with
$v_{p}>v_{0}$ experiences an acceleration and its true velocity $v_{p}$
increases, but its pseudo-velocity $\dot{x}_{p}=\beta v_{p}=v_{0}$ does not,
so the density of macroparticles throughout the beam remains uniform. The
decrease in physical density thus results from the decrease in macroparticle
weight, and not (as in PIC) from the decrease in macroparticle density.

Therefore, as $v_{p}$ increases above $v_{0}$, $\beta(v_{p})$ decreases, and
so too must $w_{p}$ decrease. Equation (10) shows that $w_{p}/\beta(v_{p})$ is
constant over any particle trajectory; we define this constant to be
$w_{j,p}\equiv w_{p}/\beta(v_{p})$. Here we choose the subscript $j$ because
$w_{j,p}$ facilitates calculation of the flux density distribution ${\bf
j}(\mathbf{x},\mathbf{v},t)$.

Simulating macroparticles following Eqs. (8–10), we can compute the phase-
space density distribution function $f$ at any time via Eq. (7). However, with
SLPIC it is often especially useful to consider particle fluxes. We can
compute the flux density distribution,
$\mathbf{j}(\mathbf{x},\mathbf{v},t)=\mathbf{v}f(\mathbf{x},\mathbf{v},t)$, as
follows:

$\displaystyle\mathbf{j}(\mathbf{x},\mathbf{v},t)=\sum_{p}\mathbf{v}_{p}w_{p}(t)S[\mathbf{x}-\mathbf{x}_{p}(t)]\delta^{3}[\mathbf{v}-\mathbf{v}_{p}(t)]=\sum_{p}\beta(v_{p})\mathbf{v}_{p}w_{j,p}S[\mathbf{x}-\mathbf{x}_{p}(t)]\delta^{3}[\mathbf{v}-\mathbf{v}_{p}(t)].$

To estimate the flux $\mathbf{j}\cdot\hat{\bf n}dA$ through some surface
element $\hat{\bf n}dA$ from the macroparticles crossing it in time $\Delta
t$, we consider that a macroparticle with $\dot{\bf x}$ will cross the surface
if it is within a distance $\dot{\bf x}\Delta t\cdot\hat{\bf n}$ of the
surface. From Eq. (II.2) we see that, since $\dot{\bf
x}=\beta(v_{p})\mathbf{v}_{p}$, the flux (integrated over time $\Delta t$)
through the surface element is simply the sum of the $w_{j,p}$ over all
macroparticles crossing the surface in $\Delta t$.

Because the $w_{j,p}$ are constant in time, a macroparticle always represents
the same physical flux whenever it crosses a surface. If a macroparticle $p$
has $w_{j,p}=100$, then every time it crosses a surface, it represents 100
physical particles crossing the surface. This statement is trivial in PIC, but
in SLPIC, where $w_{p}(t)$ varies in time, it is nontrivial and very useful.
For example, it ensures that, if a macroparticle enters into some volume $V$
and later exits that volume, the time-integrated flux through the surface (due
to that particle) is zero.

At this point readers might suspect a paradox in the SLPIC treatment arising
from a difference between $w_{p}$ and $w_{j,p}$. Consider a macroparticle with
$v_{p}=10v_{0}$, $\beta(v_{p})=0.1$, $w_{j,p}=100$, and $w_{p}=10$, and
suppose it crosses a surface to enter volume $V$ during time interval $\Delta
t$. That macroparticle represents $w_{p}=10$ physical particles in volume $V$;
however, it represented $w_{j,p}=100$ physical particles crossing the surface
into $V$. This counterintuitive behavior is correct for SLPIC. It is
counterintuitive because PIC has led us to assume that a macroparticle with
$w_{p}$ and $(\mathbf{x}_{p},\mathbf{v}_{p})$ represents a swarm of $w_{p}$
physical particles within the volume $d^{3}\\!x\,d^{3}\\!v$ around
$(\mathbf{x}_{p},\mathbf{v}_{p})$, and that those physical particles travel
roughly with the macroparticle to $(\mathbf{x}_{p}+\dot{\bf x}_{p}\Delta
t,\mathbf{v}_{p}+\dot{\bf v}_{p}\Delta t)$ over time $\Delta t$.
Fundamentally, however, a macroparticle represents a chunk of the distribution
$fd^{3}\\!x\,d^{3}\\!v$, and SLPIC decouples macroparticles from the physical
particles they represent. A speed-limited SLPIC macroparticle moves in phase
space with $\dot{\bf x}_{p}$ and $\dot{\bf v}_{p}$ slower than the physical
particles it represents. Thus the same macroparticle may represent one set of
physical particles at time $t$ and a different set of particles at time
$t+\Delta t$.

In the above example, the density represented by the SLPIC macroparticle
corresponds to $w_{p}=10$ physical particles that are near
$[\mathbf{x}_{p}(t),\mathbf{v}_{p}(t)]$ at time $t$; however, over time
interval $\Delta t$, the physical flux represented by the macroparticle
includes all the physical particles that would would be near
$[\mathbf{x}_{p}(t^{\prime}),\mathbf{v}_{p}(t^{\prime})]$ at any time
$t^{\prime}\in[t,t+\Delta t]$. Since the physical particles travel 10 times
faster than the macroparticle, the physical flux is $10w_{p}(t)/\Delta
t=w_{j,p}/\Delta t$.

It is important to remember that SLPIC is accurate only if
$f(\mathbf{x},\mathbf{v},t)$ changes sufficiently slowly—and it is as accurate
as PIC in the steady-state limit. The $f(\mathbf{x},\mathbf{v},t)$ describing
a single particle with $v>v_{0}$ is not a slowly-changing function; therefore,
one typically cannot verify SLPIC based on single-particle thought
experiments. In the above single-particle example, the continuity equation is
violated because $f(\mathbf{x},\mathbf{v},t)$ changes too rapidly for the
SLPIC approximation to be valid.

Reiterating an important point: the flux weight
$w_{j,p}=w_{p}(t)/\beta(v_{p}(t))$ of a SLPIC macroparticle is a constant in
time, and this ensures that the flux density is divergenceless in the steady-
state limit. For example, if a particle receives an abrupt kick—whether from
the electric field or from a collision—its $w_{j,p}$ remains constant. If the
impulse increases the particle’s speed $v_{p}$, its density weight $w_{p}$
must change accordingly so that $w_{p}=\beta(v_{p})w_{j,p}$. In the next
section we will see that the flux weight is conserved when simulating
collisions involving SLPIC particles.

##  III Collisions in SLPIC

Collisions can be simulated in PIC with a PIC-MCC algorithm;Birdsall (1991);
Bird (1994); Nanbu (2000) the MCC algorithm can be used for SLPIC with just a
few modifications: (1) using the flux weight $w_{j,p}$ instead of the density
weight $w_{p}$ (cf. §II.2), and (2) considering the slowing of time when
determining the collision rate. (Since most SLPIC simulations will have
macroparticles with a distribution of $w_{j,p}$, the PIC-MCC algorithm must be
able to handle macroparticles with different weights before modification for
SLPIC.)

Before discussing collisions in SLPIC, we briefly review the pertinent aspects
of collisions in PIC. In PIC-MCC, a binary macroparticle collision is not at
all the same as a binary collision between two physical particles. The
collision between two macroparticles must statistically represent collisions
between the two swarms of physical particles represented by the
macroparticles. The spatial distribution of physical particles is usually
assumed to be uniform within the same grid cell containing the macroparticle
centers, regardless of the macroparticle shape $S$ [cf. Eq. (7)], and zero
outside that cell.

For concreteness, let us consider electron impact ionization in PIC-MCC. A
primary electron macroparticle $p$ with weight $w_{p0}$ and velocity
$\mathbf{v}_{p0}$ collides with a neutral macro-atom with $w_{n0}$ and
$\mathbf{v}_{n0}$. Each primary macroparticle will be split into two parts:
one part (weights $w_{p1}$ and $w_{n1}$) represents uncollided particles with
unchanged velocities, while the other part ($w_{p1}^{\prime}$ and
$w_{n1}^{\prime}$) represents the scattered particles—in principle there can
be many scattered macroparticles with different velocities, but for simplicity
we will assume at most one per primary species, with velocities
$\mathbf{v}_{p1}^{\prime}$ and $\mathbf{v}_{n1}^{\prime}$. In addition, a
secondary electron and ion are created with $w_{se}$, $\mathbf{v}_{se}$, and
$w_{si}$, $\mathbf{v}_{si}$. Conservation of particle number requires
$w_{p0}=w_{p1}+w_{p1}^{\prime}$, $w_{n0}=w_{n1}+w_{n1}^{\prime}$, and
$w_{p1}^{\prime}=w_{n1}^{\prime}=w_{se}=w_{si}=W_{c}$, where $W_{c}$ is the
number of (physical-particle) collisions expected to take place when the two
macroparticles “collide.”

The expected number of collisions $W_{c}$ within a time interval $\Delta t$ in
a cell volume $\Delta V$ between two arbitrary particle distributions, $f_{p}$
(primary electrons) and $f_{n}$ (neutrals) is

$\displaystyle W_{c}$ $\displaystyle=$ $\displaystyle\Delta V\Delta t\int
d^{3}\mathbf{v}_{p0}\int
d^{3}\mathbf{v}_{n0}\,f_{p}(\mathbf{v}_{p0})f_{n}(\mathbf{v}_{p0})\sigma(\mathbf{v}_{n0}-\mathbf{v}_{p0})\|\mathbf{v}_{n0}-\mathbf{v}_{p0}\|$
(13)

where $\sigma(\mathbf{v}_{r})$ is the collision cross-section, which depends
on the relative velocity $\mathbf{v}_{r}$ between colliding particles. To
calculate $W_{c}$ between two PIC macroparticles in the same grid cell, we
substitute the appropriate distributions, i.e.,
$f_{p}(\mathbf{v})=(w_{p0}/\Delta V)\delta^{3}(\mathbf{v}-\mathbf{v}_{p0})$
and $f_{n}(\mathbf{v})=(w_{n0}/\Delta
V)\delta^{3}(\mathbf{v}-\mathbf{v}_{n0})$, ignoring the shape functions of the
macroparticles, treating them as uniform in one cell and zero outside. The
expected number of collisions between the two macroparticles is then

$\displaystyle W_{c}$ $\displaystyle=$ $\displaystyle\Delta
V\frac{w_{p0}}{\Delta V}\frac{w_{n0}}{\Delta
V}\sigma(\mathbf{v}_{n0}-\mathbf{v}_{p0})\|\mathbf{v}_{n0}-\mathbf{v}_{p0}\|\Delta
t.$ (14)

This can be verified by considering a single electron moving through a swarm
of atoms of density $n_{n0}=w_{n0}/\Delta V$ with velocity $\mathbf{v}_{n0}$.
In the reference frame co-moving with the atoms,
$d=\|\mathbf{v}_{n0}-\mathbf{v}_{p0}\|\Delta t$ is the distance traveled by
the electron in $\Delta t$, and $\lambda_{\rm mfp}=(n_{n0}\sigma)^{-1}$ is its
mean free path. Therefore $d/\lambda_{\rm mfp}$ is the expected number of
collisions that a single electron would experience (for small $\Delta t$ such
that $d/\lambda_{\rm mfp}\ll 1$). The total number of electron collisions is
then $W_{c}=w_{p0}d/\lambda_{\rm mfp}$, which is identical to Eq. (14).

There are multiple Monte Carlo strategies for colliding pairs of
macroparticles that all yield the same $W_{c}$ _on average_. For example, if
$W_{c}=0.01w_{p0}$, we could avoid splitting the primary macroparticles and
creating small-weight scattered/secondary particles by ignoring 99% of
colliding pairs, and colliding 1% with an enhanced $W_{c}^{\prime}=w_{p0}$. We
will not discuss this in any detail, because these strategies are identical
for PIC and SLPIC. We also omit other practical details that need no
modification for SLPIC, such as determining scattered and secondary
macroparticle velocities according to given differential cross-sections, and
dealing with exceptional cases such as $W_{c}>w_{p0}$ (indicating that $\Delta
t$ is too large). The only differences introduced by SLPIC are the
determination of $W_{c}$ and the resulting macroparticle weights.

As described in §II.2, in SLPIC it is better to verify physical accuracy by
considering particle fluxes rather than numbers or densities. Thus the basic
thought-experiment for developing SLPIC collision algorithms is not the
collision of two macroparticles, but of two narrow beams of macroparticles
intersecting in a small volume; the particle fluxes into and out of that
volume must satisfy the appropriate conservation laws. For electron impact
ionization, this means that $w_{j,p0}=w_{j,p1}+w_{j,p1}^{\prime}$,
$w_{j,n0}=w_{j,n1}+w_{j,n1}^{\prime}$, and
$w_{j,p1}^{\prime}=w_{j,n1}^{\prime}=w_{j,se}=w_{j,si}=W_{j,c}$. Because flux
weights $w_{j}$ are constant in time [unlike density weights $w_{p}(t)$], this
ensures that the fluxes of particles entering and exiting $\Delta V$ are
appropriately perserved; i.e., the flux of unscattered plus scattered
primaries leaving $\Delta V$ equals the flux of primaries entering $\Delta V$,
and the fluxes of secondary electrons and scattered primaries leaving $\Delta
V$ are equal, etc.

The rule for determining the collision rate, $W_{j,c}/\Delta t$, is that SLPIC
preserves the physical mean free path. For example, consider a steady-state,
monoenergetic beam of particles in the $+x$ direction, entering some
scattering medium with mean free path $\lambda_{\rm mfp}$. Scattering causes
the density of original beam particles to decrease as
$\sim\exp(-x/\lambda_{\rm mfp})$. Because this is a steady-state, SLPIC should
render this density profile accurately, requiring that SLPIC macroparticles
experience the physical mean free path. In other words, SLPIC macroparticle
must scatter with probability $d^{\prime}/\lambda_{\rm mfp}=\dot{x}\Delta
t/\lambda_{\rm mfp}=\beta(v)v\Delta t/\lambda_{\rm mfp}$ within time $\Delta
t$ (in contrast, a physical particle would scatter with probability
$d/\lambda_{\rm mfp}=\dot{x}\Delta t/\lambda_{\rm mfp}=v\Delta t/\lambda_{\rm
mfp}$). This can also be viewed as the reduction in collision rate due to the
slowing down of time by a factor $\beta(v)$. The collision rate (in terms of
fluxes) in $\Delta V$ must therefore be

$\displaystyle\frac{W_{j,c}}{\Delta t}$ $\displaystyle=$
$\displaystyle\frac{w_{j,p0}}{\Delta t}\frac{d^{\prime}}{\lambda_{\rm
mfp}}=\Delta V\frac{w_{j,p0}/\Delta t}{\Delta
V}n_{n0}\sigma(\mathbf{v}_{n0}-\mathbf{v}_{p0})\|\mathbf{v}_{ni}-\mathbf{v}_{pi}\|\beta_{p}(v_{p0})\Delta
t$ (15)

where $\beta_{p}(v)$ is the speed-limiting function for electrons. This merely
says that—if one were to measure fluxes into and out of $\Delta V$—the flux of
scattered electrons (i.e., $W_{j,c}/\Delta t$) would equal the flux of
incident electrons (i.e., $w_{j,p0}/\Delta t$) times $d^{\prime}/\lambda_{\rm
mfp}$.

The above expression involves the neutral density: $n_{n0}=w_{n0}/\Delta
V=\beta_{n}(v_{n0})w_{j,n0}/\Delta V$. (In most applications, atoms will be
slow and not speed-limited, hence $\beta_{n}\equiv 1$, but we want this
treatment to extend to arbitrary collisions.) Thus:

$\displaystyle W_{j,c}$ $\displaystyle=$ $\displaystyle\Delta
V\frac{w_{j,p0}}{\Delta V}\frac{\beta_{n}(v_{n0})w_{j,n0}}{\Delta
V}\sigma(\mathbf{v}_{n0}-\mathbf{v}_{p0})\|\mathbf{v}_{n0}-\mathbf{v}_{p0}\|\beta_{p}(v_{p0})\Delta
t$ (16) $\displaystyle=$ $\displaystyle\Delta V\frac{w_{j,p0}}{\Delta
V}\frac{w_{j,n0}}{\Delta
V}\sigma(\mathbf{v}_{n0}-\mathbf{v}_{p0})\|\mathbf{v}_{n0}-\mathbf{v}_{p0}\|\beta_{n}(v_{n0})\beta_{p}(v_{p0})\Delta
t.$

This expression (or rather $W_{j,c}/\Delta t$) yields the number of physical
electrons scattered per $\Delta t$ from the collision of particles represented
by two SLPIC macroparticles. It is the same as Eq. (14), except: (1) flux
weights $w_{j}$ are used (not density weights $w$), and (2) the timestep is
modified by a factor $\beta_{n}(v_{n0})\beta_{p}(v_{p0})$.

After calculating $W_{j,c}$ using Eq. (16), the rest of the collision
algorithm proceeds as in PIC, except that where a weight was used in the PIC
algorithm, the flux weight must be used in SLPIC. Once a macroparticle’s flux
weight $w_{j,p}$ is determined, its density weight is set by
$w_{p}=\beta(v)w_{j,p}$.

##  IV Simulation Setup and Methodology

We modelled an argon-filled parallel plate capacitor with one spatial
dimension and three velocity dimensions. The electrodes were treated as
particle-absorbing boundaries and the gap-distance was fixed at $d=1$ cm. We
imposed a constant electric field $E=V/d$. The cell size $\Delta x$ was set to
a quarter of the mean free path. The number of cells for each simulation is
given in Table 1. We ran both SLPIC and PIC simulations. The setup parameters
for the SLPIC and PIC simulations were identical except for the timestep. For
PIC, we used a timestep given by $\Delta t=\Delta x/v_{\rm max}$, where
$v_{\rm max}=\sqrt{2eV/m_{e}}$ was the maximum electron speed. For SLPIC, we
used a timestep given by $\Delta t=\Delta x/v_{0}$, where $v_{0}=v_{\rm
max}\sqrt{m_{e}/m_{Ar}}$ was the electron speed limit imposed by
$\beta(v_{p})$ given in Eq. (11). The large mass difference between electrons
and argon ions resulted in the SLPIC timestep being 270 times larger than the
PIC timestep. Therefore, the SLPIC simulations required 270 times fewer
timesteps than PIC. The number of timesteps for each simulation is given in
Table 1.

To initiate the discharge we injected 100 electrons at the cathode in the
first timestep with zero initial velocity. We simulated five types of
electron-neutral collisions: ionization, elastic collisions, and three
excitations. Many different types of collisions occur between electrons and
argon, but we simulated only the five most likely. The other collisions had
negligible cross-sections in the relevant electron energy range. The cross-
sections were extracted from the Biagi-v7.1 databaseBiagi (1999) and Phelps
databaseYamabe, Buckman, and Phelps (1983) on www.lxcat.net on August 15,
2020. We did not include any collisions for the ions. For excitation
collisions, the electrons were scattered isotropically with an energy reduced
by the excitation threshold. For elastic collisions, the electrons were
scattered according to the Vahedi-Surendra algorithm.Vahedi and Surendra
(1995) For ionization collisions, the products were generated according to an
algorithm developed by Kutasi and Donkó (2000). We used an energy-dependent
model of secondary electron emission due to ion impact at the cathode.Phelps
and Petrovic (1999) Secondary electrons were emitted with zero velocity with a
probability given by Eq. (17) that depends weakly on the energy of the
incident ion $\epsilon_{i}$.

$\gamma_{se}=\begin{cases}0.09\left(\frac{\epsilon_{i}}{700\,\rm{eV}}\right)^{0.05},&\epsilon_{i}<700~{}\rm{eV}\\\
0.09\left(\frac{\epsilon_{i}}{700\,\rm{eV}}\right)^{0.72},&\epsilon_{i}\geq
700~{}\rm{eV}\\\ \end{cases}$ (17)

To determine whether the simulation voltage $V$ was above or below $V_{b}$, we
tracked the ion population in the simulation over 30 ion crossing times
($30\sqrt{2m_{Ar}d/eE}$). For $V<V_{b}$, the ion population decreased after an
initial rise, eventually approaching to zero. For $V>V_{b}$, the ion
population increased exponentially. This method provided a clear binary
classification method of simulation results. For each pressure, we performed a
grid search in the vicinity of the experimental $V_{b}$. Once a simulation
above and a simulation below breakdown were found, we bracketed the interval
with the respective voltages. These brackets are indicated by the bars in Fig.
2. The evolutions of the ion population for simulations above and below
breakdown at 1 Torr cm are shown in Fig. 1, for both SLPIC and PIC.

Figure 1: The number of ions in the simulation domain as a function of time
for the simulations with $V$ nearest $V_{b}$ at $pd=1$ Torr cm. The SLPIC
results are given for $V>V_{b}$ (orange) and $V<V_{b}$ (cyan). The PIC results
are given for $V>V_{b}$ (black) and $V<V_{b}$ (purple). The number of ions in
the simulation domain increases over time when $V>V_{b}$, and decreases when
$V<V_{b}$.

We bracketed the breakdown voltage for 11 values of $pd$ ranging from 0.3 to
300 Torr cm. We performed simulations with the VorpalNieter and Cary (2004)
code distributed in VSim-11. The simulations ran on the CORI supercomputer at
the National Energy Research Scientific Computing Center (NERSC). The number
of cores for each simulation is given in Table 1.

##  V Simulation Results

Table 1: Summary of setup parameters and performance for PIC and SLPIC simulations with $V\approx V_{b}$. The gap-distance was fixed at $d=$1 cm. The number of timesteps is given by $N_{T}$.“Cores” refers to the number of CORI cores that the simulation was run on. The runtime $T$ is given in core-hours. The speed-up is given by $T_{\rm PIC}/T_{\rm SLPIC}$. $p$ (Torr) |  $V$ (V) | Cells | $N_{T,\rm SLPIC}$ | $N_{T,\rm PIC}$ | Cores | $T_{\rm SLPIC}$ | $T_{\rm PIC}$ | Speed-up
---|---|---|---|---|---|---|---|---
0.3 | 260 | 66 | 4.0$\times 10^{3}$ | 1.1$\times 10^{6}$ | 4 | 0.26 | 11 | 42
0.4 | 195 | 88 | 5.3$\times 10^{3}$ | 1.4$\times 10^{6}$ | 4 | 0.27 | 14 | 52
0.6 | 175 | 131 | 7.9$\times 10^{3}$ | 2.1$\times 10^{6}$ | 4 | 0.49 | 91 | 187
1 | 170 | 219 | 1.3$\times 10^{4}$ | 3.5$\times 10^{6}$ | 4 | 0.78 | 113 | 146
2 | 185 | 437 | 2.6$\times 10^{4}$ | 7.1$\times 10^{6}$ | 4 | 1.8 | 91 | 51
4 | 225 | 874 | 5.2$\times 10^{4}$ | 1.4$\times 10^{7}$ | 4 | 6.3 | 1308 222time based on extrapolation due to time restrictions on CORI | 208
6 | 262 | 1310 | 7.9$\times 10^{4}$ | N/A | 4 | 9.2 | N/A | N/A
10 | 340 | 2184 | 1.3$\times 10^{5}$ | N/A | 4 | 13 | N/A | N/A
30 | 620 | 6550 | 3.9$\times 10^{5}$ | N/A | 128 | 105 | N/A | N/A
100 | 1400 | 21832 | 1.3$\times 10^{6}$ | N/A | 128 | 771 | N/A | N/A
300 | 3100 | 65494 | 3.9$\times 10^{6}$ | N/A | 128 | 6400 a | N/A | N/A
Figure 2: The Paschen curve for argon. The experimental data (solid, green
circles) was extracted from Phelps and Petrovic (1999). The values of $A$,
$B$, and $\gamma_{se}$ in Paschen’s law (solid, black line) were taken from
Lieberman and Lichtenberg (2005). The SLPIC results (blue bars connected by
dashed, blue line) display the upper and lower brackets on the simulated
breakdown voltage. The 5 lowest pressures (open, red circles) were also
simulated using PIC, which yielded the same bounds as SLPIC.

SLPIC accurately computed the Paschen curve for argon over three orders of
magnitude in $pd$, 0.3 to 300 Torr cm, agreeing with PIC over the range where
PIC was feasible, 0.3 to 2 Torr cm. For each value of $pd$, we bracketed the
breakdown voltage with one simulation above and below breakdown. The Paschen
curve generated by the SLPIC and PIC simulations is shown along with Paschen’s
law and experimental data from Phelps and Petrovic (1999) in Fig. 2. For the
values of $pd$ where PIC was run, PIC and SLPIC classified each simulation
identically, thus yielding the same brackets. Below 3 Torr cm, the simulated
breakdown voltages were lower than those measured in experiment. This may have
been due to our idealized treatment of the electrodes. The electrodes in the
experiments were exposed to gases that may have caused oxidation or other
contaminating reactions,Phelps and Petrovic (1999) which decreases
$\gamma_{se}$ for ions with energies below 200 eV and therefore increases
$V_{b}$ in the low voltage regime. Paschen’s law, given in Eq. (1), is plotted
using the values of $A=11.5$ cm-1 Torr-1, $B=176$ V cm-1 Torr-1, and
$\gamma_{se}=0.07$ taken from Lieberman and Lichtenberg (2005). Below 0.3 Torr
cm, the simulation did not break down. This matches Paschen’s law [Eq. (1)]
which exhibits a singularity below 0.3 Torr cm. The experimental data,
however, does show breakdown at low pressures because at low pressures the
dominant mechanism for secondary emission is photoionization,Phelps and
Petrovic (1999) which was not simulated. At high pressures, the simulated
breakdown voltages agree well with the experimental data.

SLPIC ran 40 to 200 times faster than PIC. The simulation parameters and
performance for a subset of simulations run near the breakdown voltage are
given in Table 1. The speed-ups are the ratios of the PIC runtimes to the
SLPIC runtimes and are given in the rightmost column. On average, the SLPIC
simulations in Table 1 ran 116 times faster than PIC. Above 4 Torr cm, it
became unfeasible to run PIC simulations. SLPIC enabled us to explore these
higher pressures with a fully kinetic treatment.

##  VI Summary

Fast, accurate electrical discharge simulations are needed for the design of
plasma processing equipment, sensitive microchips, and volatile chemical
processing facilities. We have demonstrated that SLPIC is as accurate as PIC,
but faster in predicting breakdown voltages of a gas-filled capacitor. SLPIC
accurately computed the Paschen curve $V_{b}(pd)$ for argon for $pd$ ranging
from 0.3 to 300 Torr cm. SLPIC and PIC produced identical results, but SLPIC
ran 40 to 200 times faster and extended the range of feasible simulations. In
accurately computing the Paschen curve for argon, SLPIC has demonstrated that
it can accurately model collisions, including electron-neutral ionization,
electron-neutral elastic collisions, electron-neutral excitations, and ion-
induced secondary-electron emission. We expect that SLPIC will also be as
accurate as PIC, but much faster, in simulating voltage breakdown in gases in
more complicated 3D geometries as well as in the glow discharge regime where
the discharge current alters the applied electric field.

###### Acknowledgements.

This work was supported by NSF grant PHY1707430 and by U.S. Department of
Energy SBIR Phase II Award DE-SC0015762. This research used resources of the
National Energy Research Scientific Computing Center (NERSC), a U.S.
Department of Energy Office of Science User Facility operated under Contract
No. DE-AC02-05CH11231.

## AIP Publishing Data Sharing Policy

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Townsend (1900) J. Townsend, Nature 62, 340 (1900).
  * Paschen (1889) F. Paschen,  _Ueber die zum Funkenübergang in Luft: Wasserstoff und Kohlensäure bei verschiedenen Drucken erforderliche Potentialdifferenz_ (JA Barth, 1889).
  * Sugiura (1993) M. Sugiura, IEE Proceedings A (Science, Measurement and Technology) 140, 443 (1993).
  * Lieberman and Lichtenberg (2005) M. A. Lieberman and A. J. Lichtenberg,  _Principles of Plasma Discharges and Materials Processing_ (John Wiley & Sons, 2005).
  * Belasri, Boeuf, and Pitchford (1993) A. Belasri, J. Boeuf, and L. Pitchford, Journal of Applied Physics 74, 1553 (1993).
  * Gao  _et al._ (2003) X. Gao, J. J. Liou, W. Wong, and S. Vishwanathan, Solid-State Electronics 47, 1105 (2003).
  * Gillman and Le May (2007) T. Gillman and I. Le May, Engineering Failure Analysis 14, 995 (2007).
  * Ashwell  _et al._ (1966) J. Ashwell, E. Cole, A. Pratt, and D. Sartorio, in  _1958 IRE International Convention Record_ , Vol. 5 (IEEE, 1966) pp. 72–77.
  * Phelps and Petrovic (1999) A. Phelps and Z. L. Petrovic, Plasma Sources Science and Technology 8, R21 (1999).
  * Boeuf and Pitchford (1991) J.-P. Boeuf and L. C. Pitchford, IEEE Transactions on Plasma Science 19, 286 (1991).
  * Birdsall (1991) C. K. Birdsall, IEEE Transactions on Plasma Science 19, 65 (1991).
  * Bird (1994) G. A. Bird,  _Molecular gas dynamics and the direct simulation of gas flows_ (Clarendon Press, Oxford, 1994).
  * Nanbu (2000) K. Nanbu, IEEE Transactions on plasma science 28, 971 (2000).
  * Werner  _et al._ (2018) G. R. Werner, T. G. Jenkins, A. M. Chap, and J. R. Cary, Physics of Plasmas 25, 123512 (2018).
  * Note (1) This approximation is subtly different from that in Werner  _et al._ (2018), but only when $\beta$ has an explicit time-dependence—a regime that has not yet been investigated and is irrelevant for this paper.
  * Jenkins, Werner, and Cary (2021) T. G. Jenkins, G. R. Werner, and J. R. Cary, arXiv e-prints , arXiv:2102.11988 (2021), arXiv:2102.11988 [physics.plasm-ph] .
  * Biagi (1999) S. Biagi, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 421, 234 (1999).
  * Yamabe, Buckman, and Phelps (1983) C. Yamabe, S. Buckman, and A. Phelps, Physical Review A 27, 1345 (1983).
  * Vahedi and Surendra (1995) V. Vahedi and M. Surendra, Computer Physics Communications 87, 179 (1995).
  * Kutasi and Donkó (2000) K. Kutasi and Z. Donkó, Journal of Physics D: Applied Physics 33, 1081 (2000).
  * Nieter and Cary (2004) C. Nieter and J. R. Cary, Journal of Computational Physics 196, 448 (2004).