This paper discusses temporally continuous and discrete forms of the speed-limited particle-in-cell (SLPIC) method first treated by Werner *et al.* [Phys. Plasmas **25**, 123512 (2018)]. The dispersion relation for a 1D1V electrostatic plasma whose fast particles are speed-limited is derived and analyzed. By examining the normal modes of this dispersion relation, we show that the imposed speed-limiting substantially reduces the frequency of fast electron plasma oscillations while preserving the correct physics of lower-frequency plasma dynamics (e.g., ion acoustic wave dispersion and damping). We then demonstrate how the time step constraints of conventional electrostatic particle-in-cell methods are relaxed by the speed-limiting approach, thus enabling larger time steps and faster simulations. These results indicate that the SLPIC method is a fast, accurate, and powerful technique for modeling plasmas wherein electron kinetic behavior is nontrivial (such that a fluid/Boltzmann representation for electrons is inadequate) but evolution is on ion timescales.

## I. INTRODUCTION

The speed-limited particle-in-cell (SLPIC) method^{1} is a relatively new plasma modeling technique. It is most suitable for discharges in which the physics of interest occurs on relatively slow timescales (e.g., ion transport/profile relaxation) but is nevertheless tied to kinetic electron behaviors that a fluid/Boltzmann model cannot capture (e.g., distribution function modifications from neutral collisions or sheath interactions, or Landau damping). In such simulations, one is typically constrained to model both the heavy, slow ion species and the light, fast electrons using conventional particle-in-cell (PIC) techniques. The ensuing computational costs can be (possibly unaffordably) high; simulation time steps must resolve the electron plasma frequency since kinetic electrons are present, but ion timescales of interest may exceed the electron oscillation period by many orders of magnitude.

In the SLPIC approach, conventional PIC is modified to artificially slow down “fast” behaviors which are numerically troublesome, despite being physically unimportant for the physics of interest. Larger simulation time steps can thus be used while retaining the detailed physics behaviors associated with the slower, longer timescales. The specifics of the SLPIC method will be explained in a later section of this paper, but we will note here that numerical experiments using SLPIC simulations to model sheath formation in an argon plasma have shown that remarkable speedup factors (160 times faster than conventional PIC methods^{2}) can be achieved.^{1} When SLPIC can be appropriately used for modeling, it is both accurate and powerful.

A concept understood since the early days of PIC modeling is that while one may “recover more of the essence of the situation being simulated by changing the interaction laws,” such changes are accompanied by costs: “the more one meddles with the “laws” of nature the more one must understand the consequences.”^{3} While this quote in its original context refers to the various approximations used in PIC simulations (e.g., finite-sized particles, grid spacings, and time steps), in this paper, we explore its relevance to SLPIC. Although SLPIC is in many ways similar to conventional PIC, the ways in which it is different to introduce additional effects that one must understand in order to have confidence in its provided solutions—and this is true even independent of any effects imparted by finite-sized particles, grid spacings, and time steps (though such effects are not unimportant). Fundamentally, SLPIC and PIC methods both seek to statistically approximate the evolution of smooth particle distribution functions in a multidimensional phase space, in response to self-consistent fields and forces—but the underlying evolution equations of the two methods are different and will yield different physics (e.g., linear plasma wave dispersion) even before any particle-based approximations are made.

In this work, therefore, we focus on developing an understanding of the behavior of a plasma evolving with speed-limited dynamics (hereafter SLD). In SLD, the plasma is governed by continuous, “SLPIC-like” equations of motion that differ from the ones governing plasma evolution in our universe, but which approximate them in certain limits that we will quantify. For purposes of comparison, we also designate the dynamics of plasma evolution in our universe as “ordinary dynamics” (OD). SLPIC and PIC simulation methods are, respectively, the discrete numerical analogues of the SLD and OD that we will explore. [Alternatively, one can think of SLD or OD, respectively, as continuous limits of SLPIC or PIC, wherein both the time step $\Delta t$ and the grid spacing $\Delta x$ approach zero as the velocity dependence of the distribution function $f\alpha $ becomes smooth.] We will show that some well-understood physics processes from OD persist in SLD, and that other processes are substantially modified, some of them in very helpful ways.

More specifically, in this paper, we derive and analyze the analytical dispersion relation in an electrostatic, collisionless, and unmagnetized plasma evolving according to modified Vlasov–Poisson equations that govern SLD. Such a plasma will have different wave modes, with different dispersion, relative to a real plasma evolving with ordinary dynamics (OD). By comparison with the dispersive behavior that arises in the conventional (OD) Vlasov–Poisson system, we demonstrate both analytically and numerically that the speed-limiting of SLD can quantifiably modify high-frequency behaviors of this plasma (electron plasma oscillations) while leaving low-frequency motion [ion acoustic wave (IAW) decay via electron Landau damping] undisturbed.

Section II of this paper explains the SLD concept, together with its connections both to the SLPIC algorithm and to the conventional kinetic theory of OD. In Sec. III, we discuss the behavior of an electrostatic ion–electron plasma that evolves with SLD and derive the dispersion relation associated with this plasma. Section IV contains an analysis of the various waves permitted by this dispersion relation, together with the new behaviors imparted by SLD relative to known OD behaviors. We demonstrate that the speed-limiting significantly relaxes a fundamental numerical constraint associated with conventional PIC methods and makes faster numerical simulations possible. Section V then considers the spatial fluctuation spectrum associated with SLD and shows it to be the same as that of OD; we briefly compare SLPIC and PIC simulations to demonstrate this point. Finally, in Sec. VI, we summarize our findings, review additional research directions which these findings might enable, and discuss various applications for the SLPIC concept in plasma modeling more generally.

## II. PHASE SPACE EVOLUTION AND ITS CONNECTION TO PIC AND SLPIC METHODS

To kinetically model plasma with OD (i.e., using the familiar physics of the real world), we first consider the self-consistent evolution of a distribution function $f\u0302\alpha (x,v,t)$ of physical particles of species *α*. This distribution evolves according to a phase-space continuity equation

where $a(x,v,t)$ is the self-consistent Lorentz acceleration experienced by a physical particle in the distribution at position **x** with velocity **v** at time *t*, in response to local (microscopic) electromagnetic fields. For Hamiltonian systems, the additional phase-space preserving constraint $(\u2202/\u2202x)\xb7v+(\u2202/\u2202v)\xb7a=0$ of Liouville's theorem allows us to rewrite Eq. (1) in the familiar Klimontovich form

which can be formally solved by the method of characteristics. Along characteristic trajectories

the value of the distribution function $f\u0302\alpha [xj\alpha (t),vj\alpha (t),t]$ is preserved; writing the distribution function in the form

solves Eq. (2) and captures the detailed microscopic behavior of each of the $N\alpha p$ particles in response to the electromagnetic fields they (and the other species in the system) produce.

For realistic physical particle counts, this microscopic behavior is far too detailed to simulate numerically in most plasmas, and the singular nature of Eq. (5) is likewise problematic. If, instead, one passes to the continuum limit (subdividing the discrete particle charges and masses in a manner that preserves the local volumetric charge, mass, and energy content), one arrives at the Vlasov equation

which describes the evolution of a continuous (nonsingular) phase space fluid $f\alpha $ from which the effects of particle discreteness have been removed. More detailed discussion of this transition, which considers ensemble averages of Eq. (2), two-particle and higher-order correlation terms, collision operators,^{4} etc., has been considered by other authors,^{5–7} but Eq. (6) suffices for our purposes here. Like Eq. (2), it can also be solved by the method of characteristics; with trajectories evolving according to Eqs. (3) and (4), we may formally write

and can verify that it is a solution of Eq. (6). The smooth distribution is now represented as a set of $N\alpha $ discrete macroparticles which evolve along the trajectories given by Eqs. (3) and (4). The weight function $wj\alpha $ is representative (in some statistical sense) of the local value of $f\alpha $ in a region near the particle's initial point on the phase space trajectory. PIC simulation techniques build upon the fundamental concept that Eqs. (3), (4), and (7) solve Eq. (6), and that a sufficiently large number $N\alpha $ of macroparticles, distributed so as to adequately resolve the relevant regions of the phase space, can statistically represent the 6D + time evolution of the smooth distribution function $f\alpha (x,v,t)$.

In this paper, we will consider Eqs. (3), (4), and (6) as the fundamental equations describing plasma evolution under OD.

Techniques for mapping the equations of PIC onto discrete computational time steps and finite grids (broadening the spatial extent of macroparticles from delta-functions to small-but-finite widths) are discussed extensively in the existing literature;^{3,8–13} detailed explanations and/or derivations of such techniques will not be discussed here except as needed. For the present, it suffices to note (as the previously cited works discuss) that finite-sized grid cells and time steps impose a number of constraints on conventional explicit PIC simulations:

The

*Debye length resolution constraint*: that a representative grid cell size $\Delta x$ should adequately resolve the Debye length $\lambda D\alpha $ associated with any of the species in the simulation in order to avoid numerical heating effects [the Debye length of species*α*is defined as $\lambda D\alpha 2\u2261\u03f50T\alpha /(q\alpha 2n\alpha )$, where ${q\alpha ,n\alpha ,T\alpha}$ are the species charge, density, and temperature (in units of energy) and*ϵ*_{0}is the permittivity of free space];The

*cell-crossing-time constraint*: that the distance traveled by any macroparticle in the simulation during a finite time step $\Delta t$ should not exceed a representative grid cell size $\Delta x$, so that forces experienced by a particle during a single simulation time step are adequately resolved; andThe

*plasma oscillation constraint*: that the plasma frequency*ω*constrains the time step through the relation $\omega p\Delta t\u22642$; otherwise, numerical instability of these high-frequency oscillations ensues [the plasma frequency is defined as $\omega p2\u2261\u2211\alpha \omega p\alpha 2$, with the species plasma frequency $\omega p\alpha $ defined through $\omega p\alpha 2\u2261q\alpha 2n\alpha /(\u03f50m\alpha )$. Here, $m\alpha $ is the mass of species_{p}*α*and the other quantities were defined previously]. To be precise, this constraint arises from a more general requirement that every plasma mode frequency must be resolved by the simulation time step. But for a large class of problems, including the ones considered in this work, the highest mode frequencies are on the order of*ω*._{p}

These constraints can impose significant restrictions on a plasma simulation. Low temperatures and/or high densities decrease the Debye length and the allowable grid size, necessitating the use of finer grids and smaller time steps. Further, when both cold, massive ions and hot, light electrons are simulated with PIC, the time steps imposed by the plasma oscillation constraint (now dominated by fast electron motion since $\omega p\u223c\omega pe$) are so small that ions may hardly move at all in that time interval. Numerical techniques such as subcycling,^{15} in which ions are pushed less frequently and with a larger effective time step, can provide minor computational savings, but this gains one at most a factor of $\u223c2$ in speedup (for typical cases with comparable electron and ion particle counts) since the time step constraints arising from electron motion remain. For simulations where many periods of harmonic ion motion are of interest, the number of time steps required can be enormous.

The speed-limited particle-in-cell approach, and the more general speed-limiting concepts we explore in this work, is motivated by a desire to relax some of these constraints. The key idea is simple: the fastest particles and highest-frequency wave phenomena necessitate the smallest time steps, and if these fast particle and wave motions can be slowed, the time step constraints can be relaxed.

Accordingly, we introduce speed-limited dynamics (SLD), wherein equations from the derivation of the SLPIC method presented in Ref. 1 govern the plasma dynamics. Here, a distribution function $f\alpha $ of species *α* evolves as prescribed by a modified Vlasov equation

wherein a speed-limiting function $\beta (v)$ in the range $(0,1]$ has been introduced. This function transitions from values at or near unity (for “slow” particles) to values approaching $v0/|v|$ (for “fast” particles), and we can understand its effect by looking at the characteristic trajectories of the modified Vlasov equation

The product $|v|\beta (v)$, the speed at which an element of phase space changes its position $x(t)$ as it moves through the phase space, now has value $\u223c|v|$ at low velocities but is limited to value *v*_{0} (in the original vector direction of motion) at high velocities. We will hereafter refer to *v*_{0} as the “speed limit;” it is the upper bound on the rate at which motion in the position coordinate of phase space may proceed. Accordingly, some nuance is required in discussing the phase space evolution since the meaning of “velocity” is now ambiguous. An element of phase space has both a “true velocity” (the phase space coordinate **v**) and a “pseudo-velocity” $dx/dt=\beta v$ (the speed and direction at which it is permitted to move from one physical space coordinate to another).^{14} In a given time interval, elements in the phase space with large true velocity $v(t)$ [so that $\beta (v)\u226a1$] experience both smaller pseudo-velocities [Eq. (9)] as they move through the space and smaller changes to these pseudo-velocities (pseudo-acceleration) in response to applied forces [Eq. (10)]. Elements in the phase space with small true velocity experience no speed-limiting [$\beta (v)\u223c1$] and evolve in the same manner as their OD counterparts, as in Eqs. (3) and (4). Transitions across the boundary $|v|=v0$ in either direction are well-defined; this has already been demonstrated for SLPIC in Fig. 4 of Ref. 1, wherein particles are not observed to “pileup” at the boundary.

In this paper, we will consider Eqs. (8)–(10) as the fundamental equations describing plasma evolution under SLD.

Solutions to the SLD Vlasov equation, Eq. (8), can be represented statistically in the same manner as outlined above, setting

for a suitably large number $N\alpha $ of macroparticles evolving along the characteristic trajectories described by Eqs. (9) and (10). This is the speed-limited particle-in-cell method we have presented in previous work.^{1} However, this work will not focus on PIC or SLPIC implementations of Eq. (6) or Eq. (8). Instead, we will consider these equations analytically.

## III. A 1D1V ELECTROSTATIC PLASMA MODEL

In this section, we will apply the Vlasov equation of OD and the modified Vlasov equation of SLD to model dispersion in an electrostatic, unmagnetized plasma with a single ion species, in one spatial dimension and one velocity-space dimension. We will use the analytic forms of these equations (foregoing for the moment any discussion of the effects of finite time steps, grid spacings, or particle sizes) to ensure that we understand the new physics that the imposed speed-limiting of SLD imparts. Each species will use the kinetic equation

In SLD, we will use a speed-limiting function of form

wherein *H*(*x*) is the Heaviside function. An equivalent representation for this speed-limiting function is

In OD, $\beta (v)=1$ (the $v0\u2192\u221e$ limit of the SLD). The species couple via the Poisson equation

Many equilibrium solutions of Eq. (12) are possible. We will choose physically reasonable solutions that are stationary Maxwellians in each of the individual species (though we will allow the two species to have different temperatures and neglect both the collisional processes that have brought the individual species to their present state and the interspecies collision processes that would further relax the system to a single temperature). Formally, we write the distribution function of species *α* as

where *n*_{0} is a species-independent constant number density, $T0\alpha $ is the constant temperature of species *α*, and $m\alpha $ is the species mass. With these equilibrium distributions, the corresponding equilibrium potential $\varphi 0(x,t)$ is a constant that can be set to zero.

Linearizing Eqs. (12) and (13) in perturbed quantities, and Fourier transforming from spacetime coordinates {*x*, *t*} to wavenumber and frequency coordinates ${k,\omega}$, yields the result

We obtain, for the perturbed distribution functions,

which we can then substitute into Eq. (18). The ensuing integrals will be undefined for resonant velocities $v\beta (v)=\omega /k$; we will implicitly stipulate that the integrals are to be evaluated using the Landau contour (traversing below any singularity) to retain both causality and the resonant physics. (Formally, this can be shown to be equivalent to the use of a Laplace transform, rather than a Fourier transform, in the time domain; it also permits the generalization of *ω* to complex values.) We obtain the integral relation

which describes the dispersive wave behavior of the plasma in both OD [where $\beta (v)=1$] and SLD [where Eq. (14) defines $\beta (v)$].

In OD, the integral in Eq. (20) can be expressed in terms of the plasma dispersion function,^{16} defined as

for $Im(\zeta )>0$ and by its analytic continuation for $Im(\zeta )\u22640$. In SLD, this integral is more complicated; given our choice for $\beta (v)$, there are both high-velocity regions of integration wherein $v\beta (v)=\xb1v0$, a constant (the speed limit), and low-velocity regions wherein $v\beta (v)=v$. We may represent the integrals over these various regions in terms of the complementary error function and the incomplete plasma dispersion function. The latter function was introduced by Franklin^{17} and its properties have been discussed extensively by Baalrud;^{18} it takes the form [generalized from Eq. (21)]

for $Im(\zeta )>0$ and by its analytic continuation for $Im(\zeta )\u22640$. It will be useful to note that $Z(\u2212\u221e,\zeta )=Z(\zeta )$ and that $Z(\u221e,\zeta )=0$. Additional properties of this function are provided in Appendix A.

In terms of these functions, we may rewrite Eq. (20) in the form $DSLD(k,\omega ;v0)=0$, where we define

Here, we make use of the parameters $\gamma \alpha \u2261v0/(2vt\alpha )$, a measure of the relative speed-limiting of species *α*; $\zeta \alpha =\omega /(2kvt\alpha )$, the conventional (and complex) argument of the plasma dispersion function; $\lambda D\alpha 2\u2261\u03f50T0\alpha /(qi2n0)$, the Debye length of species *α*; and $vt\alpha 2\u2261T0\alpha /m\alpha $, the thermal velocity for species *α*. The complementary error function is related to the conventional error function by $erfc(x)\u22611\u2212erf(x)$. The expression $DSLD(k,\omega ;v0)=0$ is the plasma dispersion relation of SLD, and its analysis and solutions will be the topic of Sec. IV.

Recalling that $erfc(\u221e)=0$ and the limits of the incomplete plasma dispersion function discussed above, we can show that $DSLD(k,\omega ;\u221e)=DOD(k,\omega )$, where the latter function has the explicit form

a standard result from elementary plasma kinetic theory.^{19,20} This is the plasma dispersion relation of OD, against which solutions with SLD will be compared.

## IV. ANALYSIS OF THE DISPERSION RELATION

We now consider various limits of $DOD(k,\omega )$ and $DSLD(k,\omega ;v0)$, together with the waves that arise in these various limits. For clarity, we will consider the physics of OD first and will then explore the changes which the speed-limiting of SLD imparts to these familiar processes.

### A. Plasma oscillations

The OD dispersion relation, Eq. (24), admits approximate analytic solutions corresponding to cold-plasma oscillations in the $\zeta \alpha \u226b1$ limit. The asymptotic expansion of $Z(\zeta \alpha )$ in this limit

can be substituted into Eq. (24) to obtain, at the lowest order,

where $\omega p\alpha 2=qi2n0/(\u03f50m\alpha )$ is the square of the plasma frequency of species *α*. These oscillations are dominated by electron motion (the ion term is order $me/mi$ smaller than the electron term); we may write the solution as $\omega 2=\omega pe2(1+me/mi)\u2261\omega p2$. We anticipate the prospect of significant changes to these oscillations in SLD, since the speed-limiting will preferentially modify the fast electron motion.

What is the behavior of the SLD dispersion relation, Eq. (23), in the same $\zeta \alpha \u226b1$ limit? Using the asymptotic expansions for $Z(\gamma ,\zeta )$ presented in Ref. 18 and summarized in Appendix A, we can show that the SLD version of Eq. (26) takes the form

which reproduces Eq. (26) in the $v0\u2192\u221e\u2009(i.e.\u2009\gamma \alpha \u2192\u221e)$ limit. It admits the solutions

The behavior of the function $h(\gamma \alpha )$ is shown in Fig. 1. For small $\gamma \alpha ,\u2009h(\gamma \alpha )\u223c2\gamma \alpha 2$, while for large $\gamma \alpha $, the species is not appreciably speed-limited and $h(\gamma \alpha )\u22481$. Recalling that $\gamma \alpha =v0/2vt\alpha $, it will be instructive to express the argument of *h* strictly in terms of the electron slowing-down parameter *γ _{e}*; we have $\gamma i=\gamma eTemi/(Time)$. The presence of the ion–electron mass ratio suggests that unless electrons are much cooler than ions, the quantity

*γ*will always be much larger than

_{i}*γ*. Accordingly, we may generally choose

_{e}*v*

_{0}in a way that alters electron behavior but not ion behavior, consistent with the intentional slowing down of the fastest particles in SLD while leaving slower particles undisturbed. Such a

*v*

_{0}will be faster than most ions but much slower than most electrons. For such a choice ($v0=4vti$), Fig. 1 also shows values of $\gamma e,\gamma i$ and their corresponding

*h*values for a hydrogen plasma thermalized to 10 eV. In this plasma, ion motion associated with plasma oscillations does not differ appreciably between SLD and OD, but electron motion is considerably modified.

Since it is primarily the speed-limiting of electrons that affects the oscillation frequencies predicted by Eq. (28), we examine the behavior of these modified electron plasma oscillations as a function of the ratio $v0/vte$ (providing intuition as to which velocities in a typical electron distribution, e.g., a Maxwellian, are restricted by the speed-limiting). The normalized oscillation frequency of Eq. (28) is shown in Fig. 2 for a monatomic hydrogen plasma with equal electron and ion temperatures. While the oscillation frequencies of SLD and OD are identical at large values of $v0/vte$ (minimal speed-limiting), reducing this ratio, and hence increasing the corresponding fraction $erfc(\gamma e)$ of speed-limited electrons, reduces the SLD frequency monotonically. As $v0/vte\u21920$, the linear approximation to the plasma frequency approaches the heuristic estimate made in Ref. 1: $\omega /\omega p\u223cv0/vte$. The high-frequency oscillations of OD have been mapped to lower-frequency oscillations in SLD by the speed-limiting of electrons.

We conclude that in SLD, the frequency of conventional electron-dominated plasma oscillations is reduced relative to OD. In Sec. IV C, we will consider how this reduction relaxes the “plasma oscillation constraint” referred to in Sec. II. But we must first determine whether the speed-limited dynamics preserves physics associated with lower-frequency waves.

### B. Ion acoustic waves

The OD dispersion relation, Eq. (24), also admits approximate analytic solutions corresponding to ion acoustic waves. Physically, these solutions are associated with “hot” electrons ($\zeta e\u226a1$) and “cold” ions ($\zeta i\u226b1$), and in these limits, the dispersion relation can be written in the approximate form

where the multiplicative factor of $k2\lambda De2$ simplifies the algebra and where *c _{s}*, the sound speed, satisfies the relation $cs2=T0e/mi$. Assuming that the imaginary part of the (now complex) $\omega =\omega r+i\omega i$ is small, we can perform a Taylor expansion

to show that this equation permits damped wavelike solutions of the form

where *ϵ* is again the electron/ion mass ratio. These are the ion acoustic wave (IAW) modes, which are essentially longitudinal compressions of the ion mass density that decay via Landau damping.

What is the corresponding behavior of the SLD dispersion relation? In the $\zeta i\u226b1,\zeta e\u226a1$ limit, Eq. (23) can be written with the same multiplicative factor in the form

Repeating the Taylor expansion procedure above, we can show that the ion acoustic wave dispersion relation of SLD has approximate analytic solutions

wherein

This is the SLD equivalent to Eq. (31). The OD and SLD forms are approximately equivalent provided that *γ _{i}* is sufficiently large (so that $h(\gamma i)\u223c1$, see Fig. 1) and $|\omega r/k|<v0$ (so that the speed limit

*v*

_{0}exceeds the phase velocity of the ion acoustic wave). Since the approximation $h(\gamma i)\u22481$ holds to one part in 10

^{4}or better for $v0>4vti$, we will write these constraints as a single condition

When *v*_{0} is chosen to satisfy this condition, and when the other conditions we assumed in the derivation [$\zeta i\u226b1,\zeta e\u226a1,Im(\omega )\u226aRe(\omega )$] are valid, the effects of SLD on the propagation and damping of the IAW are minimal.

What happens when the phase velocity condition is violated, such that $v0<|\omega r/k|$? In this case, the explicitly imaginary terms in Eq. (32) (proportional to Heaviside functions) vanish, and the solutions admitted now only capture the real part of Eq. (33). While the IAW still propagates, its Landau damping is not correctly modeled. This is consistent with a result that we have demonstrated in the previous work,^{1} namely, that the correct dynamics of resonant wave-particle interactions cannot be captured by speed-limited particle-in-cell simulations when particles whose velocities were previously synchronous with the wave velocity are speed-limited. Particles whose speed-limiting renders them unable to keep up with the wave cannot exchange energy with it, so the dissipative effects which lead to wave damping are effectively turned off. Although nothing in principle prevents us from choosing a smaller *v*_{0} value that violates Eq. (35), the inherent advantage of the velocity-dependent speed-limiting approach (preservation of IAW physics) would be lost in doing so.

The analytic form of the IAW above is approximate. The dispersion relation can also be solved numerically to find the complex *ω* associated with a given *k*, and we have done so for a plasma with density $n0=5.0\xd71016\u2009m\u22123,\u2009Te0=10$ eV, and $Ti0=1/40$ eV. In the supplemental materials for Ref. 18, Matlab algorithms for evaluating $Z(\zeta ),\u2009Z(\gamma ,\zeta )$ and their derivatives have been provided. We have made use of these algorithms and Matlab minimization routines to compute ${\omega ,k}$ values which satisfy the OD and SLD dispersion relations within a given tolerance, according to norm-minimization criteria

Here, the tolerance parameter $\delta \u226a1$ is a small positive number which constrains the allowable error in the numerical solution of a particular dispersion relation. In these computations, we fix *k* (and, for SLD, *v*_{0}) and then numerically evaluate the function *D* for various complex values of *ω*. Exact solutions to the dispersion relation have *D* = 0; we vary *ω* in the complex plane in a manner that seeks to minimize the norm of *D* and thus to approach these exact solutions. For *ω* values near an exact solution, this norm can, in principle, be reduced to be no greater than *δ*.

Various values for the speed-limiting parameter *v*_{0} can be chosen to assess the effect of speed-limiting on IAW behavior. Depending on the value of *v*_{0} and *k*, to obtain numerical convergence, it is sometimes necessary to raise *δ _{SLD}* to values as high as $10\u22123$, but solutions can usually be found for $\delta OD<10\u22126$ and $\delta SLD<10\u22125$ (and often for

*δ*values that are several orders of magnitude smaller, when $k\u226b1$).

In Fig. 3, we show the real and negative-imaginary parts of the computed frequency *ω* for various values of wavenumber *k* spanning several orders of magnitude. In this figure, the speed limit *v*_{0} has the value $0.2vte$; the analytic approximation to the dispersion relation (valid for $k\lambda De\u22721$) is also shown. The speed-limiting does not affect the IAW behavior appreciably; for modes whose wavelengths are large compared to the Debye length (i.e., where the analytic approximation is valid), the SLD real frequency *ω* is about half a percent low relative to the OD value and the SLD damping rate also drops by about 3%. For shorter-wavelength modes, the effect of the speed-limiting on both the real frequencies and damping rates is negligible. Nevertheless, the frequency of the modified plasma oscillations (from Fig. 2) is decreased to about $0.2\omega pe$ for this case.

For more restrictive speed limits *v*_{0}, greater deviation of IAW frequencies and damping rates from their non-speed limited (OD) values is observed, though still only for modes whose wavelengths are long compared to the Debye length. Generally speaking, real frequencies are shifted downward as the speed limit is decreased, though only by a few percent (<9% for $v0=0.05vte=2.1cs$ and <2% for $v0=0.1vte=4.2cs$). Damping rates also (generally) decrease in magnitude as the speed limit is decreased, but the magnitude of the relative decrease is larger (<42% for $v0=0.05vte=2.1cs$ and <13% for $v0=0.1vte=4.2cs$). In Figs. 4 and 5, we have plotted the variation of the IAW real frequency and damping rate of SLD as a function of the speed-limiting velocity normalized to the wave phase velocity, so as to more generally quantify the effect of the speed-limiting on the IAW dispersive behavior. These figures illustrate that lowering the speed limit has relatively little effect on either the real IAW oscillation frequency or damping rate as long as the speed limit is higher than the wave phase velocity. For longer-wavelength modes $(k\lambda De\u22721)$ in SLD, speed limits that approach the wave phase velocity generally give rise to modest reductions in the mode frequency (Fig. 4) and more pronounced reductions of the IAW damping rate (Fig. 5). This failure to properly capture the Landau damping of the IAW arises because portions of the distribution function that resonate with the wave (in OD) are prevented from doing so by the imposed speed-limiting.^{1} As we transition to shorter-wavelength modes (moving up the graph legend), absolute damping rates are increased to become comparable to the real mode frequency, while phase velocities are reduced to be of order *v _{ti}*. For these waves, speed-limiting does not appreciably influence the dynamics since the smallest sensible speed limit ($v0\u223c$ $4vti$, so as not to speed limit the bulk ion distribution) still exceeds $v\varphi =\omega /k$ for large

*k*.

Accordingly, we may assert that SLD preserves the physics of IAW propagation and damping provided that the condition in Eq. (35) holds, namely, when ions are not speed-limited and when the speed limit reasonably exceeds the phase velocity of the IAWs in the system. At the same time, this speed-limiting considerably reduces the frequency of plasma oscillations (as was shown in Sec. IV A).

### C. Normal modes and the plasma oscillation constraint

Having explored the normal modes of our 1D1V plasma, we now revisit the numerical constraints which the use of a finite time step $\Delta t$ will impose on particle-in-cell simulations of this plasma.

We have noted in Sec. II that numerical instability will ensue in a PIC simulation if the frequency of any plasma mode is not resolved, and that, in particular, we must resolve the frequency associated with electron plasma oscillations. These oscillations are generally the highest-frequency modes in a PIC simulation containing electrons (because ions introduce only small corrections to the approximation $\omega p\u2248\omega pe$ when the electron–ion mass ratio is small). Accordingly, any numerical simulation method satisfying the constraint $\omega p\Delta t\u22652$ will resolve both these plasma oscillations and all other modes of lower frequency.

What is the effect of the speed-limiting on this constraint?

Because of the speed-limiting, the frequency that we must resolve is now not *ω _{p}*; rather, it is the fastest oscillation frequency that the speed-limiting permits. But we have shown in Sec. IV that ordinarily “fast” plasma oscillations [see Eq. (28)] are considerably slowed in SLD. Thus, in its most general form, SLPIC replaces the plasma oscillation constraint by the result [obtained by substituting the frequency derived in Eq. (28) for

*ω*]

_{p}where $\gamma \alpha =v0/2vt\alpha $ is the slowing-down parameter. [For example, when $v0/vte=0.1$, we have $\gamma e=0.07,\u2009h(\gamma e)\u22489.5\xd710\u22123$, and $1/h(\gamma e)\u224810$, thus relaxing the constraint tenfold with minimal effect on the ion modes (assuming IAW phase velocities are low compared to *v*_{0}).] This constraint is less restrictive than the PIC result $\Delta t\u22642/\omega pe$ and permits larger time steps to be taken in SLPIC simulations without instability or loss of accuracy in the low-frequency plasma modes. It can also be shown that by limiting the maximum speed to *v*_{0}, SLPIC trivially relaxes the cell-crossing-time constraint by nearly the same factor (see Appendix B). Both the plasma oscillation constraint and the cell-crossing-time constraint are thus modified by speed-limiting to restrict $\Delta t\u223c\Delta x/v0$.

## V. ANALYSIS OF THE FLUCTUATION SPECTRUM

From the fluctuation–dissipation theorem and the theory of linear response, the fluctuation spectrum of an electrostatic 1D plasma in thermal equilibrium can be shown^{13} to take the form

where *T*_{0} is the equilibrium temperature, $\u27e8E2\u27e9(k)$ is the time-average of the continuous spatial Fourier transform of the square of the electric field, and $D(k,\omega )$ is the dispersion relation. In OD, taking these limits of Eq. (24) yields the result

wherein $\lambda D2\u2261\u03f50T0/\u2211\alpha (q\alpha 2n\alpha )$ is the square of the plasma Debye length.

What is the behavior of the SLD fluctuation spectrum? At low frequencies, where speed-limiting is not expected to influence the wave dynamics, we likewise recover the same result as for OD,

At high frequencies, we have shown that the behavior of solutions to the dispersion relation is significantly altered by the speed-limiting. Nevertheless, terms associated with speed-limiting vanish in the high-frequency limit of the dispersion relation

[because $erfc(\u2212x)\u22612\u2212erfc(x)$]; the effects of speed-limiting should therefore have no bearing on the spatial fluctuation spectrum, Eq. (39). We recover the same result as Eq. (40) for SLD,

In a PIC or SLPIC simulation, the fluctuation spectra of Eqs. (40) and (43) are altered by finite particle size (associated with the transfer of charge and force fields between continuous particle positions and the discrete grid) as well as by finite grid spacing (associated with the wavenumber spectrum that is able to be resolved by the simulation). In addition, the introduction of a discrete grid and a finite volume leads to discrete, finite Fourier spectra and introduces the possibility of aliasing between gridded fields and subgrid particle modes. While we do not propose to discuss the effects of finite particle and grid size in detail in this work, we have used PIC and SLPIC to simulate a 1D single-ion-species hydrogen plasma in thermal equilibrium and have measured its fluctuation spectrum. For a discrete Fourier mode, this spectrum satisfies a relation of the general form^{10,12,13}

wherein *N _{p}* is the number of simulation macroparticles of either species, $|S(k)|2$ is a geometric factor associated with the particle shape, $k=2\pi l/L$ defines the discrete mode index

*l*, $K=k\u2009sinc(k\Delta /2)$ captures the effect of the 1D Laplacian operator on the discrete 1D grid with spacing Δ [with $sinc(x)\u2261\u2009sin\u2009(x)/x$], and $|El|2$ is the time-averaged norm of the

*l*-th discrete mode in the Fourier transform of the electric field.

We modeled this scenario with the VSim^{21} code, using both PIC with a small time step ($\Delta t=1.0\xd710\u221213$ s) and SLPIC with a (50×) larger time step. We used a highly resolved grid (80 cells per Debye length, with a 1D simulation length *L* = 10 Debye lengths), with equilibrium plasma density $n0=5.0\xd71016\u2009m\u22123$ and temperature $T0=10$ eV. One hundred particles per cell of each species were used; the speed limit *v*_{0} for the SLPIC simulations was set to one-half the electron thermal velocity. Particles were mapped to the grid with a three-cell (four-gridpoint) stencil using the method prescribed by Esirkepov.^{22} For this mapping, $|S(k)|2=sinc8(k\Delta /2)$. As shown in Fig. 6, both PIC and SLPIC capture the general behavior of the fluctuation spectrum when the effects of finite grid size and particle width are accounted for. Larger variation between adjacent modes in the spectrum is observed for SLPIC relative to PIC, an effect which becomes more prominent as *γ _{e}* is decreased (stronger speed-limiting) and which is perhaps a function of interparticle correlations imposed by the speed-limiting constraint (since all fast particles now move with the same pseudo-velocity in the domain). Nevertheless, both methods agree closely with the theory. SLPIC simulations with smaller time steps (identical to PIC) were not seen to differ substantially from the large-time step result (green curve) in the figure.

A more detailed consideration of the role of finite grid spacing and particle width in SLPIC is a topic of ongoing interest, and we anticipate future efforts along these lines.

## VI. CONCLUSIONS

In this paper, we have discussed the linear wave dispersion in a 1D1V unmagnetized electrostatic plasma that evolves with both ordinary (OD) and speed-limited (SLD) dynamics. We have demonstrated that speed-limiting can effectively reduce the frequency of fast electron oscillations while quantitatively preserving low-frequency ion and electron motion, e.g., the physics needed to correctly model the Landau damping of ion acoustic waves. We have also shown that this speed-limiting relaxes the “plasma oscillation constraint” of the conventional PIC method, permitting larger time steps, and have demonstrated that the spatial dependence of the ensuing fluctuation spectrum is nevertheless preserved. These findings suggest that the speed-limited particle-in-cell (SLPIC) method, as outlined in the previous work,^{1} is a fast, accurate, and powerful technique for modeling plasmas wherein electron kinetic behavior is significant (such that a fluid/Boltzmann representation for electrons is inadequate), but evolution is on ion timescales. In these cases, the use of PIC is computationally demanding, but the use of speed-limited electrons can substantially reduce computational demands without sacrificing the desired physics.

For plasmas with $vti\u226av\varphi \u226avte$ [where these velocities, respectively, are the ion thermal velocity, ion acoustic wave phase velocity $\omega /k$ ($\u223ccs$ for long-wavelength modes), and the electron thermal velocity], the speed limit *v*_{0} can be chosen with $v\varphi <v0<vte$. Choosing $v0<vte$ ensures that the speed-limited plasma oscillation frequency is reduced below the true plasma oscillation frequency *ω _{p}* by a factor $\omega /\omega p\u223cv0/vte$ and allows the time step to be increased (and the simulation sped up) by a factor of $vte/v0$. However, choosing $v0>v\varphi $ ensures that ion acoustic waves are still accurately simulated (including the Landau damping rate).

Potential applications for the SLPIC method include its use in the modeling of plasma thrusters (wherein very small electron/ion mass ratios impose especially demanding numerical constraints) and sheath formation (e.g., near a Langmuir probe^{23}). Collisional low-temperature plasma discharges are also an area of particular interest; recent efforts have demonstrated that SLPIC can be used in conjunction with standard Monte Carlo collision (MCC) techniques, in the same manner as is done in collisional PIC discharge modeling (PIC-MCC).^{24} We anticipate exploring SLPIC's capability for rapid collisional plasma discharge modeling in future publications.

## ACKNOWLEDGMENTS

This research was financially supported by the U.S. Department of Energy, SBIR Phase I/II Award No. DE-SC0015762, and by the U.S. National Science Foundation, Grant No. PHY1707430. We thank the reviewers of this manuscript for their constructive comments.

## DATA AVAILABILITY

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

### APPENDIX A: ASYMPTOTIC EXPANSIONS OF $Z(\gamma ,\zeta )$

A detailed overview of the properties of the incomplete plasma dispersion function $Z(\gamma ,\zeta )$ was given by Baalrud in Ref. 18. A number of these relations have been used in this work and are summarized here.

For $\zeta \u226b1$,

wherein *H*(*x*) is the Heaviside function, $erfc(x)$ is the complementary error function $erfc(x)=1\u2212erf(x)$, and

For $\zeta \u226a1$,

wherein

is the exponential integral.

Matlab algorithms for evaluating $Z(\zeta ),\u2009Z(\gamma ,\zeta )$, and their derivatives were also provided in the supplemental materials for Ref. 18. These algorithms were used in the numerical calculations of this work.

### APPENDIX B: CONSTRAINTS AND SPEED-LIMITING

In this appendix, we briefly consider the scaling of the numerical constraints outlined in Sec. II in SLD.

The Debye length resolution constraint is independent of the speed-limiting. When it is satisfied, the grid size $\Delta x=\delta \lambda De$ for some $\delta \u22721$.

In SLD, the cell-crossing time constraint is altered by the speed-limiting and it becomes $v0\Delta t<\Delta x$, since no particle can move faster than the speed limit *v*_{0}. Substituting the result from the Debye length resolution constraint then yields the scaling $\Delta t<\delta \lambda De/v0$.

The plasma oscillation constraint, in the limit of aggressive speed-limiting ($\gamma e\u21920$), replaces $\omega p\u223c\omega pe$ by $\omega pev0/vte$ (as shown in the small-*γ* limit of Fig. 1). Substituting the result from the Debye length resolution constraint then yields the scaling $\Delta t<2vte/(v0\omega pe)=2\lambda De/v0$.

It is significant that the time step $\Delta t$ restriction scales linearly as the ratio of Debye length to speed limit in both of the latter two constraints. If this were not so, there could be regions of parameter space where one constraint or the other prevailed, and the speed-limiting concept would be less useful. But the physical scaling for both constraints is the same—the largest time step for speed-limited particles is on the order of the time required for the fastest such particles to cross a Debye length. This restriction preserves the physics of local Debye shielding (time-independent phenomena) even while slowing the rapid plasma oscillations. In addition, this constraint is independent of the electron-ion mass ratio, suggesting that SLPIC can be used even when this ratio is small.

### APPENDIX C: COMPARISON OF SPEED-LIMITED AND RELATIVISTIC DYNAMICS

It has been noted that SLD exhibits some similarities with relativistic dynamics, wherein the speed of light *c* plays a role somewhat analogous to the SLD speed limit *v*_{0}. Although we have not explored this idea in detail in this work, it is instructive to compare the kinetic equation for nonrelativistic SLD [using a more general Lorentz acceleration term that includes the electromagnetic fields $E(x,t)$ and $B(x,t)$] with the relativistic kinetic (Vlasov) equation in the form

Here, respectively, the SLD distribution $f\alpha =f\alpha (x,v,t)$ is a function of position, velocity, and time while the relativistic distribution $g\alpha =g\alpha (x,p,t)$ is a function of position, momentum, and time. The spatial components of the four-velocity, $u\u2261(p/m\alpha )$, are related to the conventional three-velocity **v** through the relativistic Lorentz factor $\gamma =1+u\xb7u/c2$, such that $u=\gamma v$.

The structure of these equations is very similar. The relativistic **u** (whose magnitude may exceed *c*) is like the SLD “true velocity” **v** (whose magnitude may exceed *v*_{0}), and the factor $1/\gamma $ ($\u223c1$ for $|u|\u226ac$ and $\u223cc/|u|$ for $|u|\u226bc$) plays a role akin to the speed-limiting function *β* ($\u223c1$ for $|v|\u226av0$ and $\u223cv0/|v|$ for $|v|\u226bv0$). The product of these relativistic functions, $u/\gamma $ (three-velocity), can never exceed *c* just as the SLD “pseudo-velocity” can never exceed *v*_{0}.

Nevertheless, key differences appear. The term proportional to the electric field, in the relativistic case, contains no physics equivalent to the speed-limiting that occurs in SLD—in effect, relativistic physics applies the speed-limiting concept to the magnetic, but not the electric, components of the Lorentz acceleration. The ensuing trajectories thus vary from those of SLPIC, wherein the appearance of *β* in all but the first term of Eq. (C1) can be viewed as a local rescaling of time (with *β*) along a trajectory that is constant regardless of the value of *v*_{0}. So while SLD is somewhat like relativistic dynamics, in that it tracks both unbounded (SLD true velocity/relativistic momentum) and bounded (SLD pseudo-velocity/relativistic three-velocity) phase space variables, with the latter restricted by fixed speed limits (SLD *v*_{0}/relativistic *c*), the dynamics of the two systems differ enough to make intuitive comparisons difficult.