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

Townsend discharge is the conduction of electric current in a gas due to a cascade of ionizing collisions.¹ 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 toward 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 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 Sec. II A.

Electrical discharge and breakdown have been extensively studied both for a multitude of applications, including lighting,³ plasma processing,⁴ and lasers,⁵ as well as to avoid detrimental effects, including damage to electrical devices,⁶ explosions in chemical processing facilities,⁷ and damage to aerospace systems.⁸ Following Townsend's description of the mechanism behind gas discharge,¹ various models and simulation methods have been employed to study the phenomena.⁹ Analytical models have proven effective for computing breakdown and steady-state current amplification in simple geometries;⁹ 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.^4,10 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.¹⁰ 

PIC methods evolve the distribution function $f(x,v,t)$ in time via the method of characteristics and, therefore, can simulate arbitrary distribution functions.¹¹ Collisions can be included in PIC using the PIC Monte Carlo Collisions framework (PIC-MCC).^11–13 PIC-MCC will be discussed further in Sec. III. PIC-MCC requires that the smallest mean free path be resolved by the grid cell size $Δx$⁠. Further, in PIC, the time step $Δt$ is limited according to $Δt≲Δx/vmax$⁠, 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 time step imposed by the maximum electron speed. This makes PIC-MCC simulations computationally slow. Because PIC-MCC is slow, variations of PIC simulation have been applied to discharge. One approach is treating the ions as immobile to eliminate the timescale mismatch,^14,15 but this technique fails when ion motion becomes relevant as is the case for Townsend discharge. Another technique is electron sub-cycling, in which the time step for the electron is smaller than that of the ion,¹⁶ but this technique yields at most a factor of two speed-up. Hybrid methods combine PIC and fluid approaches (e.g., with fluid electrons but PIC ions),¹⁰ 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 method¹⁷ 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 at the cathode and collisions. A detailed introduction to SLPIC is provided in Sec. II B. 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 Sec. III. The simulation design and parameters will be provided along with our methodology for determining breakdown in Sec. IV. In Sec. 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.

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 traveled is termed¹ Townsend's first ionization coefficient α. The average number of secondary electrons emitted per ion impact is Townsend's second ionization coefficient γ_se. The coefficient α depends on the voltage V across the gap and the collision cross sections. The coefficient γ_se depends on the ion species and the electrode material and treatment. When $γse(eα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 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

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 α can be expressed as

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 $α/ng$ as a function of the reduced field $E/ng$⁠) 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 γ_se taken from Lieberman and Lichtenberg.⁴ 

Different gases yield different Paschen curves due to different cross sections, ionization energies, and secondary emission yields at the cathode. In this work, we simulated a neutral argon gas. Argon has been studied extensively in the context of electrical breakdown, including simulation,^10,18,19 and experimental work.^9,20,21 Research has been focused on identifying the mechanisms responsible for breakdown and the regimes in which these mechanisms are relevant. In this work, we will compare our simulation results to experimental results compiled by Phelps and Petrovic.⁹ 

Before describing how collisions can be implemented within the SLPIC method in Sec. III, we review the derivation of SLPIC, following Werner et al.¹⁷ PIC numerically evolves a particle distribution $f(x,v,t)$ according to the Vlasov equation

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

and then make the critical “speed-limiting” approximation

This approximation is valid if either (1)$∂tf(x,v,t)≈0$ or (2)$β(v)≈1$⁠.²² Thus, we arrive at the “speed-limited” approximation to the Vlasov equation

In the special case of a steady state (i.e., $∂tf≡0$⁠), the speed-limiting approximation is exact and SLPIC, though different from PIC, is as accurate as PIC. For the special choice of $β(v)≡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 $[xp(t),vp(t)]$ with weights $wp(t)$ (a macro-electron with weight w_p represents w_p electrons, i.e., it has mass $wpme$ and charge $−wpe$⁠):

where S is the particle shape function¹¹ and δ 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

where we have assumed that $∂x·v+∂v·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 $β̇≡(d/dt)β(vp(t))=βa·∂vβ$⁠.

When $β(vp)=1$⁠, these are the usual (PIC) equations of motion. When $0<β(vp)<1$⁠, a SLPIC macroparticle, representing physical particles with true velocity $vp$⁠, travels with pseudovelocity $ẋp=β(vp)vp$ and pseudoacceleration $v̇p=β(vp)a(xp(t),vp(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 $Δt$⁠, the macroparticle moves $(vp,a)βΔt$⁠, while a physical particle moves $(vp,a)Δt$⁠. For example, a SLPIC electron in a uniform magnetic field will have the correct gyroradius, but its gyrofrequency will be (unphysically) reduced by a factor $βvp$⁠.

Choosing $β(vp)<1$ invokes the SLPIC approximation [Eq. (5)], reducing accuracy; however, it limits the simulated speed of macroparticles from v_p to $ẋ=β(vp)vp$⁠, which allows larger time steps, hence faster simulation. In most cases, the PIC time step must be smaller than $Δt≲Δx/vmax$⁠, where $Δx$ is the grid cell size and $vmax$ is the maximum particle speed. With SLPIC, a larger time step is allowed: $Δt≲Δx/[β(vmax)vmax]$⁠. Typically, we choose a speed limit $v0≪vmax$ and then define $β(vp)≡v0/vp$ to ensure that for all particles, $||ẋ||≤v0$⁠, so that we can use $Δt≈Δx/v0≫Δx/vmax$⁠. For particles with $vp<v0$⁠, however, we define $β(vp)≡1$ so that they experience “true” motion. In this way, we speed-limit only the particles (with $vp>v0$⁠) that need to be speed-limited to allow a large time step, $Δt=v0/Δx$⁠. It turns out that speed-limiting also reduces the plasma frequency so that it is resolved by the large $Δt$ and poses no instability.^17,23 For simplicity, we will use this specific $β(vp)$ for the rest of the paper: $β(vp)=v0/vp$ for $vp≥v0$⁠, and $β(vp)=1$ for $vp<v0$⁠, or, equivalently

where Θ is the Heaviside step function (and $vp=||v||$⁠).

As a macroparticle accelerates so that $vp>v0$⁠, its simulated speed $ẋp=βvp$ remains at v₀, and its weight $wp(t)$ decreases. The decrease in $wp(t)$ reflects the fact that the physical particle density represented by the macroparticle at time t decreases, and the macroparticle represents fewer physical particles at time $t+Δ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 $vp>v0$ experiences an acceleration and its true velocity v_p increases, but its pseudovelocity $ẋp=βvp=v0$ 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. While SLPIC may sometimes result in a more uniform macroparticle density, this is an unintended bonus and it may not always be the case. SLPIC is orthogonal to techniques for managing macroparticle weights/numbers for efficiency and accuracy.^24,25

Therefore, as v_p increases above v₀, $β(vp)$ decreases, and so too must w_p decrease. Equation (10) shows that $wp/β(vp)$ is constant over any particle trajectory; we define this constant to be $wj,p≡wp/β(vp)$⁠. Here, we choose the subscript j because $wj,p$ facilitates calculation of the flux density distribution $j(x,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, $j(x,v,t)=vf(x,v,t)$⁠, as follows:

To estimate the flux $j·n̂dA$ through some surface element $n̂dA$ from the macroparticles crossing it in time $Δt$⁠, we consider that a macroparticle with $ẋ$ will cross the surface if it is within a distance $ẋΔt·n̂$ of the surface. From Eq. (12), we see that, since $ẋ=β(vp)vp$⁠, the flux (integrated over time $Δt$⁠) through the surface element is simply the sum of the $wj,p$ overall macroparticles crossing the surface in $Δt$⁠.

Because the $wj,p$ is constant in time, a macroparticle always represents the same physical flux whenever it crosses a surface. If a macroparticle p has $wj,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 $wp(t)$ varies in time, it is nontrivial and very useful. For example, it ensures that, if a macroparticle enters 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 $wj,p$⁠. Consider a macroparticle with $vp=10v0, β(vp)=0.1, wj,p=100$⁠, and w_p = 10 and suppose it crosses a surface to enter volume V during time interval $Δt$⁠. That macroparticle represents w_p = 10 physical particles in volume V; however, it represented $wj,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 $(xp,vp)$ represents a swarm of w_p physical particles within the volume $d3x d3v$ around $(xp,vp)$⁠, and that those physical particles travel roughly with the macroparticle to $(xp+ẋpΔt,vp+v̇pΔt)$ over time $Δt$⁠. Fundamentally, however, a macroparticle represents a chunk of the distribution $fd3x d3v$⁠, and SLPIC decouples macroparticles from the physical particles they represent. A speed-limited SLPIC macroparticle moves in phase space with $ẋp$ and $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+Δt$⁠.

In the above example, the density represented by the SLPIC macroparticle corresponds to w_p = 10 physical particles that are near $[xp(t),vp(t)]$ at time t; however, over time interval $Δt$⁠, the physical flux represented by the macroparticle includes all the physical particles that would be near $[xp(t′),vp(t′)]$ at any time $t′∈[t,t+Δt]$⁠. Since the physical particles travel 10 times faster than the macroparticle, the physical flux is $10wp(t)/Δt=wj,p/Δt$⁠.

It is important to remember that SLPIC is accurate only if $f(x,v,t)$ changes sufficiently slowly—and it is as accurate as PIC in the steady-state limit. The $f(x,v,t)$ describing a single particle with $v>v0$ 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(x,v,t)$ changes too rapidly for the SLPIC approximation to be valid.

Exactly what constitutes “sufficiently slowly” for SLPIC accuracy remains an open question. We have previously observed that for (Langmuir) wave-particle interaction in 1D, “sufficiently slowly” means that the wave phase velocity must be slower than the SLPIC speed limit for accurate simulation.^17,23 Similarly, although studying SLPIC in magnetized plasmas is beyond the scope of this work, we imagine that the electron-cyclotron resonance interaction will be too fast for accurate simulation, as long as SLPIC reduces the electron gyrofrequency—but ion-cyclotron resonance might be accurately simulated. Although disadvantageous, it is precisely because SLPIC neglects fast motions that it can offer faster simulation. Importantly, in this paper we apply SLPIC to a steady-state problem in the sense that in the very beginning of breakdown, the space charge has negligible effect on the applied constant field; in the steady-state limit, we expect SLPIC to be essentially as accurate as PIC.

Reiterating an important point: the flux weight $wj,p=wp(t)/β(vp(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 $wj,p$ remains constant. If the impulse increases the particle's speed v_p, its density weight w_p must change accordingly so that $wp=β(vp)wj,p$⁠. In Sec. III, we will see that the flux weight is conserved when simulating collisions involving SLPIC particles.

Collisions can be simulated in PIC with a PIC-MCC algorithm;^11–13 the MCC algorithm can be used for SLPIC with just a few modifications: (1) using the flux weight $wj,p$ instead of the density weight w_p (cf. Sec. II B), and (2) considering the slowing of time when determining the collision rate. (Since most SLPIC simulations will have macroparticles with a distribution of $wj,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, in which a primary electron collides with a primary neutral atom, resulting in the scattering of the primary electron, the neutral “becoming” an ion, and the “creation” of a secondary electron. (Although there are other kinds of collisions, e.g., elastic or excitation, between electrons and neutrals that do not ionize, we consider only ionizing collisions for simplicity.) A primary electron macroparticle p with weight $wp0$ and velocity $vp0$ collides with a neutral macro-atom with $wn0$ and $vn0$⁠. The primary electron macroparticle will be split into two parts: one part with weight $wp1$ represents uncollided electrons and the other with weight $wp1′$ represents scattered primary electrons (recoiling from having just ionized atoms). Similarly, $wn1$ is the weight of uncollided neutrals, and w_i the weight of “scattered” neutrals, which are really ions and are now subject to acceleration in the applied electric field. In addition, secondary electrons are created with weight w_se. In principle, there could be many scattered/secondary macroparticles with different velocities, but for simplicity we will assume at most one per primary species. The uncollided particles retain the original velocities; the other particles must be assigned new velocities according to the ionization process: $vp1′, vi$⁠, and $vse$⁠. Conservation of subatomic particles requires $wp0=wp1+wp1′, wn0=wn1+wi$⁠, and $wp1′=wi=wse=Wc$⁠, 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 $Δt$ in a cell volume $ΔV$ between two arbitrary particle distributions, f_p (primary electrons) and f_n (neutrals) is

where $σ(vr)$ is the collision cross section, which depends on the relative velocity $vr$ between colliding particles. To calculate W_c between two PIC macroparticles in the same grid cell, we substitute the appropriate distributions, i.e., $fp(v)=(wp0/ΔV)δ3(v−vp0)$ and $fn(v)=(wn0/ΔV)δ3(v−vn0)$⁠, 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

This can be verified by considering a single electron moving through a swarm of atoms of density $nn0=wn0/ΔV$ with velocity $vn0$⁠. In the reference frame co-moving with the atoms, $d=||vn0−vp0||Δt$ is the distance traveled by the electron in $Δt$⁠, and $λmfp=(nn0σ)−1$ is its mean free path. Therefore, $d/λmfp$ is the expected number of collisions that a single electron would experience (for small $Δt$ such that $d/λmfp≪1$⁠). The total number of electron collisions is then $Wc=wp0d/λ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 $Wc=0.01wp0$⁠, 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 $Wc′=wp0$⁠. 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 $Wc>wp0$ (indicating that $Δt$ is too large). The only differences introduced by SLPIC are the determination of W_c and the resulting macroparticle weights.

As described in Sec. II B, 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 $wj,p0=wj,p1+wj,p1′, wj,n0=wj,n1+wj,i$⁠, and $wj,p1′=wj,i=wj,se=Wj,c$⁠. Because flux weights w_j are constant in time [unlike density weights $wp(t)$], this ensures that the fluxes of particles entering and exiting $ΔV$ are appropriately preserved, i.e., the flux of unscattered plus scattered primaries leaving $ΔV$ equals the flux of primaries entering $ΔV$⁠, and the fluxes of secondary electrons and scattered primaries leaving $ΔV$ are equal, etc.

The rule for determining the collision rate, $Wj,c/Δ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 $λmfp$⁠. Scattering causes the density of original beam particles to decrease as $∼ exp (−x/λ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 macroparticles must scatter with probability $d′/λmfp=ẋΔt/λmfp=β(v)vΔt/λmfp$ within time $Δt$ (in contrast, a physical particle would scatter with probability $d/λmfp=ẋΔt/λmfp=vΔt/λmfp$⁠). This can also be viewed as the reduction in collision rate due to the slowing down of time by a factor $β(v)$⁠. The collision rate (in terms of fluxes) in $ΔV$ must therefore be

where $βp(v)$ is the speed-limiting function for electrons. This merely says that—if one were to measure fluxes into and out of $ΔV$—the flux of scattered electrons (i.e., $Wj,c/Δt$⁠) would equal the flux of incident electrons (i.e., $wj,p0/Δt$⁠) times $d′/λmfp$⁠.

The above expression involves the neutral density: $nn0=wn0/ΔV=βn(vn0)wj,n0/ΔV$⁠. (In most applications, atoms will be slow and not speed-limited, hence $βn≡1$⁠, but we want this treatment to extend to arbitrary collisions.) Thus

This expression (or rather $Wj,c/Δt$⁠) yields the number of physical electrons scattered per $Δ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 time step is modified by a factor $βn(vn0)βp(vp0)$⁠.

After calculating $Wj,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 $wj,p$ is determined, its density weight is set by $wp=β(v)wj,p$⁠.

We modeled 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 $Δx$ was set to a quarter of the mean free path. The number of cells for each simulation is given in Table I. We ran both SLPIC and PIC simulations. The setup parameters for the SLPIC and PIC simulations were identical except for the time step. For PIC, we used a time step given by $Δt=Δx/vmax$⁠, where $vmax=2eV/me$ was the maximum electron speed. For SLPIC, we used a time step given by $Δt=Δx/v0$⁠, where $v0=vmaxme/mAr$ was the electron speed limit imposed by $β(vp)$ given in Eq. (11). This v₀ was chosen because it provided the maximum speed up without affecting the ion motion. The large mass difference between electrons and argon ions resulted in the SLPIC time step being 270 times larger than the PIC time step. Therefore, the SLPIC simulations required 270 times fewer time steps than PIC. The number of time steps for each simulation is given in Table I.

To initiate the discharge, we injected 100 electrons at the cathode in the first time step 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 database²⁶ and Phelps database²¹ on www.lxcat.net on August 15, 2020. We did not include any ion collisions. Although ion collisions could increase V_b, γ_se depends only weakly on the ion energy, so the effect of these collisions is negligible. 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.²⁷ For ionization collisions, the products were generated according to an algorithm developed by Kutasi and Donko.¹⁹ We used an energy-dependent model of secondary electron emission due to ion impact at the cathode.⁹ Secondary electrons were emitted from the cathode with zero velocity with a probability given by Eq. (17) that depends weakly on the energy of the incident ion ϵ_i

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 (⁠$302mArd/eE$⁠). For $V<Vb$⁠, the ion population decreases with time after an initial rise, eventually returning to zero. For $V>Vb$⁠, the ion population increases with time. This criterion is equivalent to the breakdown condition given by $γse(eαd−1)>1$⁠. 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. It can be seen that for $V>Vb$⁠, the ion population increases with time, and for $V<Vb$⁠, the ion population decreases with time after the initial rise. While SLPIC and PIC give the same results for V_b, the dynamics of the simulation are clearly different. The initial rise in ion population is steeper for PIC because the seed electrons cross the simulation domain and ionize the argon gas almost immediately, while in SLPIC, the speed-limited seed electrons cross the domain more slowly and therefore take longer to ionize the argon gas. The speed limiting also causes the ion population to evolve more smoothly as electrons and ions now have comparable velocities and oscillate in and out of the simulation domain at similar frequencies, but out of phase. SLPIC cannot capture the transient behavior of the particles since the speed-limiting approximation Eq. (5) is not satisfied in this regime. However, we do see convergence of the SLPIC dynamics to PIC as we increase the speed-limit.

We bracketed the breakdown voltage for 11 values of pd ranging from 0.3 to 300 Torr cm. We performed simulations with the Vorpal²⁸ 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 I.

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⁹ 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. At low pd, the experiments approach the vacuum discharge regime where processes such as outgassing, vacuum arc by the burned cathode, and flashover becomes relevant to breakdown.^15,29 We did not simulate these effects, which may account for the discrepancy. Paschen's law, given in Eq. (1), is plotted using the values of A = 11.5 cm⁻¹Torr⁻¹, B = 176 V cm⁻¹Torr⁻¹, and $γse=0.07$ taken from Lieberman and Lichtenberg.⁴ 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.

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 I. 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 I 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.

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 $Vb(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.

This work was supported by NSF Grant No. PHY1707430 and by U.S. Department of Energy SBIR Phase II Award No. 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.

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