Dispersion and the Speed-Limited Particle-in-Cell Algorithm
DDispersion and the Speed-Limited Particle-in-Cell Algorithm
Thomas G. Jenkins, a) Gregory R. Werner, and John R. Cary
1, 2 Tech-X Corporation, 5621 Arapahoe Avenue Suite A, Boulder, Colorado 80303,USA Center for Integrated Plasma Studies, University of Colorado, Boulder,Colorado 80309, USA (Dated: 25 February 2021)
This paper discusses temporally continuous and discrete forms of the speed-limitedparticle-in-cell (SLPIC) method first treated by Werner et al. [Phys. Plasmas ,123512 (2018)]. The dispersion relation for a 1D1V electrostatic plasma whose fastparticles are speed-limited is derived and analyzed. By examining the normal modesof this dispersion relation, we show that the imposed speed-limiting substantiallyreduces the frequency of fast electron plasma oscillations while preserving the correctphysics of lower-frequency plasma dynamics (e.g. ion acoustic wave dispersion anddamping). We then demonstrate how the timestep constraints of conventional elec-trostatic particle-in-cell methods are relaxed by the speed-limiting approach, thus en-abling larger timesteps and faster simulations. These results indicate that the SLPICmethod is a fast, accurate, and powerful technique for modeling plasmas whereinelectron kinetic behavior is nontrivial (such that a fluid/Boltzmann representationfor electrons is inadequate) but evolution is on ion timescales.This article has been submitted to Physics of Plasmas. After it is published, it willbe found at https://aip.scitation.org/journal/php. a) [email protected]; https://nucleus.txcorp.com/˜tgjenkins a r X i v : . [ phy s i c s . p l a s m - ph ] F e b . INTRODUCTION The speed-limited particle-in-cell (SLPIC) method is a relatively new plasma model-ing technique. It is most suitable for discharges in which the physics of interest occurs onrelatively slow timescales (e.g., ion transport/profile relaxation) but is nevertheless tied tokinetic electron behaviors that a fluid/Boltzmann model cannot capture (e.g. distributionfunction modifications from neutral collisions or sheath interactions, or Landau damping).In such simulations, one is typically constrained to model both the heavy, slow ion speciesand the light, fast electrons using conventional particle-in-cell (PIC) techniques. The en-suing computational costs can be (possibly unaffordably) high; simulation timesteps mustresolve the electron plasma frequency since kinetic electrons are present, but ion timescalesof 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 thephysics of interest. Larger simulation timesteps can thus be used while retaining the detailedphysics behaviors associated with the slower, longer timescales. The specifics of the SLPICmethod will be explained in a later section of this paper, but we will note here that numericalexperiments using SLPIC simulations to model sheath formation in an argon plasma haveshown that remarkable speedup factors (160 times faster than conventional PIC methods )can be achieved . When SLPIC can be appropriately used for modeling, it is both accurateand powerful.A concept understood since the early days of PIC modeling is that while one may “recovermore of the essence of the situation being simulated by changing the interaction laws”, suchchanges are accompanied by costs: “the more one meddles with the ‘laws’ of nature the moreone must understand the consequences.” While this quote in its original context refers to thevarious approximations used in PIC simulation (e.g. finite-sized particles, grid spacings, andtimesteps), in this paper we explore its relevance to SLPIC. Although SLPIC is in manyways similar to conventional PIC, the ways in which it is different introduce additionaleffects 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, gridspacings, and timesteps (though such effects are not unimportant). Fundamentally, SLPICand PIC methods both seek to statistically approximate the evolution of smooth particle2istribution functions in a multidimensional phase space, in response to self-consistent fieldsand forces – but the underlying evolution equations of the two methods are different and willyield different physics (e.g. linear plasma wave dispersion) even before any particle-basedapproximations are made.In this work, therefore, we focus on developing an understanding of the behavior ofa plasma evolving with speed-limited dynamics (hereafter SLD). In SLD, the plasma isgoverned by continuous, ‘SLPIC-like’ equations of motion that differ from the ones governingplasma evolution in our universe, but which approximate them in certain limits that we willquantify. For purposes of comparison, we also designate the dynamics of plasma evolutionin 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 orPIC, wherein both the timestep ∆ t and the grid spacing ∆ x approach zero as the velocitydependence of the distribution function f α becomes smooth.] We will show that somewell-understood physics processes from OD persist in SLD, and that other processes aresubstantially modified, some of them in very helpful ways.More specifically, in this paper we derive and analyze the analytical dispersion rela-tion in an electrostatic, collisionless, unmagnetized plasma evolving according to modifiedVlasov-Poisson equations which 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-limitingof SLD can quantifiably modify high-frequency behaviors of this plasma (electron plasma os-cillations) while leaving low-frequency motion (ion acoustic wave decay via electron Landaudamping) undisturbed.Section II of this paper explains the SLD concept, together with its connections both tothe SLPIC algorithm and to the conventional kinetic theory of OD. In Section III, we discussthe behavior of an electrostatic ion-electron plasma that evolves with SLD, and derive thedispersion relation associated with this plasma. Section IV contains an analysis of the variouswaves permitted by this dispersion relation, together with the new behaviors imparted bySLD relative to known OD behaviors. We demonstrate that the speed-limiting significantlyrelaxes a fundamental numerical constraint associated with conventional PIC methods and3akes faster numerical simulation possible. Finally, in Section V, we summarize our findings,review additional research directions which these findings might enable, and discuss variousapplications for the SLPIC concept in plasma modeling more generally. II. PHASE SPACE EVOLUTION AND ITS CONNECTION TO PIC ANDSLPIC METHODS
To kinetically model plasma with OD (i.e. using the familiar physics of the real world),the self-consistent evolution of a particle distribution function f α ( x , v , t ) for each species oftype α is sought, according to a phase-space continuity equation ∂∂t f α ( x , v , t ) + ∂∂ x · [ v f α ( x , v , t )] + ∂∂ v · [ a ( x , v , t ) f α ( x , v , t )] = 0 , (1)where a ( x , v , t ) is the acceleration arising from forces acting on the portion of the distributionlocated in the phase space at position x with velocity v at time t . For Hamiltonian systems,the additional phase-space preserving constraint ( ∂/∂ x ) · v + ( ∂/∂ v ) · a = 0 of Liouville’stheorem allows us to rewrite Eq. (1) in the familiar Vlasov form ∂∂t f α ( x , v , t ) + v · ∇ f α ( x , v , t ) + a · ∂∂ v f α ( x , v , t ) = 0 , (2)which can then be solved by the method of characteristics. Along characteristic trajectories d x ( t ) dt = v ( t ) , (3) d v ( t ) dt = a [ x ( t ) , v ( t ) , t ] , (4)the value of the distribution function f α [ x jα ( t ) , v jα ( t ) , t ] is preserved, so that one may rep-resent the distribution function in the form f α ( x , v , t ) = N α (cid:88) j =1 w jα ( t ) δ [ x − x jα ( t )] δ [ v − v jα ( t )] (5)i.e. as a set of N α macroparticles which evolve along the trajectories given by Eqs. (3 – 4).The weight function w jα represents the number of real (physical) particles whose behaviorthe j -th macroparticle encapsulates; in a more general sense it is representative (in somestatistical sense) of the local value of f α in a region near the particle’s initial point on thephase space trajectory. PIC simulation techniques are built upon the fundamental concept4hat Eqs. (3 - 5) solve Eq. (2), and that a sufficiently large number N α of macroparticles, dis-tributed so as to adequately resolve the relevant regions of the phase space, can statisticallyrepresent the 6D+time evolution of the distribution function f α ( x , v , t ).In this paper we will consider Eqs. (2) - (4) as the fundamental equations describingplasma evolution under OD.Techniques for mapping the equations of PIC onto discrete computational timesteps andfinite grids (broadening the spatial extent of macroparticles from delta-functions to small-but-finite widths) are discussed extensively in existing literature ; detailed explanationsand/or derivations of such techniques will not be discussed here except as needed. For thepresent it suffices to note (as the previously cited works discuss) that finite-sized grid cellsand timesteps impose a number of constraints on conventional explicit PIC simulations: • The
Debye length resolution constraint , that a representative grid cell size ∆ x shouldadequately resolve the Debye length λ Dα associated with any of the species in thesimulation in order to avoid numerical heating effects [the Debye length of species α isdefined as λ Dα ≡ (cid:15) T α / ( q α n α ), where { q α , n α , T α } are the species charge, density, andtemperature (in units of energy) and (cid:15) is the permittivity of free space]; • The cell-crossing-time constraint , that the distance traveled by any macroparticle inthe simulation during a finite timestep ∆ t should not exceed a representative grid cellsize ∆ x , so that forces experienced by a particle during a single simulation timestepare adequately resolved; and • The plasma oscillation constraint , that the plasma frequency ω p constrains thetimestep through the relation ω p ∆ t ≤
2; otherwise, numerical instability of thesehigh-frequency oscillations ensues [the plasma frequency is defined as ω p ≡ (cid:80) α ω pα ,with the species plasma frequency ω pα defined through ω pα ≡ q α n α / ( (cid:15) m α ). Here, m α is the mass of species α and the other quantities were defined previously]. Tobe precise, this constraint arises from a more general requirement that every plasmamode frequency must be resolved by the simulation timestep. But for a large classof problems, including the ones considered in this work, the highest mode frequenciesare on the order of ω p .These constraints can impose significant restrictions on a plasma simulation. Low tem-5eratures and/or high densities decrease the Debye length and the allowable grid size, ne-cessitating the use of finer grids and smaller timesteps. Further, when both cold, massiveions and hot, light electrons are simulated with PIC, the timesteps imposed by the plasmaoscillation constraint (now dominated by fast electron motion since ω p ∼ ω pe ) are so smallthat ions may hardly move at all in that time interval. For simulations where many periodsof harmonic ion motion are of interest, the number of timesteps required can be enormous.The speed-limited particle-in-cell approach, and the more general speed-limiting conceptswe explore in this work, are motivated by a desire to relax some of these constraints. The keyidea is simple: the fastest particles and highest-frequency wave phenomena necessitate thesmallest timesteps, and if these fast particle and wave motions can be slowed, the timestepconstraints can be relaxed.Accordingly, we introduce speed-limited dynamics (SLD), wherein equations from thederivation of the SLPIC method presented in Ref. 1 govern the plasma dynamics. Here, adistribution function f α of species α evolves as prescribed by a modified Vlasov equation ∂∂t f α ( x , v , t ) + β ( v ) v · ∂∂ x f α ( x , v , t ) + β ( v ) a · ∂∂ v f α ( x , v , t ) = 0 , (6)wherein a speed-limiting function β ( v ) in the range (0 ,
1] has been introduced. This functiontransitions from values at or near unity (for “slow” particles) to values approaching v / | v | (for “fast” particles), and we can understand its effect by looking at the characteristictrajectories of the modified Vlasov equation d x ( t ) dt = β [ v ( t )] v ( t ) , (7) d v ( t ) dt = β [ v ( t )] a [ x ( t ) , v ( t ) , t ] . (8)The product | v | β ( 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 ∼ | v | at low velocities but is limited tovalue v (in the original vector direction of motion) at high velocities. We will hereafter referto v as the “speed limit”; it is the upper bound on the rate at which motion in the positioncoordinate of phase space may proceed. Accordingly, some nuance is required in discussingthe phase space evolution since the meaning of ‘velocity’ is now ambiguous. An elementof phase space has both a ‘true velocity’ (the phase space coordinate v ) and a ‘pseudo-velocity’ d x /dt = β v (the speed and direction at which it is permitted to move from onephysical space coordinate to another). In a given time interval, elements in the phase space6ith large true velocity v ( t ) [so that β ( v ) (cid:28)
1] experience both smaller pseudo-velocities[Eq. (7)] as they move through the space, and smaller changes to these pseudo-velocities(pseudo-acceleration) in response to applied forces [Eq. (8)]. Elements in the phase spacewith small true velocity experience no speed-limiting [ β ( v ) ∼
1] and evolve in the samemanner as their OD counterparts, as in Eqs. (3) and (4).In this paper we will consider Eqs. (6) - (8) as the fundamental equations describingplasma evolution under SLD.Solutions to the SLD Vlasov equation, Eq. (6), can be represented statistically in thesame manner as outlined above, setting f α ( x , v , t ) = N α (cid:88) j =1 w jα ( t ) δ [ x − x jα ( t )] δ [ v − v jα ( t )] (9)for a suitably large number N α of macroparticles evolving along the characteristic trajecto-ries described by Eqs. (7) – (8). This is the speed-limited particle-in-cell method we havepresented in previous work . However, this work will not focus on PIC or SLPIC implemen-tations of Eqs. (2) or (6). 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 equationof SLD to model dispersion in an electrostatic, unmagnetized plasma with a single ionspecies, in one spatial dimension and one velocity-space dimension. We will use the analyticforms of these equations (foregoing for the moment any discussion of the effects of finitetimesteps, grid spacings, or particle sizes) to ensure that we understand the new physicsthat the imposed speed-limiting of SLD imparts. Each species will use the Vlasov equation ∂f α ( x, v, t ) ∂t + β ( v ) v ∂f α ( x, v, t ) ∂x − β ( v ) q α m α ∂φ ( x, t ) ∂x ∂f α ( x, v, t ) ∂v = 0 . (10)In SLD, we will use a speed-limiting function of form β ( v ) = − v v + (cid:16) v v + 1 (cid:17) H ( v + v ) + (cid:16) v v − (cid:17) H ( v − v ) , (11)wherein H ( x ) is the Heaviside function. An equivalent representation for this speed-limitingfunction is β ( v ) = | v | ≤ v v / | v | ; | v | ≥ v . (12)7n OD, β ( v ) = 1 (the v → ∞ limit of the SLD). The species couple via the Poisson equation, ∂ φ ( x, t ) ∂x = − (cid:88) α q α (cid:15) (cid:90) ∞−∞ f α ( x, v, t ) dv . (13)Many equilibrium solutions of Eq. (10) are possible. We will choose physically reasonablesolutions that are stationary Maxwellians in each of the individual species (though we willallow the two species to have different temperatures, and neglect both the collisional pro-cesses that have brought the individual species to their present state and the interspeciescollision processes that would further relax the system to a single temperature). Formally,we write the distribution function of species α as f α ( x, v, t ) = n (cid:114) m α πT α exp (cid:18) − m α v T α (cid:19) , (14)where n is a species-independent constant number density, T α is the constant tempera-ture of species α , and m α is the species mass. With these equilibrium distributions, thecorresponding equilibrium potential φ ( x, t ) is a constant that can be set to zero.Linearizing Eqs. (10 – 11) in perturbed quantities, and Fourier transforming from space-time coordinates { x, t } to wavenumber and frequency coordinates { k, ω } , yields the result − iωf α ( k, v, ω ) + ikβ ( v ) vf α ( k, v, ω ) = β ( v ) q α m e ikφ ( k, ω ) ∂f α ( v ) ∂v , (15) − k φ ( k, ω ) = − (cid:88) α q α (cid:15) (cid:90) ∞−∞ f α ( k, v, ω ) dv . (16)We obtain, for the perturbed distribution functions, f α ( k, v, ω ) = kvβ ( v ) ω − kvβ ( v ) q α φ ( k, ω ) T α n (cid:114) m α πT α exp (cid:18) − m α v T α (cid:19) (17)which we can then substitute into Eq. (16). The ensuing integrals will be undefined for reso-nant velocities vβ ( v ) = ω/k ; we will implicitly stipulate that the integrals are to be evaluatedusing the Landau contour (traversing below any singularity) to retain both causality and theresonant physics. (Formally, this can be shown to be equivalent to the use of a Laplace trans-form, rather than a Fourier transform, in the time domain; it also permits the generalizationof ω to complex values.) We obtain the integral relation1 − (cid:88) α q α n k (cid:15) T α (cid:90) ∞−∞ kvβ ( v ) ω − kvβ ( v ) (cid:114) m α πT α e − m α v / (2 T α ) dv = 0 (18)8hich describes the dispersive wave behavior of the plasma in both OD [where β ( v ) = 1]and SLD [where Eq. (12) defines β ( v )].In OD, the integral in Eq. (18) can be expressed in terms of the plasma dispersionfunction , defined as Z ( ζ ) = 1 √ π (cid:90) ∞−∞ e − t t − ζ dt (19)for Im( ζ ) > ζ ) ≤
0. In SLD, this integral is morecomplicated; given our choice for β ( v ), there are both high-velocity regions of integrationwherein vβ ( v ) = ± v , a constant (the speed limit), and low-velocity regions wherein vβ ( v ) = v . We may represent the integrals over these various regions in terms of the complementaryerror function and the incomplete plasma dispersion function. The latter function wasintroduced by Franklin and its properties have been discussed extensively by Baalrud ; ittakes the form [generalized from Eq. (19)] Z ( γ, ζ ) = 1 √ π (cid:90) ∞ γ e − t t − ζ dt (20)for Im( ζ ) > ζ ) ≤
0. It will be useful to note that Z ( −∞ , ζ ) = Z ( ζ ) and that Z ( ∞ , ζ ) = 0. Additional properties of this function are providedin Appendix A.In terms of these functions, we may rewrite Eq. (18) in the form D SLD ( k, ω ; v ) = 0,where we define D SLD ( k, ω ; v ) ≡ (cid:88) α k λ Dα (cid:20) ζ α erfc( γ α ) γ α − ζ α + ζ α [ Z ( − γ α , ζ α ) − Z ( γ α , ζ α )] (cid:21) . (21)Here, we make use of the parameters γ α ≡ v / ( √ v tα ), a measure of the relative speed-limiting of species α ; ζ α = ω/ ( √ kv tα ), the conventional (and complex) argument of theplasma dispersion function; λ Dα ≡ (cid:15) T α / ( q i n ), the Debye length of species α ; and v tα ≡ T α /m α , the thermal velocity for species α . The complementary error function is related tothe conventional error function by erfc( x ) ≡ − erf( x ). The expression D SLD ( k, ω ; v ) = 0is the plasma dispersion relation of SLD, and its analysis and solutions will be the topic ofthe following section.Recalling that erfc( ∞ ) = 0 and the limits of the incomplete plasma dispersion functiondiscussed above, we can show that D SLD ( k, ω ; ∞ ) = D OD ( k, ω ), where the latter function9as the explicit form D OD ( k, ω ) = 1 + (cid:88) α k λ Dα [1 + ζ α Z ( ζ α )] , (22)a standard result from elementary plasma kinetic theory . This is the plasma dispersionrelation of OD, against which solutions with SLD will be compared. IV. ANALYSIS OF THE DISPERSION RELATION
We now consider various limits of D OD ( k, ω ) and D SLD ( k, ω ; v ), together with the wavesthat arise in these various limits. For clarity, we will consider the physics of OD first, andwill then explore the changes which the speed-limiting of SLD imparts to these familiarprocesses. A. Plasma oscillations
The OD dispersion relation, Eq. (22), admits approximate analytic solutions correspond-ing to cold-plasma oscillations in the ζ α (cid:29) Z ( ζ α ) inthis limit, Z ( ζ α ) ∼ − ζ α (cid:18) ζ α + 34 ζ α + . . . (cid:19) + O (cid:16) e − ζ α (cid:17) (23)can be substituted into Eq. (22) to obtain, at lowest order, D OD ( k, ω ) = 1 − ω pe ω − ω pi ω = 0 (24)where ω pα = q i n / ( (cid:15) m α ) is the square of the plasma frequency of species α . These oscil-lations are dominated by electron motion (the ion term is order m e /m i smaller than theelectron term); we may write the solution as ω = ω pe (1 + m e /m i ) ≡ ω p . We anticipatethe prospect of significant changes to these oscillations in SLD, since the speed-limiting willpreferentially modify the fast electron motion.What is the behavior of the SLD dispersion relation, Eq. (21), in the same ζ α (cid:29) Z ( γ, ζ ) presented in Ref. 9 and summarized inAppendix A, we can show that the SLD version of Eq. (24) takes the form D SLD ( k, ω ; v ) = 1 − (cid:88) α ω pα ω (cid:32) γ α − γ α ) − γ α e − γ α √ π (cid:33) = 0 (25)10hich reproduces Eq. (24) in the v → ∞ (i.e. γ α → ∞ ) limit. It admits the solutions ω = (cid:88) α ω pα h ( γ α ) ; h ( γ α ) ≡ (cid:32) γ α − γ α ) − γ α e − γ α √ π (cid:33) . (26)The behavior of the function h ( γ α ) is shown in Figure 1. For small γ α , h ( γ α ) ∼ γ α , whilefor large γ α , the species is not appreciably speed-limited and h ( γ α ) ≈
1. Recalling that γ α = v / √ v tα , it will be instructive to express the argument of h strictly in terms of theelectron slowing-down parameter γ e ; we have γ i = γ e (cid:112) T e m i / ( T i m e ). The presence of theion-electron mass ratio suggests that unless electrons are much cooler than ions, the quantity γ i will always be much larger than γ e . Accordingly, we may generally choose v in a way thatalters electron behavior but not ion behavior, consistent with the intentional slowing downof the fastest particles in SLD while leaving slower particles undisturbed. Such a v will befaster than most ions but much slower than most electrons. For such a choice ( v = 4 v ti ),Figure 1 also shows values of γ e , γ i , and their corresponding h values for a hydrogen plasmathermalized to 10 eV. In this plasma, ion motion associated with plasma oscillations doesnot 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 frequenciespredicted by Eq. (26), we examine the behavior of these modified electron plasma oscillationsas a function of the ratio v /v te (providing intuition as to which velocities in a typical elec-tron distribution, e.g. a Maxwellian, are restricted by the speed-limiting). The normalizedoscillation frequency of Eq. (26) is shown in Fig. (2) for a monatomic hydrogen plasma withequal electron and ion temperatures. While the oscillation frequencies of SLD and OD areidentical at large values of v /v te (minimal speed-limiting), reducing this ratio, and henceincreasing the corresponding fraction erfc( γ e ) of speed-limited electrons, reduces the SLDfrequency monotonically. As v /v te →
0, the linear approximation to the plasma frequencyapproaches the heuristic estimate made in Ref. 1: ω/ω p ∼ v /v te . The high-frequency oscil-lations of OD have been mapped to lower-frequency oscillations in SLD by the speed-limitingof electrons.We conclude that in SLD, the frequency of conventional electron-dominated plasma os-cillations is reduced relative to OD. In Sec. IV C we will consider how this reduction relaxesthe ‘plasma oscillation constraint’ referred to in Sec. II. But we must first determine whetherthe speed-limited dynamics preserves physics associated with lower-frequency waves.11 ○◇◇ γ h ( γ ) ○ ions ◇ electrons FIG. 1. Behavior of the function h ( γ α ) in Eq. (26) as a function of the species speed-limitingparameter γ α = v / ( √ v tα ). For high values of v relative to the species thermal velocity, speed-limiting does not occur in densely populated portions of the phase space and this function (amultiplicative factor in the dispersion relation) approaches unity. When v ti (cid:28) v (cid:28) v te , consider-able speed-limiting of the electron distribution can be achieved without appreciable effect on theion distribution, such that h ( γ e ) (cid:28) h ( γ i ) ≈
1. To illustrate this point, values of γ α andthe ensuing h ( γ α ) are shown for a case with v = 4 v ti , in a hydrogen plasma with both species attemperature 10 eV. B. Ion acoustic waves
The OD dispersion relation, Eq. (22), also admits approximate analytic solutions corre-sponding to ion acoustic waves. Physically, these solutions are associated with ‘hot’ electrons( ζ e (cid:28)
1) and ‘cold’ ions ( ζ i (cid:29) k λ De D OD ( k, ω ) = 1 + k λ De + ζ e i √ π − k c s ω = 0 , (27)where the multiplicative factor of k λ De simplifies the algebra and where c s , the sound speed,satisfies the relation c s = T e /m i . Assuming that the imaginary part of the (now complex)12 /v te / p SLDODv = v te FIG. 2. Frequency of electron-dominated cold plasma oscillations in SLD and OD as a functionof the ratio of speed limit to electron thermal speed ( v /v te = √ γ e ). The electron-dominatedcold plasma oscillation frequency is constant in OD (red dashed line), but in SLD it unphysicallydecreases (blue solid line) as v /v te is reduced. This unphysical frequency means that plasmaoscillations are not correctly simulated, but this can greatly speed up simulation when plasmaoscillations are unimportant to the physics of interest. When v (cid:28) v te , ω/ω p ∼ v /v te . ω = ω r + iω i is small, we can perform a Taylor expansion, D ( k, ω ) ≈ D ( k, ω r ) + iω i ∂D ( k, ω ) ∂ω (cid:12)(cid:12)(cid:12)(cid:12) ω = ω r (28)to show that this equation permits damped wavelike solutions of the form ω = ± kc s (cid:112) k λ De − i (cid:114) (cid:15)π kc s (1 + k λ De ) (29)13here (cid:15) 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 Landaudamping.What is the corresponding behavior of the SLD dispersion relation? In the ζ i (cid:29) , ζ e (cid:28) k λ De D SLD ( k, ω ; v ) = 1 + k λ De − k c s ω h ( γ i ) + i √ πζ e [ H ( ζ e − γ e ) − H ( ζ e + γ e )] = 0 . (30)Repeating the Taylor expansion procedure above, we can show that the ion acoustic wavedispersion relation of SLD has approximate analytic solutions ω = ± ω r − i (cid:114) (cid:15)π kc s h ( γ i )(1 + k λ De ) (cid:104) H (cid:16) ω r k + v (cid:17) − H (cid:16) ω r k − v (cid:17)(cid:105) , (31)wherein ω r ≡ kc s (cid:112) h ( γ i ) (cid:112) k λ De . (32)This is the SLD equivalent to Eq. (29). The OD and SLD forms are approximately equivalentprovided that γ i is sufficiently large (so that h ( γ i ) ∼
1, see Fig. 1) and | ω r /k | < v (sothat the speed limit v exceeds the phase velocity of the ion acoustic wave). Since theapproximation h ( γ i ) ≈ or better for v > v ti , we will write theseconstraints as a single condition v > max(4 v ti , | ω r /k | ) . (33)When v is chosen to satisfy this condition, and when the other conditions we assumed in thederivation [ ζ i (cid:29) , ζ e (cid:28) , Im( ω ) (cid:28) Re( ω )] are valid, the effects of SLD on the propagationand damping of the IAW are minimal.What happens when the phase velocity condition is violated, such that v < | ω r /k | ? Inthis case, the explicitly imaginary terms in Eq. (30) (proportional to Heaviside functions)vanish, and the solutions admitted now only capture the real part of Eq. (31). While the IAWstill propagates, its Landau damping is not correctly modeled. This is consistent with a resultthat we have demonstrated in previous work , namely, that the correct dynamics of resonantwave-particle interactions cannot be captured by speed-limited particle-in-cell simulationswhen particles whose velocities were previously synchronous with the wave velocity arespeed-limited. Particles whose speed-limiting renders them unable to keep up with the wave14annot exchange energy with it, so the dissipative effects which lead to wave damping areeffectively turned off. Although nothing in principle prevents us from choosing a smaller v value that violates Eq. (33), the inherent advantage of the velocity-dependent speed-limitingapproach (preservation of IAW physics) would be lost in doing so.The analytic form of the IAW above is approximate. The dispersion relation can alsobe solved numerically to find the complex ω associated with a given k , and we have doneso for a plasma with density n = 5 . × m − , T e = 10 eV, and T i = 1 /
40 eV. Inthe supplemental materials for Ref. 9, Matlab algorithms for evaluating Z ( ζ ), Z ( γ, ζ ), andtheir derivatives have been provided. We have made use of these algorithms and Matlabminimization routines to compute { ω, k } values which satisfy the OD and SLD dispersionrelations within a given tolerance, according to norm-minimization criteria (cid:112) | D SLD ( k, ω ; v ) | ≤ δ SLD (34) (cid:112) | D OD ( k, ω ) | ≤ δ OD , (35)for each of the dispersion relations. In these computations we fix k (and, for SLD, v ), andthen numerically evaluate the function D for various complex values of ω . Exact solutionsto the dispersion relation have D = 0; we vary ω in the complex plane in a manner thatseeks to minimize the norm of D and thus to approach these exact solutions.Various values for the speed-limiting parameter v can be chosen to assess the effect ofspeed-limiting on IAW behavior. Depending on the value of v and k , to obtain numericalconvergence it is sometimes necessary to raise δ SLD to values as high as 10 − , but solutionscan usually be found for δ OD < − and δ SLD < − (and often for δ values that are severalorders of magnitude smaller, when k (cid:29) ω for various values of wavenumber k spanning several orders of magnitude. In this figure,the speed limit v has the value 0 . v te ; the analytic approximation to the dispersion relation(valid for kλ De (cid:46)
1) is also shown. The speed-limiting does not affect the IAW behaviorappreciably; 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 percentlow relative to the OD value and the SLD damping rate also drops by about 3%. Forshorter-wavelength modes, the effect of the speed-limiting on both the real frequencies anddamping rates is negligible. Nevertheless, the frequency of the modified plasma oscillations15 -3 -2 -1 k De -4 -2 / p i Re( )/ pi , OD-Im( )/ pi , ODAnalytic Approx. ODRe( )/ pi , v = 0.2 v te -Im( )/ pi , v = 0.2 v te FIG. 3. Real and negative-imaginary parts of the normalized IAW frequency, computed nu-merically from the OD (red) and SLD (blue) dispersion relations, as a function of normalizedwavenumber. Values from the approximate analytic expression, Eq. (29), are also shown (dashedblack curves) and agree well with the exact solutions when kλ De (cid:46)
1. In the SLD case the speedlimit v = 0 . v te . Even though all electrons are restricted to move through the domain with speedsno greater than v , the frequency (good to within 0.5%) and damping rate (good to within 3%) forthe ion acoustic wave do not differ appreciably from the OD values. Good agreement between SLDand OD is maintained even in the kλ De > [from Fig. (2)] is decreased to about 0 . ω pe for this case.For more restrictive speed limits v , greater deviation of IAW frequencies and dampingrates from their non-speed limited (OD) values is observed, though still only for modes whosewavelengths are long compared to the Debye length. Generally speaking, real frequencies16 k v /Re( ) 10 -1 R e () / p i k De =100k De =50k De =20k De =10k De =5k De =2k De =1k De =0.5k De =0.2k De =0.1OD values FIG. 4. Variation of the normalized real part of the frequencies which solve the ion acoustic branchof the SLD dispersion relation, as a function of the speed limit v normalized to the wave phasevelocity v φ = Re( ω ) /k , for various values of kλ De . The OD frequencies, which are independent of v , are also shown for each kλ De value (red dashed lines). Although speed-limiting does not affectthe mode frequencies for v /v φ (cid:29)
1, modes whose wavelengths are large compared to the Debyelength ( kλ De (cid:46) .
0) are reduced in frequency as the speed limit is reduced (moving to the left onthe graph) to be of the same order as the phase velocity. This frequency variation is minimal for v > v φ and also for modes whose wavelengths are short relative to the Debye length. For thiscase v ti /v te = 0 . (cid:112) m e /m i . are shifted downward as the speed limit is decreased, though only by a few percent ( < v = 0 . v te = 2 . c s , and <
2% for v = 0 . v te = 4 . c s ). Damping rates also (generally)decrease in magnitude as the speed limit is decreased, but the magnitude of the relativedecrease is larger ( <
42% for v = 0 . v te = 2 . c s , <
13% for v = 0 . v te = 4 . c s ). In Figures17 k v /Re( ) 10 -6 -4 -2 - I m () / p i k De =100k De =50k De =20k De =10k De =5k De =2k De =1k De =0.5k De =0.2k De =0.1OD values FIG. 5. Variation of the normalized negative-imaginary part (damping rate) of the frequencieswhich solve the ion acoustic branch of the SLD dispersion relation, as a function of the speed-limit v normalized to to the wave phase velocity v φ = Re( ω ) /k , for various values of kλ De . The ODdamping rates, which are independent of v , are also shown for each kλ De value (red dashed lines).Although speed-limiting does not affect the damping rates when v /v φ (cid:29)
1, damping rates formodes whose wavelengths are large compared to the Debye length ( kλ De (cid:46) .
0) may be reducedby several orders of magnitude as the speed limit is reduced (moving to the left on the graph) tobe of the same order as the phase velocity. As is also true for the real frequencies (Fig. 4), thedamping rate variation is minimal for v > v φ and also for modes whose wavelengths are shortrelative to the Debye length. For this case v ti /v te = 0 . (cid:112) m e /m i . kλ De (cid:46)
1) in SLD, speed limitsthat approach the wave phase velocity generally give rise to modest reductions in the modefrequency (Fig. 4) and more pronounced reductions of the IAW damping rate (Fig. 5). Thisfailure to properly capture the Landau damping of the IAW arises because portions of thedistribution function that resonate with the wave (in OD) are prevented from doing so bythe imposed speed-limiting . As we transition to shorter-wavelength modes (moving up thegraph legend), absolute damping rates are increased to become comparable to the real modefrequency, 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( v ∼ v ti , so as not to speed-limit the bulk ion distribution) still exceeds v φ = ω/k for large k . Accordingly, we may assert that SLD preserves the physics of IAW propagation anddamping provided that the condition in Eq. (33) holds, namely, when ions are not speed-limited and when the speed limit reasonably exceeds the phase velocity of the IAWs in thesystem. At the same time, this speed limiting considerably reduces the frequency of plasmaoscillations (as was shown in Section IV A). C. Normal modes and the plasma oscillation constraint
Having explored the normal modes of our 1D1V plasma, we now revisit the numericalconstraints which the use of a finite timestep ∆ t will impose on particle-in-cell simulationsof this plasma.We have noted in Section II that numerical instability will ensue in a PIC simulation ifthe frequency of any plasma mode is not resolved, and that in particular, we must resolve thefrequency associated with electron plasma oscillations. These oscillations are generally thehighest-frequency modes in a PIC simulation containing electrons (because ions introduceonly small corrections to the approximation ω p ≈ ω pe when the electron-ion mass ratio issmall). Accordingly, any numerical simulation method satisfying the constraint ω p ∆ t ≥ ω p ; rather,it is the fastest oscillation frequency that the speed-limiting permits. But we have shownin Sec. IV that ordinarily ‘fast’ plasma oscillations [see Eq. (26)] are considerably slowed inSLD. Thus, in its most general form, SLPIC replaces the plasma oscillation constraint bythe result [obtained by substituting the frequency derived in Eq. (26) for ω p ]∆ t ≤ (cid:113) ω pe h ( γ e ) + ω pi h ( γ i ) ≈ ω pe (cid:112) h ( γ e ) (36)where γ α = v / √ v tα is the slowing-down parameter. [For example, when v /v te = 0 .
1, wehave γ e = 0 . h ( γ e ) ≈ . × − , and 1 / (cid:112) h ( γ e ) ≈
10, thus relaxing the constraint tenfoldwith minimal effect on the ion modes (assuming IAW phase velocities are low compared to v ).] This constraint is less restrictive than the PIC result ∆ t ≤ /ω pe , and permits largertimesteps to be taken in SLPIC simulations without instability or loss of accuracy in thelow-frequency plasma modes. It can also be shown that by limiting the maximum speed to v , SLPIC trivially relaxes the cell-crossing-time constraint by nearly the same factor (seeAppendix B). Both the plasma oscillation constraint and the cell-crossing-time constraintare thus modified by speed-limiting to restrict ∆ t ∼ ∆ x/v . V. CONCLUSIONS
In this paper we have discussed the linear wave dispersion in a 1D1V unmagnetizedelectrostatic plasma that evolves with both ordinary (OD) and speed-limited (SLD) dynam-ics. We have demonstrated that speed-limiting can effectively reduce the frequency of fastelectron 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. Wehave also shown that this speed-limiting relaxes the “plasma oscillation constraint” of theconventional PIC method, permitting larger timesteps. These findings suggest that thespeed-limited particle-in-cell (SLPIC) method, as outlined in previous work , is a fast, ac-curate, and powerful technique for modeling plasmas wherein electron kinetic behavior issignificant (such that a fluid/Boltzmann representation for electrons is inadequate) but evo-lution is on ion timescales. In these cases the use of PIC is computationally demanding, butthe use of speed-limited electrons can substantially reduce computational demands withoutsacrificing the desired physics. 20or plasmas with v ti (cid:28) v φ (cid:28) v te [where these velocities respectively are the ion thermalvelocity, ion acoustic wave phase velocity ω/k ( ∼ c s for long-wavelength modes), and theelectron thermal velocity], the speed limit v can be chosen with v φ < v < v te . Choosing v < v te ensures that the speed-limited plasma oscillation frequency is reduced below thetrue plasma oscillation frequency ω p by a factor ω/ω p ∼ v /v te , and allows the timestep tobe increased (and the simulation sped up) by a factor v te /v . However, choosing v > v φ ensures that ion acoustic waves are still accurately simulated (including the Landau dampingrate).Potential applications for the SLPIC method include its use in the modeling of plasmathrusters (wherein very small electron/ion mass ratios impose especially demanding nu-merical constraints) and sheath formation (e.g. near a Langmuir probe ). Collisional low-temperature plasma discharges are also an area of particular interest; recent efforts havedemonstrated that SLPIC can be used in conjunction with standard Monte Carlo collisiontechniques, in the same manner as is done in collisional PIC discharge modeling (PIC-MCC) . We anticipate exploring SLPIC’s capability for rapid collisional plasma dischargemodeling in future publications. ACKNOWLEDGMENTS
This research was financially supported by the U.S. Department of Energy, SBIRPhase I/II Award DE-SC0015762, and by the U.S. National Science Foundation, GrantPHY1707430. The data that support the findings of this work (numerical solutions to thedispersion relations) are available from the corresponding author upon reasonable request.
Appendix A: Asymptotic expansions of Z ( γ, ζ )A detailed overview of the properties of the incomplete plasma dispersion function Z ( γ, ζ )was given by Baalrud in Ref. 9. A number of these relations have been used in this workand are summarized here.For ζ (cid:29) Z ( γ, ζ ) ∼ iσ √ πH ( ζ − γ ) e − ζ − erfc( γ )2 ζ − e − γ √ πζ − ζ (cid:32) γe − γ √ π + erfc( γ )4 (cid:33) − . . . (A1)21herein H ( x ) is the Heaviside function, erfc( x ) is the complementary error function erfc( x ) =1 − erf( x ), and σ = 1 − sign[Im( ζ )] . (A2)For ζ (cid:28) Z ( γ, ζ ) ∼ i √ πH ( ζ − γ ) e − ζ + E ( γ )2 √ π + ζ (cid:32) e − γ γ √ π − erfc( γ ) (cid:33) + ζ (cid:32) e − γ √ πγ − E ( γ )2 √ π (cid:33) + . . . (A3)wherein E ( x ) = (cid:90) ∞ e − xt t dt (A4)is the exponential integral.Matlab algorithms for evaluating Z ( ζ ), Z ( γ, ζ ), and their derivatives were also providedin the supplemental materials for Ref. 9. These algorithms were used in the numericalcalculations of this work. Appendix B: Constraints and speed-limiting
In this appendix we briefly consider the scaling of the numerical constraints outlined inSection II in SLD.The Debye length resolution constraint is independent of the speed-limiting. When it issatisfied, the grid size ∆ x = δλ De for some δ (cid:46) v ∆ t < ∆ x , since no particle can move faster than the speed limit v . Substituting theresult from the Debye length resolution constraint then yields the scaling ∆ t < δλ De /v .The plasma oscillation constraint, in the limit of aggressive speed-limiting ( γ e → ω p ∼ ω pe by ω pe v /v te (as shown in the small- γ limit of Fig. 1). Substituting the re-sult from the Debye length resolution constraint then yields the scaling ∆ t < v te / ( v ω pe ) =2 λ De /v .It is significant that the timestep ∆ t restriction scales linearly as the ratio of Debye lengthto speed limit in both of the latter two constraints. If this were not so, there could be regionsof parameter space where one constraint or the other prevailed, and the speed-limitingconcept would be less useful. But the physical scaling for both constraints is the same – thelargest timestep for speed-limited particles is on the order of the time required for the fastest22uch particles to cross a Debye length. This restriction preserves the physics of local Debyeshielding (a time-independent phenomena) even while slowing the rapid plasma oscillations.In addition, this constraint is independent of the electron-ion mass ratio, suggesting thatSLPIC can be used even when this ratio is small. REFERENCES G. R. Werner, T. G. Jenkins, A. M. Chap, and J. R. Cary, Phys. Plasmas , 123512(2018). Subsequent code development has enabled speedup factors of greater than 250 for thisdischarge, relative to conventional PIC. Detailed discharge properties are provided in Ref.1. A. B. Langdon and C. K. Birdsall, Phys. Fluids , 2115 (1970). Collisional effects, which would replace the zero on the right-hand side of Eq. (1) withsource or sink terms, could also be included in this equation; SLPIC is compatible withconventional PIC-MCC techniques for modeling collisional plasmas. We will not considercollisional effects in this work, but future publications demonstrating collisional SLPICdischarges are anticipated. C. K. Birdsall and A. B. Langdon,
Plasma Physics via Computer Simulation , CRC Press,2004. R. W. Hockney and J. W. Eastwood,
Computer Simulation Using Particles , CRC Press,1988. B. D. Fried and S. D. Conte,
The Plasma Dispersion Function , Academic Press, 1961. R. N. Franklin, in
Proceedings of the Tenth International Conference on Phenomena inIonized Gases , page 269, Donald Parsons and Company, Ltd., 1971. S. D. Baalrud, Phys. Plasmas , 012118 (2013). R. J. Goldston and P. H. Rutherford,
Introduction To Plasma Physics , IoP, 1995. F. F. Chen,
Introduction to Plasma Physics and Controlled Fusion , Springer, 3rd edition,2016. G. R. Werner, S. Robertson, T. G. Jenkins, A. M. Chap, and J. R. Cary, “AcceleratedSteady-State Electrostatic Particle-in-Cell Simulation of Langmuir Probes”, to be submit-ted to Phys. Plasmas. 233