Direct implicit and explicit energy-conserving particle-in-cell methods for modeling of capacitively coupled plasma devices

Simulation and modeling of plasma devices play an important role in the development and optimization of modern plasma reactors due to the high cost of experimentation over a wide range of operating conditions.^1,2 Traditional methods for modeling plasma devices include the fluid approach,^3–10 hybrid fluid/kinetic models,¹¹ and kinetic methods, such as the particle-in-cell (PIC) technique.^12,13 In the fluid approach, fluid equations are solved for all plasma species under the assumption that the particle velocity distribution function (VDF) is Maxwellian. The greatest advantage of this approach is a smaller computational cost comparing with kinetic simulations,^4,8,14,15 nominally because the assumed VDF reduces the number of equations to solve. However, important kinetic effects are neglected in this approach, when the VDF cannot be the presumed, and therefore, the accuracy of such a model is questionable in industrial processing devices operating at low pressures (where the mean free path is on the order of or larger than the sheath size).^16–18 PIC methods are often used to study the kinetic effects by solving the Vlasov equation,^19–29 with the most common method being the traditional explicit momentum-conserving algorithm. For the electrostatic explicit momentum-conserving algorithm, the time step is limited by the fastest electron dynamics (⁠ $ω p edt≤0.2$⁠), while the cell size is limited by the Debye radius (⁠ $dx≤ λ D e$⁠). Here, $ω p e( = n 0 e 2 / ε 0 m e)$ and $λ D e( = ε 0 T e / n 0 e)$ are the electron plasma frequency and Debye radius; $dt$ and $dx$ are the time step and cell size, respectively. The requirement for such a fine resolution results in very high computational cost.

The current trend in plasma processing is to operate at lower gas pressures (⁠ $∼$ several mTorr)^30–34 to achieve atomic-scale precision. The capacitively coupled plasma (CCP) reactor, which operates in the MHz range of radio-frequencies (rf),^35,36 is the most widely used contemporary plasma device for etching of silicon wafers, critical to the semiconductor industry. There have been a large number of modeling and experimental work exploring various aspects of CCP discharges. In single-frequency capacitively coupled plasma discharges, the ion flux and ion energy cannot be controlled independently since they depend simultaneously the input power for a fixed neutral gas pressure and system geometry.^31,37–46 Several other approaches have been suggested to circumvent this restriction, such as dual-frequency capacitively coupled plasma (DF-CCP) discharges,^47–54 electrical asymmetric effects,^55–59 and non-sinusoidal tailored voltage/current waveform excitations.^60–68 There are also several fascinating physical phenomena in very high frequency (VHF) driven CCP discharges, such as the generation of an electron beam accelerated by the sheath electric fields^69–71 as well as the production of higher harmonics in voltage and current.^72–79 

These devices require modeling of plasma chemistry with numerous plasma and radical species, meaning that even 1D-3V PIC simulations can become computationally expensive. Modern plasma devices, on the other hand, often operate with plasma density up to the order of $10 17 m − 3$⁠, and the size of a typical CCP discharge chamber is at least $∼10 cm$⁠. For realistic 2D simulations of a plasma discharge, this requires $∼ 10 6$ cells and usually $∼ 10 7$ time steps to achieve convergence. Furthermore, in order to reduce numerical noise, the number of macro-particles in each cell should be higher than 100,^80,81 requiring hundreds of millions or billions of simulation particles evolved over tens of millions of time steps. To achieve reasonable simulation time, one faces a difficult dilemma, whether to reduce the simulation size of the modeled system, which could alter some plasma properties, or spend significant computational time, often weeks or sometimes even months, to achieve a steady state solution. Therefore, it is essential to remove the grid size and time step restrictions on the PIC method if we are to model the whole plasma device and reduce the actual simulation run time, so that these new CCP reactors can be modeled for prototyping.

The direct implicit particle-in-cell (DI-PIC) method was developed in 1980s to provide such a reduction in computational cost, with the added benefit of stability even when under-resolving the Debye length and plasma frequency,^82–85 allowing for a larger time step and cell size. Within a CCP discharge resolution of the fast plasma electron dynamics with timescale $t∼ ω p e − 1$ is not crucial, meaning that accuracy can be maintained even when the electron time and length scales are under-resolved.

Similarly, the explicit energy conserving particle-in-cell (EC-PIC) algorithm proposed by Lewis⁸⁶ alleviates the constraint on the cell size (⁠ $dx$⁠) by making small modifications to the traditional PIC method. Despite several efforts to study the energy conserving and direct implicit methods, most studies have focused on uniform plasmas (using periodic boundary conditions). Other early DI-PIC modeling efforts were limited to 1D simulations of rf discharges.^87,88 Applying the DI-PIC approach in two-dimensions requires a more complex field solver and particle pusher. Other advanced implicit algorithms^85,89–91 have also been proposed. However, these approaches are more computationally expensive per time step and are complex to implement when compared to the much simpler DI-PIC and EC-PIC studied in this paper.

DI-PIC and EC-PIC are implemented into the EDIPIC-2D code⁹² and EC-PIC is implemented into the LTP-PIC code⁹³ (see a full description in Sec. III D), and their accuracy and performance are tested through simulations of a 2D capacitively coupled plasma discharge. It is found that both approaches allow for cell sizes larger than the Debye length, and a corresponding reduction in runtime, with only a small penalty in accuracy. The direct implicit method can also allow one for time steps larger than the electron plasma period; however, this can lead to numerical heating or cooling. In these cases, the tuning of the cell size to time step ratio can reduce this heating/cooling to tolerable levels.

To model the collision operator, we use the well-known null-collision algorithm for Monte Carlo collisions (MCCs) and approach which has been proven to be accurate and efficient at reducing numerical cost when modeling many applications of low-pressure devices.^21,94,95

We performed a large number of CCP simulations with both the DI and EC methods, as well as benchmark them against results using explicit momentum conserving PIC simulations. We determined that the computational efficiency is significantly improved while still maintaining sufficient accuracy. For DI-PIC, the accuracy can be further improved if the initial time step and cell size are chosen according to an “optimum path” recipe to be described further in the paper. An estimate of numerical heating in the direct implicit method is also provided, which enables trustworthy modeling of large plasma devices.

The paper is organized as follows: The simulation algorithm and model are introduced in Sec. II. Code tests and benchmarks are presented in Sec. III. Finally, the discussion of results and conclusions are given in Sec. IV.

The direct implicit (DI) particle-in-cell method predicts the charge density at the next time step using a linear approximation of the new electric field.⁹⁶ This approach yields an elliptic equation for the potential whose coefficients depend directly on particle data accumulated on the spatial grid in the form of an effective linear susceptibility. The rank of the matrix equation is determined by the number of field quantities defined on the spatial grids, so that the potential profile can be solved using standard numerical linear algebra techniques. The particle update is calculated for each particle concurrently.^83,85 The DI time advance loop is comprised of a predictive push for the particles, followed by solving the electrostatic field equation and then by a second, final particle push. In the presence of an external magnetic field, the electrostatic field equation takes the form of a Poisson equation for an anisotropic dielectric. On a Cartesian grid, the corresponding finite-difference approximation involves a nine-point “box” type stencil, compared to the five-point “star” appearing in the explicit scheme.⁹⁷ Also, since the effective dielectric tensor in the DI scheme is a function of the local charge density, the matrix supplied to the field solver must be re-initialized at each time step. In the un-magnetized case, the box stencil reduces to the five-point “star” stencil. A flow chart for the DI-PIC method is shown in Fig. 1.

Note that $q=−e$ in the case of electron species to which the DI treatment is being applied. In the expressions above, $Ω$ is the particle cyclotron frequency, $θ$ is a dimensionless notion of $Ω$ in the algorithm, $A − 1$ and $ϵ α β$ are defined as the field matrix, $δ α β$ is the unitary tensor, $e α β γ$ is the anti-symmetric tensor, and $α, β, γ$ denote different spatial dimensions.

Equation (8) is discretized with the help of Gauss's law with the finite-difference stencil shown in Fig. 2. It is convenient to identify the coefficients of the linear equation arising in the differencing scheme according to the points of the compass. The four-unit normal vectors involved in evaluating the flux of the $D$ vectors are pointing North (N), South (S), East (E), and West (W). The Northwest (NW), Northeast (NE), Southeast (SE), and Southwest (SW) directions are also shown in the figure. The $x$ and $y$ components of the electric field $E$ are calculated at the half-nodes through the second-order central differencing scheme.

Then, the matrix coefficients contributed by the given locally numbered $(i,j)$ node where the charge density is specified on the right-hand side are given in Table I.

The resulting properties of the DI permittivity tensor were taken into account to calculate the respective elements of the matrix fed to the field solver.

Note that in the DI scheme, the $K$ matrix for the final push is evaluated at the particle position corresponding to the previous time step.⁸³ In order to evaluate $K n$⁠, the electron particle positions are back-tracked in order to reduce memory consumption.

One of the important properties that limits the performance of the DI algorithm is the presence of numerical heating or cooling. This arises since the DI algorithm does not guarantee exact energy conservation, which is due to the dissipative nature of the implicit time integration scheme.^83,84 However, previous works has demonstrated that these effects can nearly vanish when the time step (⁠ $dt$⁠) and cell size (⁠ $dx$⁠) follow an “optimal path,” which is a linear relation between time step and cell size.^82–84 Therefore, the numerical heating/cooling of the DI algorithm is controllable provided we choose the time step and cell size appropriately. Here, we verify this important property using our in-house 2D-3V code by simulating a double periodic system with initial plasma density $n e 0=1.1× 10 14 m − 3$ and initial temperature $T e 0=1 eV$⁠. A double periodic system typically refers to a 2D system where all the four boundaries are periodic, where particles moving out from one boundary are added to the other side with the same velocity. This system also allows for the Poisson equation to be solved via Fourier techniques. In each cell, we have $N=100$ macro-particles. We alter the cell size and time step and compare how the numerical heating/cooling affects the temperature of the electrons. The normalized numerical heating $ΔE/ E 0N$ for different cases at $t=591.2 ω p e − 1$ is shown in Table II. To show the numerical heating more clearly, Fig. 3(a) shows the numerical heating as a function of cell size, different lines are for different time steps. When the time step is increased, a near-zero numerical heating can be achieved if we increase cell size simultaneously. Figure 3(b) plots the data in the table, where a much smaller numerical heating can be clearly seen, illustrated by the black region.

As can be seen from Table II and Fig. 3, a near-zero numerical heating is possible when the ratio between $ω p edt$ and $dx/ λ D e$ is approximately $0.1$⁠. We refer to this as the “optimal path” of the DI algorithm, since using these time steps and cell sizes provides near energy-conservation, which is beneficial for long timescale modeling of plasma devices.

The explicit energy conserving PIC method (EC-PIC) here refers to the scheme originally reported by Lewis in 1970,⁸⁶ whereby a variation approach is adopted for deriving the appropriate interpolation functions and Poisson equation/electric field finite difference discretizations so as to conserve system energy in the limit $dt→0$⁠. In practice, energy remains very well-conserved when adhering to the standard MC-PIC time constraint (⁠ $ω p edt<0.2$⁠), making the scheme resilient to the finite-grid-instability, even for $dx≫ λ D e.$⁸⁶ EC-PIC schemes for both electrostatic and electromagnetic PIC codes have been derived in Refs. 99 and 100 and placed on a more sound theoretical footing through the language of discrete exterior calculus.¹⁰¹ Initial investigations into the scheme revealed a deleterious cold-beam, finite-grid instability which quickly resulted in unphysical behavior.¹⁰² However, Ref. 97 demonstrated that explicit EC-PIC could safely be applied to warm plasmas primarily driven by ambipolar electric fields, such as the system investigated in this paper. Therefore, the EC-PIC method could prove to be a useful method to model large scale plasma systems without the need to resolve the Debye length, an issue this paper will aim to explore.

The EC-PIC scheme is very similar to the traditional MC-PIC approach with a few important exceptions. A flow chart for EC-PIC scheme is shown in Fig. 4.

The algorithm proceeds as follows:

1. Particles are updated via the standard Boris algorithm.¹⁰⁰ 

2. Particle charge density is interpolated to the grid via multi-linear (in this instance bi-linear) or cloud-in-cell interpolation.

3. The Poisson equation in EC-PIC is solved using exactly the same method as in explicit MC-PIC code, where the standard a five-point stencil is implemented.

4. Modification: The electric field is computed at the edge centers between cell nodes via the central difference scheme. For example, consider discrete electric potential $ϕ i , j$ on a 2D lattice with $i$ and $j$ representing the discrete $x$ and $y$ coordinates, respectively. The electric field components are computed as

   $E x , i + 1 / 2 , j=− ϕ i + 1 , j − ϕ i , j d x,$

   (15)

   $E y , i , j + 1 / 2=− ϕ i , j + 1 − ϕ i , j d y.$

   (16)
5. Modification: Interpolation of the electric field to the particles is performed using 0th order (nearest grid point) in the same direction as the electric field component and first order (linear) in all other directions. For a particle $p$ located within cell $i,j$ with node $i,j$ at the South–West corner, the electric field components are defined as

   $E x , p= ∑ i , j E x , i + 1 / 2 , j S 0 ( x p − x i + d x / 2 ) S 1( y p − y j),$

   (17)

   $E y , p= ∑ i , j E y , i , j + 1 / 2 S 1 ( x p − x i ) S 0 ( y p − y j + d y / 2 ),$

   (18)
   where

   $S 0( x)={ 1 , − d x / 2 ≤ x ≤ d x / 2 , 0 , otherwise ,$

   (19)

   $S 1( x)={ 1 + x / d x , − d x ≤ x ≤ 0 , 1 − x / d x , 0 ≤ x ≤ d x .$

   (20)

   Here, $S 0$ and $S 1$ are the zero order and first order particle shape functions. Note that other shape functions can also be used to provide higher order interpolation.⁹⁰ 

To perform code testing and benchmarking, we model a two-dimensional (2D) capacitively coupled (CCP) discharge using the traditional explicit momentum conserving PIC (MC-PIC), and our newly developed direct implicit algorithm (DI-PIC) and explicit energy conserving PIC (EC-PIC) in EDIPIC-2D. The EC-PIC results are further benchmarked our in-house low-temperature plasma particle-in-cell (LTP-PIC) code. A schematic diagram of the simulation domain is shown in Fig. 5, with $L x× L y=108×36 mm$². A powered electrode with voltage amplitude $V r f=100 V$ and $f r f=27.12 MHz$ is placed at the bottom edge, while the other three boundaries are grounded electrodes. A small vacuum gap of 2.25 mm is implemented between the powered electrode and the grounded electrodes. The initial plasma temperature and density are $T e 0=2 eV$ and $n e 0=5× 10 15 m − 3$⁠, respectively. Argon is used as the working neutral gas, with pressure $P g=5 mTorr$ for all simulations. We include electron-neutral elastic, excitation (at 11.15 eV only) and ionization collisions as well as ion-neutral charge exchange. Tables for electron-neutral cross sections are from Refs. 103 and 104 and the semi-empirical model from Ref. 105 is used to compute ion charge-exchange cross sections. We initialize all simulations with 100 macro-particles per-cell for all cases and run for $102 μs$⁠, corresponding to more than $2500$ rf cycles, at which point steady state is reached. We performed two simulations using MC-PIC with $d x explicit=7.5× 10 − 5 m≈ λ D e , 0/2$⁠, $d t explicit= 10 − 11 s$ as wall as with $d x explicit=1.5× 10 − 4 m≈ λ D e , 0$⁠, $d t explicit= 10 − 11 s$ to ensure numerical convergence. We use the most highly resolved of these cases as a reference to illustrate the accuracy of explicit EC-PIC and implicit DI-PIC methods (this reference case is referred to as “explicit” in figure captions), when increasing systematically both the time step and cell size, such that $dx$ may be greater than the Debye length $λ D e$ and time step larger than the electron plasma frequency, i.e., $ω p e − 1$⁠. All simulation parameters are shown in Table III, where we set $d x 0= λ D e , 0$⁠, $d t 0=d t explicit$ as the reference cell size and the time step, respectively.

Both the direct implicit (DI) and energy conserving (EC) algorithms can offer improvements in simulation time. We show a comparison in computational hours between DI-PIC simulations and MC-PIC simulations in Fig. 6. Each of these simulations is performed on the ANTYA HPC facility at the Institute for Plasma Research in India with 1440 cores. As we can see, when comparing explicit MC-PIC, DI-PIC, and EC-PIC with the same time step and cell size, the computational cost is similar. Improvements become apparent when considering larger cell sizes for EC-PIC and DI-PIC, where adopting cell sizes $d x / d x 0 > 1$ can allow for order of magnitude reductions in compute time. This would not be possible with MC-PIC due to the emergence of the finite-grid instability under these conditions. DI-PIC also allows for adopting a larger time step $d t / d t 0 > 1$⁠, allowing for more than an order of magnitude reduction in runtime. Therefore, using larger time step (⁠ $dt$⁠) and cell size (⁠ $dx$⁠) in DI-PIC and EC-PIC allows for significantly faster simulations, which is of benefit to future industrial applications. However, the computational efficiency does not increase proportionally with the increase in time step and cell size. This is mainly due to the equal need for outputting diagnostics, which impose some computational time nearly independent of time step and cell size.

To demonstrate the accuracy of EC-PIC and DI-PIC, we compare important physical parameters with the reference simulations performed with MC-PIC. Figures 7(a)–7(d) show profiles of time averaged electron density, electron temperature, electron power absorption, and ion current, respectively, of the 2D reference simulations. All data are averaged over more than 30 rf cycles. The electron density profile (⁠ $〈 N e〉$⁠) is peaked at the center (⁠ $∼1.6× 10 15 m − 3$⁠) although the electron temperature (⁠ $〈 T e〉$⁠) is low at the center of discharge chamber (⁠ $∼2 eV$⁠) and increases toward the boundaries. Temperature peaks in the vicinity of the sheath region (⁠ $∼3.3 eV$⁠) and further decreases toward the wall. The electron power absorption (⁠ $〈 E . J e〉$⁠) is higher near to the powered electrode (bottom edge) compared to grounded electrode (top and side edges) and is asymmetric in nature [see Fig. 7(c)]. Figure 7(d) clearly shows that the ion current is increasing toward the bottom electrode (i.e., powered) and the top electrode (i.e., grounded). In the following figures, we compare the algorithms via 1D slices, with all data taking along the vertical centerline (from powered to grounded electrode) as shown in Fig. 7(a).

Figure 8 shows the comparison of time averaged ion current for different cases. In Figs. 8(a)–8(c), for DI-PIC, the values of $dx$ is varied from $d x 0$⁠, $2d x 0$⁠, and $3d x 0$⁠, respectively, for specific values of $dt$ (i.e., $d t 0$⁠, $2d t 0$⁠, $5d t 0$⁠, and $10d t 0$⁠). It is important to note that the instances involving 0.5  $d x 0$ for both DI-PIC and EC-PIC are excluded from display because, although they offer greater accuracy, they are not computationally more efficient. Figures 8(d)–8(f) further show the time averaged ion current for EC-PIC cases for the same values of $Δx$ and $Δt$ mentioned above. We found excellent agreement comparing the current at the cathode in DI-PIC and EC-PIC with the reference MC-PIC case, with an average deviation of about $7.7%$ (all simulation cases, including DI-PIC and EC-PIC). This indicates that our algorithms could correctly capture the profiles of ion flux, highly relevant to the performance in CCPs used for silicon etching. This could greatly enhance and impact the modeling of large volume plasma devices in the future. Figure 9 presents the comparison of time averaged energy deposition to the electrons $〈 E ⋅ J e 〉 t$⁠. The integrated energy deposition in the plot for the explicit case is $〈 E ⋅ J e 〉 t=1.92× 10 4 W/ m 2$ (integrated over the 1D lineout). The average deviation of DI-PIC and EC-PIC is 21.1% and 12.9%, respectively. Almost all EC-PIC cases [i.e. (d)–(f)] provide reasonable agreement with the MC-PIC (shown in blue), while DI-PIC cases [i.e. (a)–(c)] with large time steps show poorer performance. This is because at larger time steps the DI-PIC algorithm poorly resolves the timescale of the sheath motion.

Figure 10 gives a comparison of time averaged steady state electron density for different cases, i.e., for DI-PIC [(a)–(c)] and for EC-PIC [(d)–(f)]. As we can see, almost all the cases provide good agreement with the reference MC-PIC case, with an average deviation of electron density at the center of simulation domain of only $10.3%$ and $8.7%$ for DI-PIC and EC-PIC, respectively, strongly indicating the reliability of these methods. Figure 11 further shows the comparison of time averaged electron temperature. Figures 11(a)–11(c) are for the DI-PIC cases and (d)–(f) are for EC-PIC cases. The temperature at center of discharge for the explicit case is nearly $2 eV$⁠. In this case, the temperature appears to deviate more markedly, the average deviation at the center of simulation domain is $32.1%$⁠. However, excellent agreement is found for cases with $( d t , d x)=(2d t 0,d x 0)$⁠, i.e., in Fig. 11(a) for DI-PIC and (d) for EC-PIC. We speculate that the lower agreement for cases with higher $(dt,dx)$ is caused by artificial numerical heating in the sheath boundary.^106–108 On the other hand, however, previous work has shown that the deviation of temperature profile does not significantly affect the accuracy of the model when compared to experiments,^109,110 this is because the electron temperature profile depends sensitively on the details of cross-sectional description and dynamics of fast electrons, which may not be important for the slower ion dynamics. This explains why we have an excellent match of ion current profile despite a relatively poorer agreement in electron temperature.

The time evolution of different physical quantities is also compared via probe diagnostics. For each simulation, the probe is placed at the center of the simulation domain. Figures 12 and 13 show the probe diagnostics for ion density (⁠ $N i$⁠) and electron temperature (⁠ $T e$⁠), respectively. As we can see, the ion density (see Fig. 12) matches well with the explicit MC-PIC case (shown in blue), while time evolution of electron temperature (see Fig. 13) deviates, showing how numerical cooling or heating can play a role in determining the electron temperature in the simulations. For similar reasons, the deviation of temperature evolution is acceptable and may not affect plasma parameters relevant to real plasma devices.^109,110

The error analysis mentioned above is further summarized in Fig. 14, where not only the average error for each physical quantity but also the error in each simulation case is shown. Close to the optimal path for DI-PIC (corresponding to $dx/d x 0=1$⁠, $dt/d t 0=2$⁠), the error is significantly smaller. It is also shown that the error, in general, increases with the increase in time step and cell size; however, the increment is not significant. Therefore, our codes not only maintain stability under large time step and cell size but also have acceptable numerical accuracy, which strongly benefits the future industrial usage. Figure 15 shows the electron velocity distribution function (EVDF) for several representative simulations. We can see that all of them are non-Maxwellian, and the EVDFs for both DI-PIC and EC-PIC do not deviate significantly from the explicit case. Despite some small deviations in the distribution of energetic electrons shown in Fig. 15(b), the similarity in ion flux shown in Fig. 8 indicates that this has a small effect on the ions impacting the electrode. Therefore, although there is a small compromise of numerical accuracy within our code when modeling kinetic effects, this compromise is generally deemed acceptable. A more detailed analysis for such a numerical distortion will be investigated in future publications.

As was established in Sec. II, the numerical heating of DI has a linear dependence on the ratio of time step and cell size, i.e., $dt/dx$⁠. Such a dependence has only been tested in a double periodic system, not in a practical configuration with electrodes.⁸⁴ Here, in Fig. 16, we show that such a linear dependence persists when time the step is small [shown in Fig. 16(a)] but fails to hold when $dt>5d t 0$ [shown in Fig. 16(b)]. However, as shown in Fig. 16, the linear dependence persists when the rf voltage amplitude is increased to $V r f=1000 V$⁠, although analyzing such complex behavior of numerical heating is beyond the scope of this paper. Given a certain time step and cell size, one could estimate the numerical heating using a simple linear curve shown as shown in Figs. 16(a) and 17, with a caveat that this will no longer be valid for very large time steps. On the other hand, as shown in Fig. 6, a significant improvement in computational efficiency compared to traditional explicit MC-PIC is already achieved even when using a relatively small time step increase, say, $dt=2d t 0$⁠.

This implies that, despite the relative inaccuracy in estimating electron temperature compared with more accurate, but expensive, algorithms (for instance, fully implicit algorithms^90,91), the numerical heating of the direct implicit algorithm is acceptable and could be directly implemented into simulations for modeling practical plasma devices or be used to gain an insight into the physical trends before conducting full resolution simulations.

To gain more confidence in our results, we benchmarked our simulations for different cases against a 2D version of the LTP-PIC code. The low-temperature plasma particle-in-cell (LTP-PIC) code is a general-purpose high-performance 3D-3V kinetic plasma simulation software designed and developed at the Princeton Plasma Physics Laboratory (PPPL).^94,111,112 The energy conserving version of LTP-PIC (implementing the same algorithm as discussed in Sec. II B) can also simulate plasma devices with cell size much greater than the electron Debye length. At the same time, LTP-PIC uses the similar collision algorithms as EDIPIC-2D. Therefore, a direct comparison between simulation results from the two codes is possible.

Here, we consider four representative cases for the benchmark exercise, which are the explicit cases (⁠ $dt=0.5d t 0$ and $dx=0.5d x 0$⁠) (both EDIPIC-2D and LTP-PIC) and then using EC-PIC (EDIPIC-2D and LTP-PIC) and DI-PIC (EDIPIC-2D only) with $dt=5d t 0$ and $dx=1, 2, 3d x 0$⁠. Figures 18–21 show the profile benchmark results for time average ion current (⁠ $〈 J i y〉$⁠), time average energy deposition (⁠ $〈 E ⋅ J e〉$⁠), time average electron density (⁠ $〈 N e〉$⁠), and time average electron temperature (⁠ $〈 T e〉$⁠), respectively. For all four physical quantities, the deviations of explicit MC-PIC from LTP-PIC [shown by subfigures (a) in Figs. 18–21] are within $8%$⁠. For ion current at the cathode, DI-PIC and EC-PIC deviate from LTP-PIC by about $8.9%$⁠. For the sum of 1D energy deposition, DI-PIC and EC-PIC deviate from LTP-PIC by about $18.1%$⁠. For electron density at the simulation domain center, the average deviation of DI-PIC and EC-PIC from LTP-PIC is about $6.6%$⁠. For electron temperature at the domain center, the average deviation of DI-PIC and EC-PIC from LTP-PIC is about $6.9%$⁠. Therefore, all the cases exhibit a satisfactory level of agreement and deviations are believed to arise from subtle disparities in the Monte Carlo collision algorithms employed by each code.

Figure 22 further benchmarked the probe diagnostics for ion density (⁠ $N i$⁠), and excellent agreement is found. In conclusion, both DI-PIC and EC-PIC codes are successfully benchmarked and validated, which could be directly used in future large volume plasma device modeling.

In this paper, we described the direct implicit (DI) and explicit energy conserving (EC) particle-in-cell algorithms, which offer a pathway to reduce computational cost of 2D simulations of modern low-temperature large volume plasma devices. We provide a comprehensive explanation of the derivations of the two algorithms and implemented them into two different PIC codes, i.e., EDIPIC-2D and the LTP-PIC code. These algorithms were tested against fully resolved traditional MC-PIC simulations of a CCP discharge operated with radio frequency. Good agreement was found between the algorithms when the time step was kept below that required to resolve the electron plasma period, even when the Debye length is under-resolved. Furthermore, these simulations were accomplished in significantly less time due to the smaller number of total cells and macro-particles. When the time step was increased, agreement was still favorable; however, numerical heating/cooling led to larger deviations in electron temperature. In order to address this issue, we demonstrate a way to estimate this numerical heating, which might be applicable in even in more complex CCP configurations. The outcome of this investigation has shown that the two approaches could be used to perform significantly faster simulations, with only a minimal cost in accuracy.

We noticed that the DI algorithm, despite being very useful, does not accurately conserve energy; therefore, the future work will consider the use of fully implicit PIC algorithms. In this paper, we also directed our focus to the CCP discharge; however, both algorithms are amenable to more complex configurations, since their implementations are able to handle magnetic fields. In future endeavors, we will also explore the application and limitations of these approaches when applied to other low-temperature plasma devices that are pertinent to material processing.

The PPPL group was funded by the U.S. Department of Energy through PPPL Laboratory Directed Research & Development (LDRD) program. The PPPL group used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231. S.S. would like to thank ANTYA HPC facility at the Institute for Plasma Research, Gandhinagar, India.

The authors have no conflicts to disclose.

Haomin Sun: Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Soham Banerjee: Data curation (equal); Formal analysis (equal); Validation (equal). Sarveshwar Sharma: Formal analysis (equal); Investigation (equal); Writing – review & editing (equal). Andrew Tasman Powis: Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Alexander V. Khrabrov: Methodology (lead); Writing – original draft (supporting). Dmytro Sydorenko: Methodology (equal). Jian Chen: Writing – review & editing (equal). Igor D. Kaganovich: Funding acquisition (lead); Supervision (lead).

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