The process of numerical thermalization in particle-in-cell (PIC) simulations has been studied extensively. It is analogous to Coulomb collisions in real plasmas, causing particle velocity distributions (VDFs) to evolve toward a Maxwellian as macroparticles experience polarization drag and resonantly interact with the fluctuation spectrum. This paper presents a practical tutorial on the effects of numerical thermalization in 2D PIC applications. Scenarios of interest include simulations, which must be run for many thousands of plasma periods and contain a population of cold electrons that leave the simulation space very slowly. This is particularly relevant to many low-temperature plasma discharges and materials processing applications. We present numerical drag and diffusion coefficients and their associated timescales for a variety of grid resolutions, discussing the circumstances under which the electron VDF is modified by numerical thermalization. Though the effects described here have been known for many decades, direct comparison of analytically derived, velocity-dependent numerical relaxation timescales to those of other relevant processes has not often been applied in practice due to complications that arise in calculating thermalization rates in 1D simulations. Using these comparisons, we estimate the impact of numerical thermalization in several examples of low-temperature plasma applications including capacitively coupled plasma discharges, inductively coupled plasma discharges, beam plasmas, and hollow cathode discharges. Finally, we discuss possible strategies for mitigating numerical relaxation effects in 2D PIC simulations.

The particle-in-cell (PIC) method of plasma simulation is frequently used for simulations of plasma discharges in plasma processing.1–5 The importance of avoiding artificial numerical effects caused by the increased density fluctuations inherent in the PIC method is often insufficiently analyzed. The main focus of this paper is to provide accurate theoretical estimates of the effect of numerical noise on the electron velocity distribution function (EVDF) and give examples of such analysis for widely used discharges in plasma processing: capacitively and inductively coupled plasmas (ICPs), hollow cathode plasmas, and electron beam-generated plasmas.

In most PIC simulations of low-temperature plasmas, the constraint on the number of macroparticles required to achieve accuracy using the particle-in-cell (PIC) method of plasma simulation is often expressed more qualitatively than quantitatively, prescribing that the number of macroparticles per cell should be “large” (several hundred or several thousand). Convergence tests are recommended to ensure that the results of a simulation do not depend on the number of macroparticles used. While such convergence tests are important, they should be aided by a more rigorous analysis of the nature of the error caused by numerical noise. The theory describing the effect of numerical noise on the EVDF has been developed decades ago. However, it is not often used due to the so-called “kinetic blocking” effects that reduce the effect of numerical noise in 1D PIC simulations.6,7 The rate of numerical thermalization is often significantly slower in 1D than in 2D, and analytical prediction of the relaxation time is difficult in 1D. By contrast, the theory is relatively simple in multidimensional PIC simulations and the thermalization rates are often much faster.8,9 Until recently, most PIC simulations were 1D, and hence, the effects of noise were not of great concern. The situation is drastically different in 2D where, as we show below, the effects of noise and the corresponding thermalization can strongly affect the EVDF. Hence, these effects need to be carefully evaluated.

The results presented here are a pedagogical synthesis of the available information for practical application to 2D PIC simulations. We are particularly focused on PIC simulations of low-temperature plasma discharges, where fluid models are frequently used under the assumption that particle velocity distribution functions (VDFs) are already Maxwellian. The use of a kinetic simulation technique such as the PIC method implies that there are possible kinetic effects to be investigated. However, without adequate attention to numerical thermalization errors, the potential features of interest in the electron velocity distribution function (EVDF)—such as an excess or depleted population of high-velocity electrons—will not be accurately represented.

Numerical thermalization drives particle VDFs toward a Maxwellian, with EVDFs thermalizing more rapidly than ion velocity distribution functions (IVDFs). This thermal relaxation often happens on a timescale ( τ R) much smaller than the numerical heating timescale ( τ H) as can be seen in the example 2d3v (two spatial dimensions; three velocity dimensions) PIC simulation of a spatially uniform plasma of density 1.24 × 10 14 m 3 in Fig. 1. Here, we used the explicit leapfrog time integration method with an explicit momentum conserving scheme to simulate the plasma in a 128 × 128 cell Cartesian grid, with periodic boundary conditions, a grid spacing of 0.23 mm, and 40 particles per cell. We observe that an initial waterbag EVDF
f v x , v y , v z = Π v x Π v y Π v z , Π v = 1 2 v c , v c < v < v c , 0 , v > v c | v < v c
(1)
with a cutoff velocity v c corresponding to the velocity of an electron with 4.5 eV of kinetic energy experiences numerical thermalization within a few hundred plasma frequency periods. By comparison, the average electron kinetic energy only increases by ∼2.5% after ten thousand plasma frequency periods. The numerical thermalization results in a rapid increase in the kinetic entropy calculated from the one-dimensional EVDFs in the x and y directions [see Fig. 1(b)]. Although the difference between the kinetic entropy of a Maxwellian distribution and the kinetic entropy of a waterbag distribution is small relative to their absolute value, it provides a very clear signal as it is obtained from the global EVDF of the entire simulation domain. Furthermore, information on thermalization measurement techniques can be found in Sec. III A.
FIG. 1.

Numerical thermalization of an EVDF for a 2d3v electrostatic PIC simulation run using the PIC code EDIPIC-2D.10 (a) The evolution of the one-dimensional EVDF in the x direction as it approaches a Maxwellian, and (b) the “one-dimensional” electron temperature T e , j and kinetic entropy S j in the jth dimension calculated from the one-dimensional EVDFs as shown in Eqs. (2) and (3). The dashed blue line and solid blue line representing the quantities in the spatially unrepresented z dimension lie on top of each other.

FIG. 1.

Numerical thermalization of an EVDF for a 2d3v electrostatic PIC simulation run using the PIC code EDIPIC-2D.10 (a) The evolution of the one-dimensional EVDF in the x direction as it approaches a Maxwellian, and (b) the “one-dimensional” electron temperature T e , j and kinetic entropy S j in the jth dimension calculated from the one-dimensional EVDFs as shown in Eqs. (2) and (3). The dashed blue line and solid blue line representing the quantities in the spatially unrepresented z dimension lie on top of each other.

Close modal
The quantities shown in Fig. 1(b), denoted T e , j and S j, respectively, are one-dimensional measurements of electron temperature and a quantity related to the entropy of the normalized electron distribution function, f ( v ),
T e , j = m e d v j v j 2 f ( v j ) ,
(2)
S j = k B d v j f ( v j ) ln f ( v j ) .
(3)
They are defined while considering the EVDF only as a function of one velocity dimension j x , y , z. This means that the various S j cannot be summed to produce an actual measurement of the kinetic entropy of the full multidimensional EVDF. However, separately tracking this measurement in each dimension illustrates that without Monte Carlo collisions or magnetic fields, which might redirect particle motion out of the simulated spatial plane, the evolution of the EVDF in a 2d3v simulation only occurs in velocity dimensions for which the electric field is defined. Thus, for the 2d3v simulations depicted in Fig. 1, the one-dimensional distribution function in the z direction does not evolve even after twenty thousand plasma frequency periods have passed. This of course will not be the case for 3D PIC, where the ability of particles to extract energy from fields in all three spatial dimensions makes it an attractive tool for modeling kinetic turbulence.11 

Numerical thermalization occurs due to a process, which is analogous to Coulomb collisions in real plasmas and has alternately been described as “numerical collisions.” The basic form of a numerical collision operator is plainly stated in Birdsall and Langdon's textbook.12 It has also recently been rederived including the effects of multiple particle species with lucid and detailed discussion by Touati et al.13 

Because numerical collisions are caused by the inherent graininess of particle-based simulations, they do not disappear in the limit of infinitely small grid spacing and time step. The numerical collision operator is a correction to the Vlasov evolution of plasma in phase space, and to first order, it scales with the inverse of the numerical plasma parameter 1 / N Dmac, where N Dmac is the number of macroparticles per Debye volume. For our purpose, this is defined as N Dmac = n mac λ D in a 1D simulation, N Dmac = n mac λ D 2 in a 2D simulation, and N Dmac = n mac λ D 3 in a 3D simulation with λ D representing the Debye length and n mac the macroparticle number density. Cancellation of a large portion of this first-order term in 1D PIC plasmas can cause the rate of numerical thermalization to depend on a second-order term, which scales as 1 / N Dmac 2, a process known as kinetic blocking.7 The kinetic blocking can be disrupted by the addition of Monte Carlo collisions to the simulation, recovering a first-order term.14–16 

The numerical collision operator acting on species α can be written in a Fokker–Planck form representing drag A α v and diffusion D α v in velocity space
f α t c = v · A α v f α v + 1 2 v v : D α v f α v .
(4)
A α v is a vector parallel to the direction of v, while D α v can be written as a diagonal matrix in a velocity space defined by a basis set of vectors parallel and perpendicular to v. Drag is caused both by the polarization, which occurs when Debye shielding around a moving macroparticle lags behind it and also by the transfer of kinetic energy from the moving macroparticle into the energy of the fluctuating electric field. Diffusion is similarly caused by resonant interactions with these fluctuations; the diffusion coefficient is proportional to an integral over the fluctuation spectrum. This is identical to a description of Coulomb collisions between charged particles in real plasmas, and this is because numerical thermalization can be thought of as Coulomb collisions between macroparticles. Velocity-dependent relaxation rates arise directly from the drag and diffusion coefficients of the Fokker–Planck form.

There are, however, three main differences that distinguish numerical thermalization from the effects of Coulomb collisions between real particles. First, the fluctuation spectrum is modified by aliasing effects due to finite grid resolutions and timesteps. Second, interactions between particles are reduced by the use of shape functions, which give macroparticles finite width. The final possible difference is the reduced number of dimensions, which has a significant impact on 1D simulations as described above.

This paper is divided into five sections. In Sec. II, we will present the drag and diffusion coefficients for a test particle electron in a Maxwellian electron background for 2D electrostatic explicit PIC simulations using a cloud-in-cell (CIC) particle weighting scheme on a uniform Cartesian grid. These coefficients will be compared to the drag and diffusion coefficients of the Landau collision operator. In Sec. III, we will then discuss the relaxation rates related to these coefficients and examine their resemblance to a measured empirical scaling of numerical thermalization time. We will also discuss possible strategies for reducing the rate of numerical thermalization. In Sec. IV, we will assess the relaxation rates in application to various 2D PIC simulations of low-temperature plasmas, including capacitively coupled plasma (CCP) discharges, inductively coupled plasma (ICP) discharges, electron beam-generated plasmas, and hollow cathode discharges. Finally, we will summarize the results and give final recommendations for addressing numerical thermalization in 2D PIC simulations.

We first present the drag and diffusion coefficients from the Fokker–Planck form of the 2D electrostatic PIC numerical collision operator for test particle electrons in a homogeneous quasineutral background plasma with one singly charged ion species. This analysis is carried out for simulations using a uniform, effectively infinite Cartesian grid. Both background electrons and ions are assumed to be Maxwellian, with respective temperatures T e and T i and masses m e and m i. The numerical collision operator was calculated from the results derived by Touati et al.13 for explicit PIC simulations, where we have neglected the effects of finite time differencing. This is done under the assumption that the time step is sufficiently well-resolved with respect to the plasma frequency so that time aliasing modes other than the zeroth mode are negligible. We have also assumed equal macroparticle weighting factors for electrons and ions, resulting in an equal density of electron and ion macroparticles. Furthermore, details regarding the calculation of these Fokker–Planck [Eq. (4)] drag and diffusion coefficients can be found in the  Appendix
A e = β A e / β , D e = β D e / β , D e = V ̂ T · D e · V ̂ , D e = V ̂ T · D e · V ̂ ,
(5)
A e / β = ω p e v T e N Dmac I A e / β , D e / β = ω p e v T e 2 N Dmac I D e / β ,
(6)
I A e / β = 1 + m e / m β 2 5 / 2 π 3 / 2 R v β 3 d k p S k S k p ϵ k · V , k K 2 2 K · k p k · V K · V ̂ k p 3 e k · V 2 / k p 2 R v β 2 ,
(7)
I D e / β = 1 2 5 / 2 π 3 / 2 R v β d k p S k S k p ϵ k · V , k K 2 2 K K k p e k · V 2 / k p 2 R v β 2 ,
(8)
ϵ k · V , k 1 p S k p 2 K · k p 2 K 2 k p 2 Z k · V k p + T e T i Z k · V k p R v β .
(9)
These coefficients are proportional to dimensionless integrals presented as functions of scaled velocity V = v / v T e, where v T e = 2 k B T e / m e is the thermal velocity of the background electrons. The test particle interactions with species β = ( i , e ) are also given as a function of the ratio of the species' thermal velocity with that of the background electrons, so R v β = v T β / v T e = T β m e / T e m β. Wavenumbers and similar quantities have been scaled by the electron Debye length (i.e., k = k λ D e with λ D e = k B T e ε 0 / n q 2) so that the integral over the volume d k = d k x d k y is a dimensionless quantity. The spatial aliasing effects are accounted for through sums over all integers p = ( p x , p y ), which couple aliased wavelengths with corresponding wavenumbers dependent on the spatial grid resolution, k p = k p 2 π / R g where R g = Δ x / λ D e is the grid spacing Δ x scaled by the electron Debye length. Here, we assume Δ x to be equal in both directions. K is a quantity dependent on the particular Maxwell solver used, but here assumed to be equal to the scaled wavenumber k as would be the case for a spectral Maxwell solver. For our calculations, we have chosen the commonly used CIC method for assigning charge to the grid, yielding a shape function with the Fourier transform S k = sinc k x R g / 2 2 sinc k y R g / 2 2. Here, sinc x = sin x x .

As in a real plasma, the contribution to the dielectric function ϵ k · V , k from each species γ is proportional to its plasma frequency squared: ω p γ 2 = n γ q γ 2 / ε 0 m γ. Also the dielectric function is the plasma dispersion function Z ( x ), where Z ( x ) is its derivative evaluated at x. We choose not to approximate ϵ k · V , k with its static approximation at V = 0 based on the results of Okuda and Birdsall.17 They showed that the static shielding approximation (Landau collision operator analog) is not a good approximation of the full Balescu–Lenard collision operator analog for collisions between clouds of charge or for particles in a reduced number of dimensions.

A key feature to note is that both the drag and diffusion coefficients are inversely proportional to N Dmac = n mac λ D 2 = N ppc / R g 2, where N ppc is the number of particles per cell. Unlike real Coulomb collisions in 3D, there is no logarithmic dependence on N Dmac in the numerator. This is because, thanks to the shape functions, there is no need to truncate the integrals at large k values (small wavelengths) in order to achieve convergence.

The numerical drag and diffusion coefficients for a Maxwellian electron background can be compared (see Fig. 2) to drag and diffusion coefficients for the Landau collision operator for a test electron in a Maxwellian electron background
A L e / e = ω p e v T e ln Λ e e 2 N D I A , L e / e , D L e / e = ω p e v T e 2 ln Λ e e 2 N D I D , L e / e ,
(10)
I A , L e / e = 1 2 1 / 2 π erf V e 2 π 1 / 2 V e e V e 2 2 V e 2 ,
(11)
I D , L e / e = 1 2 3 / 2 π erf V e 2 π 1 / 2 V e e V e 2 2 V e 3 ,
(12)
I D , L e / e = 1 2 5 / 2 π 2 V e 2 erf V e erf V e + 2 π 1 / 2 V e e V e 2 2 V e 3 .
(13)
FIG. 2.

Scaled numerical drag and diffusion coefficients at a variety of grid mesh sizes (solid color lines) and drag and diffusion coefficients for the Fokker–Planck form of the Landau collision operator (dashed black lines) for test particle electron interactions with Maxwellian background electrons. Numerical drag coefficients are scaled by ω p e v T e / N Dmac, and the corresponding diffusion coefficients are scaled by ω p e v T e 2 / N Dmac while the drag coefficients of the Landau collision operator are scaled by ω p e v T e ln Λ e e / 2 N D and the corresponding diffusion coefficients are scaled by ω p e v T e 2 ln Λ e e / 2 N D.

FIG. 2.

Scaled numerical drag and diffusion coefficients at a variety of grid mesh sizes (solid color lines) and drag and diffusion coefficients for the Fokker–Planck form of the Landau collision operator (dashed black lines) for test particle electron interactions with Maxwellian background electrons. Numerical drag coefficients are scaled by ω p e v T e / N Dmac, and the corresponding diffusion coefficients are scaled by ω p e v T e 2 / N Dmac while the drag coefficients of the Landau collision operator are scaled by ω p e v T e ln Λ e e / 2 N D and the corresponding diffusion coefficients are scaled by ω p e v T e 2 ln Λ e e / 2 N D.

Close modal

Here, V e = v / v Te is the velocity scaled by the electron thermal velocity, N D = n λ D 3 is the number of real particles per Debye cube, and ln Λ e e = ln 4 π N D is the Coulomb logarithm.

The numerical drag and diffusion coefficients depend not only on the magnitude of V e but also its orientation with respect to the grid. However, this dependence was found to be negligible when grid spacings were set equal in both directions, and in Fig. 2, the orientation of the test particle velocity was taken to be along one of the grid axes.

The coefficients shown in Fig. 2 can be compared to those derived by Abraham-Shrauner18 for single-species point particles in 2D with a Lorentz background distribution, those derived by Okuda and Birdsall17 for plasmas of single-species Gaussian charge clouds in 2D, and those derived by Reynolds et al.19 for multispecies point particles in 2D. Similar to real Coulomb collisions, the drag coefficient is linear at v v T e. The parallel and perpendicular diffusion coefficients are equal at v = 0, with perpendicular diffusion (corresponding to pitch angle scattering) falling off more slowly as the test particle velocity increases.

As shown by Okuda and Birdsall, increasing the width of the particle shape functions (here accomplished by increasing the size of the grid spacing) reduces numerical drag and diffusion. Decreasing the width of the particle shape function causes the numerical drag and diffusion coefficients to approach that of a point particle.17 However, at v v T e, the effect of the finite macroparticle width is reduced. Okuda and Birdsall attribute this to the fact that fast-moving charge clouds resonantly interact with background fluctuations at wavelengths which scale linearly with the test particle velocity; thus, there is some v at which the wavelength of the resonant interaction is larger than the width of the shape function and the behavior of the charge cloud is indistinguishable from that of a point particle.17 

We can see in Fig. 2 that drag and diffusion coefficients of the Landau collision operator appear to fall off much more rapidly with increasing velocity than the 2D numerical drag and diffusion coefficients. However, the behavior of the 2D coefficients at very large velocities has been shown to have the same dependence on v as the coefficients in 3D in the limit that v v T.19 

From these coefficients, we can define the three relaxation processes and their corresponding rates as they are usually written: slowing down, parallel diffusion, and perpendicular diffusion
ν s e / β = A e / β v , ν e / β = D e / β v 2 , ν e / β = D e / β v 2 .
(14)
It is immediately apparent that any thermalization time based on these processes will scale with τ R N D , mac / ω pe. While there is no relaxation timescale which universally applies to every initial electron distribution function, there are nevertheless many methods, which have been used to measure some definition of τ R in PIC plasmas.

Many of these measurement methods involve setting up an initial waterbag EVDF and tracking its evolution via a parameter, which measures the distribution's resemblance to a Maxwellian. Examples of these parameters include the value of the velocity distribution function at zero velocity,20 a measure of the kinetic entropy of the system,9,21,22 the coefficients of an expansion of the distribution function in eigenfunctions of the Lenard–Bernstein collision operator,14,15,23 and a simple standard statistical chi-squared fit to a Maxwellian.15,16

The thermalization time has also often been equated to the slowing-down timescale ( τ s) associated with the drag on a slow-moving test particle.8,24 This resulted in initial errors in estimating the thermalization time in a 1D single-species plasma24 due to the necessity of distinguishing between drag and diffusion timescales for a test particle ( τ s , τ , τ N Dmac) and processes, which evolve the total velocity distribution function ( τ R N Dmac 2). There is no such distinction in 2D, and so, the scaling of this test particle slowing-down time will accurately reflect the scaling of the thermalization time ( τ s , τ R N Dmac).8 

As mentioned before, complications arise in 2D PIC simulations if there is an applied magnetic field that points out of the plane of simulation. Numerical collisionality is enhanced as the gyromotion of particles can cause Coulomb interactions between particles to repeat in a reinforcing manner. The inability of the particles to travel away from each other along a nonexistent third dimension parallel to the magnetic field lines can cause the scaling of the thermalization time under these circumstances to be τ R N Dmac .21,22,25

We can compare our analytical estimates for the average test particle drag timescale in a 2D PIC plasma with a CIC shape function to Hockney's8 empirical estimate of the same
τ ¯ s Hockney = τ p e N Dmac K 1 1 + Δ x λ D 2 , K 1 = 0.98 ± 0.20 .
(15)
Considering only drag due to electron–electron interactions, we can calculate this average drag timescale by integrating over velocity
τ ¯ s e / e = d v ν s e / e 1 π v T e 2 e v 2 / v T e 2 .
(16)
The comparison in Fig. 3 shows that this is reasonably close to Hockney's empirical estimate, which provides an order of magnitude estimate for the timescale of numerical thermalization.
FIG. 3.

Numerical collision timescales for 2D PIC with a CIC scheme. Comparison of timescales derived from analytical theory (blue solid curve) with the empirical scaling given by Hockney (black dashed curve).8 Hockney's estimates accounted for large grid spacings up to 32 times the Debye length.

FIG. 3.

Numerical collision timescales for 2D PIC with a CIC scheme. Comparison of timescales derived from analytical theory (blue solid curve) with the empirical scaling given by Hockney (black dashed curve).8 Hockney's estimates accounted for large grid spacings up to 32 times the Debye length.

Close modal

We should note that Hockney's method8 of analyzing numerical collisional time scales by conducting “numerical experiments” to measure average electron slowing-down and deflection times is an invaluable test for any new particle-in-cell code. In cases where a numerical collision operator cannot be easily analytically defined direct measurement is especially useful. Simulations with unstructured grids fall into this category, and measurement of slowing-down and deflection times in these cases offer the best insight into numerical thermalization rates. Examples of this type of analysis for unstructured grids can be found in work by Gatsonis and Spirkin26 and work by Averkin and Gatsonis.27 

Given that numerical thermalization happens on timescales roughly comparable to N Dmac τ p e, the most obvious way to delay the onset of numerical thermalization is to increase the number density of macroparticles. Increasing the grid spacing size with respect to the Debye length produces effectively larger macroparticle shape functions, reducing the Coulomb interactions between macroparticles and delaying numerical thermalization. Significant reductions in thermalization time are only achieved by under-resolving the Debye length, however, and may result in numerical instabilities in explicit schemes.

Implicit or energy-conserving schemes28 can be used to create larger grid sizes while avoiding numerical heating due to grid aliasing, but one should note that the thermalization rate depends linearly on the density of the macroparticles, so that simulations with large grid length scales should use larger numbers of macroparticles per cell to fully take advantage of the reduced numerical collisions offered by broad shape functions. The numerical collision operators shown in Sec. II strictly speaking only apply to the explicit schemes for which they were derived, but as can be seen in Fig. 4(a), the numerical thermalization rates are comparable for a direct implicit and explicit scheme with a well-resolved grid resolution, time step, and number of particles per cell. There we can also see that increasing the grid length scale above the Debye length while maintaining the same density of macroparticles will result in a reduced rate of numerical thermalization. Numerical heating (or cooling) in the direct implicit scheme with an under-resolved Debye length is not ignorable, but it can be made negligible using a judicious choice of time step; according to Sun et al.,29 this “optimal path” of near-zero numerical heating for direct implicit algorithms occurs when the ratio between ω pe d t and Δ x / λ D e is roughly 0.1. As can be seen in Fig. 4(b), the direct implicit simulation with a grid spacing of Δ x = 3 λ D e does experience a modest ∼6% increase in the kinetic energy of the electrons over ten thousand plasma periods—if stricter control of heating is required, the time step can be refined further.

FIG. 4.

The evolution of the “one-dimensional” electron kinetic entropy (a) and electron temperature (b) in the x direction [see Eqs. (2) and (3)] for explicit and direct implicit PIC simulations in EDIPIC-2D with initial waterbag EVDFs.

FIG. 4.

The evolution of the “one-dimensional” electron kinetic entropy (a) and electron temperature (b) in the x direction [see Eqs. (2) and (3)] for explicit and direct implicit PIC simulations in EDIPIC-2D with initial waterbag EVDFs.

Close modal

The simulated densities and initial EVDFs are identical to those described in Fig. 1. The explicit and direct implicit simulations with a well-resolved Debye length ( Δ x = 0.2 λ D) contain 40 particles per cell, while the direct implicit case with an unresolved Debye length ( Δ x = 3 λ D) contains 9000 particles per cell, leaving the density of macroparticles constant between the cases.

The setup and the quantities calculated are identical to those discussed in Sec. I, though the one-dimensional kinetic entropy in Fig. 4(a) is presented as a fractional value of the entropy of a Maxwellian with the average kinetic energy of the EVDF at that point in time. In this way, we compare increases in entropy due to numerical thermalization more directly, ignoring increases in entropy due to numerical heating.

Aside from using implicit schemes with large grid spacings, another option that effectively increases the width of the macroparticle charge clouds and thus delays thermalization is to use an interpolation scheme of a higher order than the CIC scheme analyzed in this work. There is also some indication that unstructured grids may result in reduced thermalization times; the slowing-down and deflection times measured by Averkin and Gatsonis for the 3D electrostatic PIC code EUPIC with unstructured tetrahedral grids were in fact larger than the measured heating times.27 This should prompt further investigation into the use of such grids as a potential avenue to minimize numerical thermalization.

Finally, for some applications, it is reasonable to reduce the effective plasma frequency by scaling the permittivity of free space, ε 0, by some large factor. The latter strategy will, of course, only be useful in situations where there is certainty that other important phenomena such as instabilities are of little interest. Many low-temperature plasma sources fall into this category, and in some of these cases, it is possible to sufficiently reduce the numerical collision rate below that of real Coulomb collisions by increasing ε 0. See, for example, the hollow cathode simulations detailed in Sec. IV B.

A simple glance at the scaling in Fig. 2 is sufficient to infer that reducing the effects of numerical thermalization below that of real Coulomb collisions by simply increasing the number of macroparticles is often computationally expensive to attempt. To make the numerical processes in 2D PIC roughly comparable to the drag and diffusion caused by true Coulomb collisions (at velocities small compared to the electron thermal velocity), it would require that a macroparticle represent only a handful of real particles: N Dmac 2 N D / ln Λ e e. At higher velocities, reducing the numerical drag and diffusion coefficients below that of real Coulomb collisions may require an even larger number of macroparticles.

While the situation may seem hopeless for schemes with a resolved Debye radius, if relaxation or partial thermalization due to real Coulomb collisions is an important process in the plasma which is to be simulated, the numerical thermalization can to some degree act as a substitute for real Coulomb collisions. However, it is important to compare the rate of this rapid thermalization to other relevant timescales to assess whether it will affect the results.

In Sec. IV, we will provide several examples of exactly this analysis for a variety of low-temperature plasma simulations in 2D PIC.

Electron beam-generated plasmas are of interest as a source of plasmas with large densities ( 10 15 10 17 m 3) and low-electron temperatures ( 1 eV) suitable for materials processing applications.2,30,31 This can produce a large flux of low-energy ( < 5 eV) ions potentially useful for atomic layer etching or deposition.

Here, we examine the effects of numerical thermalization on a plasma generated by a beam of electrons, such as that described by Rauf et al.2 The explicit 2d3v open-source PIC code EDIPIC-2D10 was used to simulate a 2 keV beam passing through the center of a 9 cm by 12 cm chamber with grounded conducting boundaries. A symmetry plane passes through the center of the beam, seen in Fig. 5(a). The beam ionizes background neutral argon gas and an applied magnetic field runs parallel to the direction of the beam. The simulations were run for 800 μs, at which point a steady state was achieved.

FIG. 5.

2d3v PIC simulation of electron beam-generated plasma. (a) Colormap of electric potential in eV. The spatial region for EEDF averaging is outlined by a black rectangle. (b) The one-dimensional EEDFs for electrons within the outlined spatial region, plotted against the kinetic energy in the corresponding direction times the sign of the velocity: sgn ( v j ) m e 2 v j 2 with j { x , y , z }. The black dashed line shows the theoretical EEDF of a Maxwellian distribution with a temperature of 1.1 eV. (c) Electron density profiles along field lines at five different distances from the center of the beam. (d) Electron temperature profiles along field lines at five different distances from the center of the beam [left boundary of (a)].

FIG. 5.

2d3v PIC simulation of electron beam-generated plasma. (a) Colormap of electric potential in eV. The spatial region for EEDF averaging is outlined by a black rectangle. (b) The one-dimensional EEDFs for electrons within the outlined spatial region, plotted against the kinetic energy in the corresponding direction times the sign of the velocity: sgn ( v j ) m e 2 v j 2 with j { x , y , z }. The black dashed line shows the theoretical EEDF of a Maxwellian distribution with a temperature of 1.1 eV. (c) Electron density profiles along field lines at five different distances from the center of the beam. (d) Electron temperature profiles along field lines at five different distances from the center of the beam [left boundary of (a)].

Close modal

Taking as an example a background gas pressure of 20 mTorr and an applied magnetic field strength of 100 Gauss, we can examine the EVDF within a region a few centimeters away from the beam. From Fig. 5(b), it is clear to see that the electrons have a Maxwellian distribution function, though the EEDF is depleted somewhat at energies sufficient to overcome the large sheath potential. The density rapidly falls off with increasing distance from the beam [Fig. 5(c)], while the electron temperature decreases steadily and remains essentially constant along field lines [Fig. 5(d)].

1. Timescale analysis: Electron beam-generated plasma

It should be emphasized that Coulomb collisions were not implemented in this simulation; the thermalization process is entirely numerical. To understand how it has impacted the results of the simulation, we should compare the timescale for numerical thermalization to the time required for electrons to diffuse away from the beam. These in turn should be compared to timescales associated with thermalization via real Coulomb collisions.

Based on the electron fluid velocity in the x-direction v ¯ x observed in the simulations, we can estimate the time required for an electron to travel 2 cm away from the beam center along the x direction to be t x = 250 μ s. On the other hand, an electron with sufficient kinetic energy (say 8 eV) to immediately escape along the magnetic field lines to the upper or lower wall will do so in a time less than a few dozen nanoseconds ( t y = 27 ns), assuming it does not experience an elastic collision with a neutral which redirects it. For electrons with a kinetic energy between 1 and 8 eV, the electron–neutral elastic collision time ranges from t e Ar = 6.75 ns to t e Ar = 180 ns. Here, t e Ar = 1 / ν e Ar is the inverse of the electron elastic collision rate, where ν e Ar = n Ar σ e l E 2 E / m e depends on the background gas density n Ar and the collision cross section σ e l E for an electron with kinetic energy E.

The numerical thermalization timescale falls between the two travel times. With an electron temperature of 1.3 eV and plasma density of 2 × 10 15 m 3, the electron Debye length near the beam is 0.19 mm and the plasma frequency is 2.5 × 10 9 s 1. The grid resolution is 0.21 mm, and there are roughly 200 macroparticles per cell. Using Hockney's empirical estimate of the thermalization time, we would have τ R num  1.1  μ s. Clearly, the thermalization occurs well before the electrons have diffused beyond the beam region, as confirmed in Fig. 4. However, it occurs slowly enough that the fast electrons with a kinetic energy larger than the sheath potential will almost certainly have time to escape before they are significantly slowed down.

How does this numerical thermalization compare to real thermalization timescales via Coulomb collisions? Taking as our estimate for the Coulomb thermalization timescale the electron–electron slowing-down time for a test particle traveling at the electron thermal velocity in a Maxwellian plasma, we have τ R real 4.5 μ s. While this is larger than the numerical thermalization time, it is still quite small compared to the time required for the electrons to diffuse a few centimeters away from the beam. It is also large compared to the time required for fast electrons to escape to the upper and lower walls along magnetic field lines. This analysis can be summed up by the following relation:
t x > τ R real > τ R num > t e A r > t y .
(17)
Because the thermalization caused by numerical noise happens on a timescale that can be ordered with respect to travel/escape timescales in a manner similar to the real thermalization time due to Coulomb collisions, we can conclude that the numerical thermalization in the beam region will more or less reproduce the effects of real Coulomb collisions, though it occurs at a more rapid rate.

2. Possible effects of numerical noise: Electron beam-generated plasma

There are some complications when we consider the flux of electrons to the wall in regions far from the beam. In these areas, the cold trapped electrons contained by the large sheath potentials must escape to the wall by diffusing outwards in velocity space, as there is no other mechanism for them to gain kinetic energy. The fact that the EEDF is depleted above the energy required to overcome the sheath potential indicates that the flux to the wall is limited by the flux of electrons into the untrapped region of phase space. (This is analogous to the loss cone of a magnetic mirror with an additional confining electric potential; here, the region is bounded by two planes in velocity space since the magnetic field strength is constant.) The only mechanism for electrons to gain energy is through parallel diffusion in velocity space—which occurs more rapidly due to numerical noise than the analogous process caused by real Coulomb collisions.

If we were to hypothetically increase the number of macroparticles per cell and thus decrease the rate of numerical collisions, we may at one point achieve a rate of refilling the untrapped region of phase space comparable to the rate caused by real Coulomb collisions. However, without implementing Coulomb collisions into the code, there would be no indication that this state had been achieved—and further increasing the number of particles per cell would cause the outward diffusion of electrons in velocity space to occur ever more slowly.

Here, we analyze the impact of numerical thermalization on a 2D PIC simulation of a hollow cathode discharge in an argon plasma, also performed using EDIPIC-2D.10 Electrons are emitted by the cathode with 50 A of current into a narrow channel 1 cm in diameter and 10 cm in length, ionizing a background argon gas at a pressure of 93 mTorr. Electrons produced via ionization expand out of the channel and into the larger bulk plasma, accelerating toward the anode above the channel opening, which is held at 30 V. In this simulation, the permittivity of free space was scaled by a factor of 1059, reducing the plasma frequency by a factor of about 32.5. The electron velocity distribution function within the channel after a total simulation time of 200  μ s is shown in Fig. 6(b), averaged over the channel area shown in Fig. 6(a).

FIG. 6.

2d3v PIC simulation of a hollow cathode plasma. (a) Colormap of electric potential in eV. The spatial region for EEDF averaging is outlined by a black rectangle. The left boundary is a reflective symmetry axis. (b) The EEDF for electrons within the outlined spatial region. The black dashed line shows the theoretical EEDF of a Maxwellian distribution with a temperature of 2.7 eV. (c) x-averaged electron density along the length of the channel, within the spatial region outlined by the black rectangle in (a). (d) x-averaged electron temperature (mean kinetic energy) along the length of the channel within the spatial region outlined by the black rectangle in (a).

FIG. 6.

2d3v PIC simulation of a hollow cathode plasma. (a) Colormap of electric potential in eV. The spatial region for EEDF averaging is outlined by a black rectangle. The left boundary is a reflective symmetry axis. (b) The EEDF for electrons within the outlined spatial region. The black dashed line shows the theoretical EEDF of a Maxwellian distribution with a temperature of 2.7 eV. (c) x-averaged electron density along the length of the channel, within the spatial region outlined by the black rectangle in (a). (d) x-averaged electron temperature (mean kinetic energy) along the length of the channel within the spatial region outlined by the black rectangle in (a).

Close modal

Within the channel, a population of cold electrons at energies less than 10 eV is well described by a Maxwellian distribution with a temperature of 2.7 eV. Above this energy, the EEDF falls off faster than a Maxwellian until the peaks representing a population of fast electrons injected from the cathode is visible at roughly 30 eV. These peaks are especially pronounced in the one-dimensional EEDF obtained from the velocity in the x direction [see Fig. 6(b)] as the strong electric field from the sheath potential at the cathode wall accelerates injected electrons across the channel. However, the one-dimensional EEDFs in the y and z directions are also enhanced at high energies due to elastic scattering of the accelerated electrons.

1. Timescale analysis: Hollow cathode plasma

Based on the electron fluid velocity in the y-direction v ¯ y observed in the simulations, we can estimate the time required for an electron to travel the length of the channel (10 cm) to be t y = 13 μ s. A negligible number of electrons return to the cathode wall, as they would require nearly 35 eV of kinetic energy in the x direction. Such an electron unimpeded by collisions would travel to the wall in t x = 1.4 ns. The elastic collision timescale for a thermal electron at 2.7 eV in this simulation is t e A r = 9 ns, while the same for 35 eV electron is only t e Ar = 1.5 ns.

With an electron density of roughly 10 20 m 3 and a temperature of 2.7 eV in the channel region, the plasma frequency is about 5.6 × 10 11 s 1 and the Debye length is 1.2 μ m. However, due to the scaling of the permittivity of free space, the “simulated” plasma frequency and Debye length are 1.7 × 10 10 s 1 and 38 μ m, respectively. Based on a macroparticle density of n mac = 3.12 × 10 10 m 2 and a cell width of 125  μ m, this yields a thermalization time estimate of τ R num 200 ns. In comparison, the thermalization timescale from (hypothetical, not implemented) Coulomb collisions (estimated as in Sec. IV A) is τ R real 1.6 ns. We can once again order these timescales:
t y > τ R num > t e Ar t x τ R real .
(18)
In this case, the scaling of the permittivity of free space has actually reduced the rate of thermalization below that caused by real Coulomb collisions. Both simulated and real thermalization rates occur faster than physical transport of electrons along the channel.

2. Inclusion of real electron–electron collisions

As shown in Sec. IV B 1, scaling the permittivity of free space reduces numerical collisions below the rate of real electron–electron collisions. Thus, simulations that include Coulomb collisions will no longer be overwhelmed by numerical thermalization and depend wholly on real processes. This can be seen in the channel EEDF of the following hollow cathode simulation, using a similar geometry to the setup described previously, identical background gas pressure, and with an emitted current density at the cathode of 2.06  kA m 2. Thermalization here is driven by implemented electron–electron Coulomb collisions using the Nanbu method,32 and numerical collisions are reduced through scaling of the permittivity of free space. The electron–electron collision frequency in this case has been made faster than the estimated numerical thermalization rate by a factor of 20.

This EEDF in Fig. 7 is clearly non-Maxwellian; the high-energy electrons accelerated across the sheath from the cathode maintain large peaks in the wings, even more pronounced than in Fig. 6(b). Though the cold electrons have begun to thermalize, we can attribute this to the actual impact of Coulomb collisions rather than the numerical thermalization, which would otherwise dominate the process if permittivity of free space had not been scaled to mitigate it.

FIG. 7.

The one-dimensional EEDFs for electrons within the channel of the hollow cathode. The black dashed line shows the theoretical EEDF of a Maxwellian distribution with a temperature of 0.6 eV.

FIG. 7.

The one-dimensional EEDFs for electrons within the channel of the hollow cathode. The black dashed line shows the theoretical EEDF of a Maxwellian distribution with a temperature of 0.6 eV.

Close modal

The following 2D PIC simulation is that of an inductively coupled plasma operating in argon at 5 mTorr background gas pressure. The following electromagnetic simulation was performed using a direct implicit scheme and Darwin algorithm in an adaptation of the EDIPIC-2D code. The antennae within the dielectric carry a 120 A current that oscillates at a frequency of 2 MHz, and the non-dielectric boundaries of the simulation space are grounded metal walls. The simulation ran for 100  μ s before the snapshots depicted in Fig. 8 were taken.

FIG. 8.

2d3v PIC simulation of an inductively coupled plasma. (a) Colormap of electric potential in eV. Note that the bounding conductive walls are grounded, and the sheath is not resolved. The spatial region for EEDF averaging is outlined by a black rectangle. (b) The EEDF for electrons within the outlined spatial region. The black dashed line shows the theoretical EEDF of a Maxwellian distribution with a temperature of 1.5 eV. (c) Electron density profiles at three different distances from the dielectric. (d) Electron temperature profiles at three different distances from the dielectric.

FIG. 8.

2d3v PIC simulation of an inductively coupled plasma. (a) Colormap of electric potential in eV. Note that the bounding conductive walls are grounded, and the sheath is not resolved. The spatial region for EEDF averaging is outlined by a black rectangle. (b) The EEDF for electrons within the outlined spatial region. The black dashed line shows the theoretical EEDF of a Maxwellian distribution with a temperature of 1.5 eV. (c) Electron density profiles at three different distances from the dielectric. (d) Electron temperature profiles at three different distances from the dielectric.

Close modal

1. Timescale analysis: Inductively coupled plasma

A thermal 1.5 eV electron would take t x , y 210 ns to travel 15 cm without collisions, while the electron–neutral elastic collision time for such an electron is t e Ar = 410 ns. With an electron density of 8.7 × 10 16 m 3 and electron temperature of about 1.5 eV, the plasma frequency within the outlined region in Fig. 8 is 1.7 × 10 10 s 1 and the Debye length is 31 μ m. Considering a cell size of 1.66 mm and 1000 particles per cell, this yields a thermalization time based on Hockney's estimate of τ R 380 ns. Thermalization due to real Coulomb collisions is estimated to occur on a timescale of τ R 590 ns
τ R rea l > τ R num t e Ar > t x , y .
(19)
The numerical thermalization thus reasonably approximates the effects of real Coulomb collisions, and for trapped cold electrons will occur roughly on the same timescale as it takes for electrons to travel the length of the system. As with the electron beam simulation, complications arise when we consider the fast-moving electrons in the depleted wings of the distribution at kinetic energies sufficient to overcome the sheath potential. For those electrons, numerical diffusion in velocity space is likely occurring faster than it would be if it were caused by real Coulomb collisions.

Capacitively coupled plasmas are frequently used in materials processing applications. In the following 2d3v simulation of a capacitively coupled argon plasma with a background gas pressure of 10 mTorr, we examine the potential and EEDF distribution averaged over the last 1.5  μ s of a 295- μ s simulation time (see Fig. 9). A powered electrode at the bottom of the simulation space [colored gray in Fig. 9(a)] is pulsed at a frequency of 13.56 MHz and is separated from the grounded outer walls by a dielectric (in blue). The amplitude of the applied voltage is 100 V. The simulation was run using the low-temperature particle-in-cell code (LTP-PIC), a massively parallel, GPU accelerated particle-in-cell code specifically designed for modeling low-temperature plasma devices.29 The code has been benchmarked in the collisionless33 and collisional regime.3 

FIG. 9.

2d3v simulation of a capacitively coupled plasma, averaged over 1.5  μ s. (a) Colormap of time-averaged electric potential in eV. The bounding walls are grounded, the powered electrode is in gray, and the dielectric is in blue. The spatial region for EEDF averaging is outlined by a black rectangle. (b) The EEDF for electrons within the outlined spatial region. The black dashed line shows the theoretical EEDF of a Maxwellian distribution with a temperature of 0.5 eV, while the black dotted line shows an EEDF with a temperature of 3.5 eV. (c) Electron density profiles at four different distances from the powered electrode. (d) Electron temperature profiles at four different distances from the powered electrode.

FIG. 9.

2d3v simulation of a capacitively coupled plasma, averaged over 1.5  μ s. (a) Colormap of time-averaged electric potential in eV. The bounding walls are grounded, the powered electrode is in gray, and the dielectric is in blue. The spatial region for EEDF averaging is outlined by a black rectangle. (b) The EEDF for electrons within the outlined spatial region. The black dashed line shows the theoretical EEDF of a Maxwellian distribution with a temperature of 0.5 eV, while the black dotted line shows an EEDF with a temperature of 3.5 eV. (c) Electron density profiles at four different distances from the powered electrode. (d) Electron temperature profiles at four different distances from the powered electrode.

Close modal
The time-averaged EEDF from the spatial region outlined in black in Fig. 9(a) shows a characteristic bi-Maxwellian character of low-pressure CCP discharges34 [see Fig. 7(b)]. Cold electrons are trapped by the relatively flat potential, while the hotter population experiences heating through interactions with the fluctuating sheath near the electrode. Unlike the EEDFs described by Godyak et al.,34 the EEDF produced in this simulation does not sharply drop near the excitation energy of argon at 11.5 eV. The densities and temperatures of the high-temperature and low-temperature electron populations can be obtained from the best fit in the high-energy and low-energy regions as in He et al.35 
f E = f low E + f high E = A e E / T e \_ low + B e E / T e \_ high .
(20)
T e \_ high can be immediately extracted from the fit to the distribution at higher energies, where e E / T e \_ high e E / T e \_ low, while the low-temperature population can then be defined as f low E = f E f high E and the relative density and temperature extracted accordingly. Thus, the distribution shown in Fig. 7(b) yields a high-temperature electron population with T e \_ high = 3.5 eV and a low-temperature population with T e \_ low = 0.6 eV, the former population having a number density roughly 80 times that of the latter.

1. Timescale analysis: Capacitively coupled plasma

The numerical thermalization clearly does not occur so quickly that electrons have been thermalized to a uniform temperature. Electrons with an energy of 3.5 eV, unimpeded by collisions, will travel 5 cm in t x , y = 45 ns. The electron–neutral elastic collision time for an electron with 3.5 eV of kinetic energy is about t e Ar = 56 ns. One oscillation of the applied voltage occurs over a time period of about τ RF = 460 ns.

Considering the simulated time-averaged electron density in the region of interest is 1.1 × 10 15 m 3, the plasma frequency is 1.9 × 10 9 s 1. Using a temperature of 3.5 eV, the Debye length is 0.42 mm. With a grid resolution of 0.1  mm and a macroparticle weighting of fifty thousand simulated electrons per macroparticle, we might estimate the thermalization time as τ R num 1.4 μ s. In contrast, the timescale of thermalization due to real Coulomb collisions is τ R real 130 μ s. This results in the ordering
τ R r eal > τ R num > τ R F > t e Ar t x , y .
(21)
Numerical thermalization occurs two orders of magnitude faster than real thermalization. Though the bi-Maxwellian EEDF is produced, it is possible that the more rapid numerical thermalization has smoothed out the very high-energy portion of the EEDF so that there is no rapid drop off beyond the electron excitation energy at 11.5 eV.

It has been observed by Vass et al. that for some (one-dimensional) PIC simulations of low-pressure RF discharges, regardless of how many particles per cell are used, the simulations will not converge unless Coulomb collisions are incorporated into the code.4 (Notably, the dependence of the EEDF at low energies on electron–electron Coulomb collisions in RF discharges has been discussed analytically by Kaganovich and Tsendin,36 and the generation of cold trapped electrons in this scenario is sensitive to the rate of Coulomb collisions.37) As discussed for beam electrons, this is likely due to the fact that the diffusion of electrons in velocity space will at some point depend on numerical collisions unless the effective numerical collision rate is made smaller than that of implemented electron–electron Coulomb collisions.

We have calculated the drag and diffusion coefficients for the Fokker–Planck form of the numerical collision operator for 2D electrostatic PIC simulations employing the commonly used CIC particle weighting scheme. It is explained that these drag and diffusion coefficients can be used to assess numerical thermalization timescales, which generally depend on the number of macroparticles per Debye volume. Via simple scaling, it is easily demonstrated that the number of macroparticles required to suppress the rate of numerical thermalization below the rate of thermalization caused by Coulomb collisions is quite large, necessitating that each macroparticle represent only a handful of real particles. This requirement may be achievable if the plasma that is being simulated has a sufficiently small plasma parameter.

Grid spacings larger than the Debye length can reduce the numerical drag and diffusion coefficients by increasing the size of the macroparticle shape functions, but these reductions will be negligible for electrons traveling at many times the electron thermal velocity, v T e. If the simulated permittivity of free space is increased by some large factor, then the effective simulated plasma frequency and thus the thermal relaxation rate will also be reduced.

The slowing-down timescale associated with the drag coefficient is shown to qualitatively match the empirical formula developed by Hockney in 1970.8 While Hockney recommended that the thermalization timescale should be used to set an upper limit on the length of the simulation, this need not be the case if one can ensure that the timescale of numerical thermalization can be ordered with other timescales that affect the EVDF in the same way that the real thermalization time due to Coulomb collisions would be ordered. Practical examples of this thermalization timescale investigation have been given for a variety of low-temperature plasma simulations.

One feature of note for both the electron beam-generated plasma and ICP simulations is that the EEDF is depleted at energies larger than the sheath potential. In order for cold electrons to escape to the wall, they must be able to diffuse outward in velocity space. Without implemented Coulomb collisions, the only mechanism by which cold electrons can gain energy is through numerical parallel diffusion in velocity space.

For electrons with a kinetic energy just below that of the sheath potential, this numerical parallel diffusion in two-dimensional simulations often happens much more rapidly than an analogous process caused by Coulomb collisions. Because of this, the loss region in velocity space is more rapidly refilled, and the electron wall fluxes are likely higher than they otherwise would be if the parallel diffusion were caused by Coulomb collisions between real particles. This is something to bear in mind when using PIC codes to simulate plasmas, which contain populations of trapped cold electrons.

The research described in this paper was supported by the U.S. Department of Energy under Laboratory Directed Research & Development (LDRD) program at Princeton Plasma Physics Laboratory (PPPL) and the PPPL CRADA agreement with Applied Materials entitled “Two-Dimensional Modeling of Plasma Processing Reactors.” PPPL is a national laboratory operated by Princeton University for the U.S. Department of Energy under Prime Contract No. DE-AC02-09CH11466.

Resources of the National Energy Research Scientific Computing Center (NERSC) were used for high-performance computation.

The authors have no conflicts to disclose.

Sierra Jubin: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Andrew Tasman Powis: Investigation (equal). Willca Villafana: Investigation (equal); Writing – review & editing (equal). Dmytro Sydorenko: Investigation (equal). Shahid Rauf: Investigation (equal). Alexander V. Khrabrov: Investigation (supporting). Salman Sarwar: Formal analysis (supporting). Igor D. Kaganovich: Conceptualization (equal); Supervision (equal); Writing – review & editing (equal).

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

Here, we give an expression for the numerical collision operator for species α in an infinite homogeneous PIC plasma with a 2D Cartesian grid, as derived by Touati et al.13 We assume that the time step is sufficiently small in comparison with the electron plasma period that aliasing effects due to the time step can be neglected, that the grid spacing is equal in both spatial directions, and that a CIC scheme has been used for particle shape functions
f α t c = β n β m α v · d v Q α β v , v · δ N β m α v δ N α m β v f α v f β v ,
(A1)
Q α β v , v d k 2 π 2 K K p S ( k ) S ( k p ) q α q β ε 0 ϵ k · v , k K 2 2 π δ k · v k p · v ,
(A2)
ϵ ω , k 1 + γ ω p γ 2 K 2 p S ( k p ) 2 d v f γ v · K ω k p · v ,
(A3)
k p = k p k g , k g = 2 π Δ x , S k = sinc k x Δ x / 2 2 sinc k y Δ x / 2 2 .
(A4)
This collision operator is written such that particle density is factored out of the velocity distribution functions, so 1 = d v f α v. The macroparticle weighting factor, δ N α, scales the macroparticle densities ( n α , mac), masses ( M α), and charges ( Q α) to that of real particles such that δ N α = n α / n α , mac = M α / m α = Q α / q α.

Vector quantities have two directions; volume integrals are defined as one would expect; for example, d k represents d k x d k y and similarly for d v . Aliasing effects are accounted for through the sums over all integers p = ( p x , p y ), which couple various aliased wavelengths and frequencies, and Δ x is the grid resolution.

The vector K depends on the specific Maxwell solver used. Touati et al. give the following expressions for its jth component when a spectral Maxwell solver or alternately the second-order finite-difference time-domain (FDTD) method of Yee38 coupled with the charge conserving scheme of Villasenor and Buneman39 is implemented
K j = k j for the spectral Maxwell solver , k j sinc k j Δ j / 2 for the second - order FDTD Maxwell solver .
We assume quasineutrality and equal macroparticle weighting of the electrons and ions limiting our consideration to one singly charged ion species: n β = n and δ N β = δ N = n / n mac. We can then write
f α t c = v · A α 1 v f α v + 1 2 v · D α ( v ) · f α v v ,
(A5)
A α 1 v = β n 2 n mac m α m β d v Q α β v , v · f β v v ,
(A6)
D α v = β 2 n 2 n mac m α 2 d v Q α β v , v f β v .
(A7)
This is a common way of writing the collision operator frequently referred to as the Fokker–Planck equation. It allows for better comparison with previous results regarding drag and diffusion in two-dimensional plasmas.17–19 However, while A α 1 ( v ) is the largest portion of the electron–electron drag term at small velocities (reflecting polarization drag), the dynamical friction drag term has been neglected
A α 2 v = 1 2 v · D α v , A α v = A α 1 v + A α 2 v ,
(A8)
A α v = β n 2 n mac m α 2 1 + m α m β d v Q α β v , v · β v v .
(A9)
When dynamical friction is incorporated, we can write the proper Fokker–Planck form with A α v representing the drag coefficient
f α t c = v · A α v f α v + 1 2 v v : D α v f α v .
(A10)
We neglect the effects of discrete time steps, assuming that the time step is sufficiently small to make a negligible impact. We then evaluate the drag and diffusion coefficients for a test electron in a background of Maxwellian electrons and ions with distributions: f β v = 1 π v T β 2 e v 2 / v T β 2.
Writing A e 1 v = β A 1 e / β v and D e v = β D e / β v, we have
A e / β v = ω p e 4 2 π 3 / 2 n mac 1 + m e m β 1 v T β 3 d k p S k S k p ϵ k · v , k K 2 2 K K · k p k · v k p 3 e k · v 2 / k p 2 v T β 2 ,
(A11)
D e / β v = ω p e 4 2 π 3 / 2 n mac 1 v T β d k p S k S k p ϵ k · v , k K 2 2 K K k p e k · v 2 / k p 2 v T β 2
(A12)
with plasma frequency ω p e 2 = n q 2 / m e ε 0.
Scaling velocities by the thermal electron velocity ( v = V v T e and v T β = R v β v T e) and wavenumber quantities by the electron Debye length ( k = k / λ D e, k p = k p / λ D e, and K = K / λ D e with λ D e = k B T e ε 0 / n q 2), we arrive at expressions for the drag and diffusion with dimensionless integrals
A e / β v = ω p e v T e 2 5 / 2 π 3 / 2 N D , mac 1 + m e m β 1 R v β 3 d k p S k S k p ϵ k · V , k K 2 2 K K · k p k · V k p 3 e k · V 2 / k p 2 R v β 2 ,
(A13)
D e / β v = ω p e v T e 2 2 5 / 2 π 3 / 2 N D , mac 1 R v β d k p S k S k p ϵ k · V , k K 2 2 K K k p e k · V 2 / k p 2 R v β 2 .
(A14)
The numerical dielectric function evaluated for a quasineutral plasma with Maxwellian electrons and singly charged ions at temperatures T e and T i, respectively, can be written as
D k · V , k 1 p S k p 2 K · k p 2 K 2 k p 2 Z k · V k p + T e T i Z k · V k p R v i .
(A15)
Here, Z x is the derivative of the plasma dispersion function evaluated at x, which was calculated using PlasmaPy.40 We chose to use this full Balescu–Lenard collision operator analog rather than simplifying the dielectric function to obtain a Landau collision operator analog for the numerical collision operator based on the findings of Okuda and Birdsall who demonstrated that, unlike with point charges in three dimensions, the full collision operator is not well approximated by a Landau collision operator analog in a reduced number of dimensions or with particles of finite width.17 

The drag and diffusion coefficients were found to have little dependence on the orientation of the test particle velocity with respect to the grid, and so for most of the following calculations, V was chosen to point along the k x axis. This would clearly not be the case for unequal grid spacing, in which case we would recommend using the smaller spacing size to obtain a larger estimate of the drag and diffusion coefficients scaled by macroparticle density. The number of particles per cell can then be adjusted to account for the larger cell size.

Defining a grid size R g = Δ x / λ D e, the k x and k y integrals were evaluated numerically at velocities up to 3  v T e. We would expect the grid aliasing effects to be most pronounced at larger grid sizes, so a comparison of the R g = 1 case calculated with the first aliasing modes included p x , p y [ 1 , 1 ] was compared to the case with only the zeroth mode included p x = p y = 0. For both the spectral Maxwell solver and the FDTD Maxwell solver, under the grid scalings considered, the aliasing modes have a negligible effect on the drag and diffusion coefficients.

1.
Z.
Donko
,
A.
Derzsi
,
M.
Vass
,
B.
Horvath
,
S.
Wilczek
,
B.
Hartmann
, and
P.
Hartmann
, “
eduPIC: An introductory particle based code for radio-frequency plasma simulation
,”
Plasma Sources Sci. Technol.
30
,
095017
(
2021
).
2.
S.
Rauf
,
D.
Sydorenko
,
S.
Jubin
,
W.
Villafana
,
S.
Ethier
,
A.
Khrabrov
, and
I.
Kaganovich
, “
Particle-in-cell modeling of electron beam generated plasma
,”
Plasma Sources Sci. Technol.
32
(
5
),
055009
(
2023
).
3.
M. M.
Turner
,
A.
Derzsi
,
Z.
Donko
,
D.
Eremin
,
S. J.
Kelly
,
T.
Lafleur
, and
T.
Mussenbrock
, “
Simulation benchmarks for low-pressure plasmas: Capacitive discharges
,”
Phys. Plasmas
20
,
013507
(
2013
).
4.
M.
Vass
,
P.
Palla
, and
P.
Hartman
, “
Revisiting the numerical stability/accuracy conditions of explicit PIC/MCC simulations of low-temperature gas discharges
,”
Plasma Sources Sci. Technol.
31
,
064001
(
2022
).
5.
D.-Q.
Wen
,
J.
Krek
,
J. T.
Gudmundsson
,
E.
Kawamura
,
M. A.
Lieberman
, and
J. P.
Verboncoeur
, “
Particle-in-cell simulations with fluid metastable atoms in capacitive argon discharges: Electron elastic scattering and plasma density profile transition
,”
IEEE Trans. Plasma Sci.
50
(
9
),
2548
2557
(
2022
).
6.
O. C.
Eldridge
and
M.
Feix
, “
Numerical experiments with a plasma model
,”
Phys. Fluids
6
,
398
(
1963
).
7.
J.-B.
Fouvry
,
B.
Bar-Or
, and
P.-H.
Chavanis
, “
Kinetic theory of one-dimensional homogeneous long-range interacting systems sourced by 1/N2 effects
,”
Phys. Rev. E
100
,
052142
(
2019
).
8.
R. W.
Hockney
, “
Measurements of collision and heating times in a two-dimensional thermal computer plasma
,”
J. Comput. Phys.
8
,
19
44
(
1971
).
9.
D.
Montgomery
and
C. W.
Nielson
, “
Thermal relaxation in one- and two-dimensional plasma models
,”
Phys. Fluids
13
,
1405
(
1970
).
10.
D.
Sydorenko
,
A.
Khrabrov
,
W.
Villafana
,
S.
Ethier
, and
S.
Janhunen
, see https://github.com/PrincetonUniversity/EDIPIC-2D for “
EDIPIC-2D online repository
” (
2022
).
11.
S. P.
Gary
,
Y.
Zhao
,
R. S.
Hughes
,
J.
Wang
, and
T. N.
Parashar
, “
Species entropies in the kinetic range of collisionless plasma turbulence: Particle-in-cell simulations
,”
Astrophys. J.
859
,
110
(
2018
).
12.
C.
Birdsall
and
A.
Langdon
,
Plasma Physics via Computer Simulation
(
Taylor and Francis
,
2004
).
13.
M.
Touati
,
R.
Codur
,
F.
Tsung
,
V. K.
Decyk
,
W. B.
Mori
, and
L. O.
Silva
, “
Kinetic theory of particle-in-cell simulation plasma and the ensemble averaging technique
,”
Plasma Phys. Controlled Fusion
64
,
115014
(
2022
).
14.
M. M.
Turner
, “
Kinetic properties of particle-in-cell simulations compromised by Monte Carlo collisions
,”
Phys. Plasmas
13
,
033506
(
2006
).
15.
P. Y.
Lai
,
T. Y.
Lin
,
Y. R.
Lin-Liu
, and
S. H.
Chen
, “
Numerical thermalization in particle-in-cell simulations with Monte-Carlo collisions
,”
Phys. Plasmas
21
,
122111
(
2014
).
16.
P. Y.
Lai
,
L.
Chen
,
Y. R.
Lin-Liu
, and
S. H.
Chen
, “
Study of discrete-particle effects in a one-dimensional plasma simulation with the Krook type collision model
,”
Phys. Plasmas
22
,
092127
(
2015
).
17.
H.
Okuda
and
C. K.
Birdsall
, “
Collisions in a plasma of finite-size particles
,”
Phys. Fluids
13
,
2123
(
1970
).
18.
B.
Abraham-Shrauner
, “
Test particle in a two-dimensional plasma
,”
Physica
43
,
95
104
(
1969
).
19.
M. A.
Reynolds
,
B. D.
Fried
, and
G. J.
Morales
, “
Velocity-space drag and diffusion in a model, two-dimensional plasma
,”
Phys. Plasmas
4
(
5
),
1286
1296
(
1997
).
20.
J. M.
Dawson
, “
Thermal relaxation in a one-species, one-dimensional plasma
,”
Phys. Fluids
7
,
419
(
1964
).
21.
J.
Hsu
,
G.
Joyce
, and
D.
Montgomery
, “
Thermal relaxation of a two-dimensional plasma in a d.c. magnetic field. Part 2. Numerical simulation
,”
J. Plasma Phys.
12
,
27
31
(
1974
).
22.
I.
Jechart
,
T.
Katsouleas
, and
J.
Dawson
, “
Anomalous thermal relaxation of a two-dimensional magnetized plasma
,”
Phys. Fluids
30
(
1
),
65
(
1987
).
23.
J.
Virtamo
and
H.
Tuomisto
, “
Verification of a simple collision operator for one-dimensional plasma by simulation experiments
,”
Phys. Fluids
22
,
172
(
1979
).
24.
J.
Dawson
, “
One-dimensional plasma model
,”
Phys. Fluids
5
,
445
(
1962
).
25.
J.-Y.
Hsu
,
D.
Montgomery
, and
G.
Joyce
, “
Thermal relaxation of a two-dimensional plasma in a d.c. magnetic field. Part 1. Theory
,”
J. Plasma Phys.
12
,
21
26
(
1974
).
26.
N.
Gatsonis
and
A.
Spirkin
, “
A three-dimensional electrostatic particle-in-cell methodology on unstructured Delaunay-Voronoi grids
,”
J. Comput. Phys.
228
,
3742
3761
(
2009
).
27.
S.
Averkin
and
N.
Gatsonis
, “
A parallel electrostatic Particle-in-Cell method on unstructured tetrahedral grids for large-scale bounded collisionless plasma simulations
,”
J. Comput. Phys.
363
,
178
199
(
2018
).
28.
A. T.
Powis
and
I. D.
Kaganovich
, “
Accuracy of the explicit energy-conserving particle-in-cell method for under-resolved simulations of capacitively coupled plasma discharges
,” arXiv:2308.13092 (
2023
).
29.
H.
Sun
,
S.
Banarjee
,
S.
Sharma
,
A. T.
Powis
,
A. V.
Khrabrov
,
D.
Sydorenko
,
J.
Chen
, and
I. D.
Kaganovich
, “
Direct implicit and explicit energy-conserving particle-in-cell methods for modeling of capacitively-coupled plasma devices
,”
Phys. Plasmas
30
,
103509
(
2023
).
30.
S. G.
Walton
,
D. R.
Boris
,
S. C.
Hernandez
,
E. H.
Lock
,
T. B.
Petrova
,
G. M.
Petrov
,
A. V.
Jagtiani
,
S. U.
Engelmann
,
H.
Miyazoe
, and
E. A.
Joseph
, “
Electron beam generated plasmas: Characteristics and etching of silicon nitride
,”
Microelectron. Eng.
168
,
89
96
(
2017
).
31.
S. G.
Walton
,
D. R.
Boris
,
S. G.
Rosenberg
,
H.
Miyazoe
,
E. A.
Joseph
, and
S. U.
Engelmann
, “
Etching with electron beam-generated plasmas: Selectivity versus ion energy in silicon-based films
,”
J. Vac. Sci. Technol., A
39
,
033002
(
2021
).
32.
K.
Nanbu
, “
Theory of cumulative small-angle collisions in plasmas
,”
Phys. Rev. E
55
,
4642
4652
(
1997
).
33.
T.
Charoy
,
J.-P.
Boeuf
,
A.
Bourdon
,
J.
Carlsson
,
P.
Chabert
,
B.
Cuenot
,
D.
Eremin
,
L.
Garrigues
,
K.
Hara
,
I. D.
Kaganovich
,
A. T.
Powis
,
A.
Smolyakov
,
D.
Sydorenko
,
A.
Tavant
,
O.
Vermorel
, and
W.
Villafana
, “
2D axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas
,”
Plasma Sources Sci. Technol.
28
,
105010
(
2019
).
34.
V. A.
Godyak
,
R. B.
Piejak
, and
B. M.
Alexandrovich
, “
Probe diagnostics of non-Maxwellian plasmas
,”
J. Appl. Phys.
73
,
3657
3663
(
1993
).
35.
Y.
He
,
Y.-M.
Lim
,
J.-H.
Lee
,
J.-H.
Kim
,
M.-Y.
Lee
, and
C.-W.
Chung
, “
Effect of parallel resonance on the electron energy distribution function in a 60 MHz capacitively coupled plasma
,”
Plasma Sci. Technol.
25
,
045401
(
2023
).
36.
I. D.
Kaganovich
and
L. D.
Tsendin
, “
The space-time-averaging procedure and modeling of the RF discharge, Part II: Model of collisional low-pressure RF discharge
,”
IEEE Trans. Plasma Sci.
20
(
2
),
66
75
(
1992
).
37.
S. V.
Berezhnoi
,
I. D.
Kaganovich
, and
L. D.
Tsendin
, “
Generation of cold electrons in a low-pressure RF capacitive discharge as an analogue of a thermal explosion
,”
Plasma Phys. Rep.
24
,
556
563
(
1998
).
38.
K.
Yee
, “
Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media
,”
IEEE Trans. Antennas Propag.
14
(
3
),
302
307
(
1966
).
39.
J.
Villasenor
and
O.
Buneman
, “
Rigorous charge conservation for local electromagnetic field solvers
,”
Comput. Phys. Commun.
69
,
306
316
(
1992
).
40.
PlasmaPy Community,
PlasmaPy, Version 2023.5.1
(
Zenodo
,
2023
).