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.

## I. INTRODUCTION

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.

^{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 \lambda D$ in a 1D simulation, $ N Dmac= n mac \lambda D 2$ in a 2D simulation, and $ N Dmac= n mac \lambda D 3$ in a 3D simulation with $ \lambda 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}

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.

## II. NUMERICAL THERMALIZATION DRAG AND DIFFUSION COEFFICIENTS IN 2D PIC

*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

As in a real plasma, the contribution to the dielectric function $\u03f5 k \xb7 V , k$ from each species $\gamma $ is proportional to its plasma frequency squared: $ \omega p \gamma 2= n \gamma q \gamma 2 / \epsilon 0 m \gamma $. Also the dielectric function is the plasma dispersion function $Z(x)$, where $Z\u2032(x)$ is its derivative evaluated at $x$. We choose not to approximate $\u03f5 k \xb7 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 \lambda 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.

Here, $ V e=v/ v Te$ is the velocity scaled by the electron thermal velocity, $ N D=n \lambda D 3$ is the number of real particles per Debye cube, and $ ln \Lambda e e= ln 4 \pi 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-Shrauner^{18} for single-species point particles in 2D with a Lorentz background distribution, those derived by Okuda and Birdsall^{17} 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\u226a 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\u226b 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\u226b v T$.^{19}

## III. NUMERICAL THERMALIZATION RATES AND MITIGATION STRATEGIES

### A. Relaxation processes and rates

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 ( $ \tau 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 plasma^{24} due to the necessity of distinguishing between drag and diffusion timescales for a test particle ( $ \tau s, \tau \u2225, \tau \u22a5\u223c N Dmac$) and processes, which evolve the total velocity distribution function ( $ \tau R\u223c 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 ( $ \tau s, \tau R\u223c\u2009 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 $ \tau R\u223c N Dmac.$^{21,22,25}

^{8}empirical estimate of the same

We should note that Hockney's method^{8} 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 Spirkin^{26} and work by Averkin and Gatsonis.^{27}

### B. Mitigation strategies

Given that numerical thermalization happens on timescales roughly comparable to $ N Dmac \tau 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 schemes^{28} 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 $ \omega pedt$ and $ \Delta x/ \lambda D e$ is roughly 0.1. As can be seen in Fig. 4(b), the direct implicit simulation with a grid spacing of $ \Delta x=3 \lambda 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.

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 ( $ \Delta x=0.2 \lambda D$) contain 40 particles per cell, while the direct implicit case with an unresolved Debye length ( $ \Delta x=3 \lambda 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, $ \epsilon 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 $ \epsilon 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\u2248 2 N D / ln \Lambda 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.

## IV. EXAMPLE APPLICATIONS TO 2D PIC SIMULATIONS OF LOW-TEMPERATURE PLASMA SOURCES

### A. Electron beam-generated plasma in a magnetic field

Electron beam-generated plasmas are of interest as a source of plasmas with large densities ( $ 10 15$– $ 10 17\u2009 m \u2212 3$) and low-electron temperatures ( $\u22721\u2009eV$) suitable for materials processing applications.^{2,30,31} This can produce a large flux of low-energy ( $<5\u2009eV$) 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-2D^{10} 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.

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 \xaf 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\u2009\mu 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\u2009ns$), 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 \u2212 Ar=6.75\u2009ns$ to $ t e \u2212 Ar=180\u2009ns$. Here, $ t e \u2212 Ar=1/ \nu e \u2212 Ar$ is the inverse of the electron elastic collision rate, where $ \nu e \u2212 Ar= n Ar \sigma e l E 2 E / m e$ depends on the background gas density $ n Ar$ and the collision cross section $ \sigma 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\xd7 10 15\u2009 m \u2212 3$, the electron Debye length near the beam is 0.19 mm and the plasma frequency is $2.5\xd7 10 9\u2009 s \u2212 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 $ \tau R num\u2248$ 1.1 $\mu 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.

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

### B. Hollow cathode plasma

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 $\mu s$ is shown in Fig. 6(b), averaged over the channel area shown in Fig. 6(a).

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 \xaf 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\u2009\mu 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\u2009ns$. The elastic collision timescale for a thermal electron at 2.7 eV in this simulation is $ t e \u2212 A r=9\u2009ns$, while the same for 35 eV electron is only $ t e \u2212 Ar=1.5\u2009ns$.

#### 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\u2009 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.

### C. Inductively coupled plasma

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 $\mu s$ before the snapshots depicted in Fig. 8 were taken.

#### 1. Timescale analysis: Inductively coupled plasma

### D. Capacitively coupled plasma

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 $\mu s$ of a 295- $\mu 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 collisionless^{33} and collisional regime.^{3}

^{34}[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}

#### 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\u2009ns$. The electron–neutral elastic collision time for an electron with 3.5 eV of kinetic energy is about $ t e \u2212 Ar=56\u2009ns$. One oscillation of the applied voltage occurs over a time period of about $ \tau RF=460\u2009ns$.

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.

## V. CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX: CALCULATION OF THE DRAG AND DIFFUSION COEFFICIENTS OF THE NUMERICAL COLLISION OPERATOR

*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

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

*et al.*give the following expressions for its

*j*th component when a spectral Maxwell solver or alternately the second-order finite-difference time-domain (FDTD) method of Yee

^{38}coupled with the charge conserving scheme of Villasenor and Buneman

^{39}is implemented

^{17–19}However, while $ A \alpha 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

^{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= \Delta x/ \lambda 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\u2208[\u22121,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.

## REFERENCES

_{2}effects

*Plasma Physics via Computer Simulation*

*PlasmaPy, Version 2023.5.1*