{'fraction_non_alphanumeric': 0.04084665434974697, 'fraction_numerical': 0.007189332476504137, 'mean_word_length': 4.508628318584071, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computing the Paschen curve for argon with speed-limited particle-in-cell simulation

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. 1 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. 2 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 A.
Electrical discharge and breakdown have been extensively studied both for a multitude of applications, including lighting, 3 plasma processing, 4 and lasers, 5 as well as to avoid detrimental effects, including damage to electrical devices, 6 explosions in chemical processing facilities, 7 and damage to aerospace systems. 8 Following Townsend's description of the mechanism behind gas discharge, 1 various models and simulation methods have been employed to study the phenomena. 9 Analytical models have proven effective for computing breakdown and steady-state current amplification in simple geometries; 9 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. 10 PIC methods evolve the distribution function f (x, v, t) in time via the method of characteristics and, therefore, can simulate arbitrary distribution functions. 11 Collisions can be included in PIC using the PIC Monte-Carlo Collisions framework (PIC-MCC). [11][12][13] PIC-MCC will be discussed further in §III. PIC-MCC requires that the smallest mean free path be resolved by the grid cell size ∆x. Further, in PIC, the timestep ∆t is limited according 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, 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 timestep for the electron is smaller than that of the ion, 16 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), 10 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 17 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 §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 §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.

A. 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 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 α/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 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 9 .

B. SLPIC
Before describing how collisions can be implemented within the SLPIC method in §III, we review the derivation of SLPIC, following Werner et al. 17 . 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) Thus we arrive at the "speed-limited" approximation to the Vlasov equation In the special case of a steady state (i.e., ∂ t f ≡ 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 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): where S is the particle shape function 11 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-spacevolume preserving (Liouville's theorem proves this for any Hamiltonian motion). We have When β(v p ) = 1, these are the usual (PIC) equations of motion. When 0 < β(v p ) < 1, a SLPIC macroparticle, representing physical particles with true velocity v p , travels with 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 (v p , a)β∆t, while a physical particle moves (v p , 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 β(v p ) < 1 invokes the SLPIC approximation [Eq. (5)], reducing accuracy; however, it limits the simulated speed of macroparticles from v p toẋ = β(v p )v p , which allows larger timesteps, hence faster simulation. In most cases, the PIC timestep must be smaller where ∆x is the grid cell size and v max is the maximum particle speed. With SLPIC, a larger timestep is allowed Typically, we choose a speed limit v 0 ≪ v max , and then define β(v p ) ≡ v 0 /v p to ensure that for all particles, we define β(v p ) ≡ 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, ∆t = v 0 /∆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 β(v p ) for the rest of the paper: where Θ is the Heaviside step function (and v p = v ).
As a macroparticle accelerates so that v p > v 0 , its simulated speedẋ p = β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 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 v p > v 0 experiences an acceleration and its true velocity v p increases, but its pseudo-velocityẋ p = β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. 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 macroparticles weights/numbers for efficiency and accuracy. 24,25 Therefore, as v p increases above v 0 , β(v p ) decreases, and so too must w p decrease. Equation (10) shows that w p /β(v p ) is constant over any particle trajectory; we define this constant to be w j,p ≡ w p /β(v p ). Here we choose the subscript j because w j,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 ·ndA through some surface elementndA 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ẋ = β(v p )v p , the flux (integrated over time ∆t) through the surface element is simply the sum of the w j,p over all macroparticles crossing the surface in ∆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 , β(v p ) = 0.1, 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 > 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 (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 w j,p = w p (t)/β(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 = β(v p )w j,p .
In Sec. III, 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; 11-13 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 B), 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, 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 w p0 and velocity v p0 collides with a neutral macroatom with w n0 and v n0 . The primary electron macroparticle will be split into two parts: one part with weight w p1 represents uncollided electrons and the other with weight w ′ p1 represents scattered primary electrons (recoiling from having just ionized atoms). Similarly, w n1 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: v ′ p1 , v i , and v se . Conservation of subatomic particles requires w p0 = w p1 + w ′ p1 , w n0 = w n1 + w i , and w ′ p1 = w i = w se = 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 ∆t in a cell volume ∆V between two arbitrary particle distributions, f p (primary electrons) and f n (neutrals) is where σ(v r ) is the collision cross-section, which depends on the relative velocity v r between colliding particles. To calculate W c between two PIC macroparticles in the same grid cell, This can be verified by considering a single electron moving through a swarm of atoms of density n n0 = w n0 /∆V with velocity v n0 . In the reference frame co-moving with the atoms, ∆t is the distance traveled by the electron in ∆t, and λ mfp = (n n0 σ) −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 W c = w p0 d/λ 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 = 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 ∆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 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 w j,p0 = w j,p1 + w ′ j,p1 , w j,n0 = w j,n1 + w j,i , and w ′ j,p1 = w j,i = w j,se = 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 ∆V are appropriately perserved; 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, W j,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 macroparticle 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., W j,c /∆t) would equal the flux of incident electrons (i.e., w j,p0 /∆t) times d ′ /λ mfp .
(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 W j,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 timestep is modified by a factor β n (v n0 )β 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 = β(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 ∆x was set to a quarter of the mean free path. The number of cells for each simulation is given in Table I 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 28 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 in Fig. 2. For the values of pd where PIC was run, PIC and SLPIC classified each simula-

AIP PUBLISHING DATA SHARING POLICY
The data that support the findings of this study are available from the corresponding author upon reasonable request.