The Particle-in-Cell (PIC) method, a cornerstone in plasma modeling, is widely employed for its ability to simulate kinetic phenomena in device-scale domains. Part of what makes this possible is that computational macroparticles represent many physical particles. It converges under certain constraints, including a grid spacing that resolves the Debye length and a time step small enough to respect the Courant–Friedrichs–Lewy condition and plasma frequency stability limit. Here, we introduce a new constraint necessary to avoid Artificial Correlation Heating (ACH). This requires that the macroparticle coupling strength be smaller than one, $\Gamma w<1$, where $\Gamma w\u2261\Gamma w2/3$, $\Gamma =Z2e2/(4\pi \epsilon oakBT)$ is the physical coupling strength, and *w* is the macroparticle weight. This is particularly relevant to 3D simulations of dense plasmas, which are becoming common with modern computing power. If this condition is violated, the finite macroparticle weight artificially enhances the coupling strength and causes the plasma to heat until the macroparticle coupling strength is near unity, depending on the grid resolution. A comprehensive model of ACH is developed that incorporates electron density, temperature, macroparticle weight, and grid resolution. It is then tested using PIC simulations, delineating the boundaries of the method's applicability and offering a predictive framework for ACH. Moreover, the research explores a runaway heating process induced by ACH in the presence of ionization, which can lead to numerical instability. A conclusion of this study is that the onset of ACH can impose a more stringent constraint on the macroparticle weight and average number of macroparticles per cell than what is typically expected, particularly in 3D simulations of dense plasmas.

## I. INTRODUCTION

The Particle-in-Cell (PIC) method is a cornerstone in plasma modeling, offering a versatile and robust computational tool widely used across various plasma regimes. The method's universality stems from its ability to effectively simulate plasmas under diverse conditions, from low-temperature plasmas^{1–5} to high-energy astrophysical phenomena.^{6,7} Applications of the PIC method are extensive, including space propulsion,^{8–10} fusion research,^{11–14} and laser–plasma interactions.^{15} Given the transformative impact of these applications in numerous fields, understanding and accurately modeling plasma behavior becomes paramount. The PIC method has been instrumental in advancing our understanding of plasma dynamics, thus facilitating the development and optimization of plasma-based technologies and applications.

The PIC method is particularly effective in bridging the gap between microscopic particle dynamics and macroscopic plasma behavior. In PIC simulations, macroparticles are used to represent groups of physical particles, enabling the method to capture large-scale dynamics of plasmas with computational efficiency.^{16–18} As with any simulation, certain constraints must be met to ensure convergence. The Courant–Friedrichs–Lewy condition should be satisfied to accurately resolve the plasma frequency.^{17} In momentum conserving schemes, it is usually expected that the Debye length must be resolved to prevent grid heating [i.e., the finite grid instability (FGI)],^{17,19} though grid heating can be controlled using energy conserving integration schemes, or higher order shape functions and filtering methods.^{19–25} A third condition is that there is a sufficient number of macroparticles per cell to ensure a statistically representative representation of the plasma behavior. By adhering to these principles, the PIC method is expected to provide an accurate solution of the plasma kinetic equation. Here, we describe an additional constraint called artificial correlation heating (ACH) that should also be considered, particularly in 3D simulations of dense plasmas.

^{26}Here, we propose that in order to avoid ACH, a PIC simulation must obey a criterion that the macroparticle coupling parameter be less than one, $\Gamma ss\u2032w<1$, for each combination of interacting species

*s*and $s\u2032$. The macroparticle coupling parameter is defined as

*s*and $s\u2032$, $Tss\u2032$ is a mean temperature, and

*w*is the macroparticle weight. Here, we will concentrate on like-particle interactions ( $s=s\u2032$), so $ass\u2032=as=(3/4\pi ns)1/3$ and $Tss\u2032=Ts$ is the temperature of species

*s*. The $\Gamma ss\u2032w<1$ condition limits how high a macroparticle weight can be for a given density and temperature. If this condition is not met, then the weakly coupled particles will be represented numerically as strongly coupled macroparticles, which induces ACH and can significantly raise the temperature of that species on a short timescale characteristic of the plasma period of that species $\omega ps\u22121=(\epsilon oms/Zs2e2ns)1/2$.

Artificial correlation heating is similar to the physical effect of disorder-induced heating (DIH) that arises in strongly coupled plasmas.^{27,28} DIH occurs when a gas is ionized, causing interparticle interactions to transition from short to long range. This causes the charged particles to separate more, on average, leading to a reduction of the potential energy and a corresponding increase in kinetic energy (heating). As an example, recent work showed that ions in atmospheric pressure plasmas can be strongly coupled if the ionization fraction is large enough.^{29} It was shown that for this reason, PIC simulations are generally not well suited to simulate these conditions, missing effects such as DIH, unless they are run in very restrictive conditions like using a macroparticle weight of unity ( $w=1$) and a grid spacing that is much smaller than the average interparticle spacing.^{26} With these parameters, PIC effectively operates as a molecular dynamics simulation but becomes impractical for simulating macroscopic systems. Like DIH, ACH is also caused by the Coulomb repulsion between particles and the conversion of potential to kinetic energy associated with a reconfiguration to a lower potential energy state. However, in the case of ACH, the strong coupling is induced artificially by a large macroparticle weight. For this reason, it can apply to a species that is physically weakly coupled ( $\Gamma \u226a1$) if the macroparticle weight is sufficiently large that $\Gamma w>1$. For example, electrons at 2 eV and a density of $1\xd71025$ m^{−3}, corresponding to that of an initial atmospheric pressure gas, have $\Gamma e\u22480.1$. Thus, a macroparticle weight $w>10$, which is common, will lead to the $\Gamma w>1$ conditions where ACH is possible.

Acknowledging the significance of ACH, a comprehensive model is developed that incorporates density, temperature, macroparticle weight, and grid resolution. Validated against PIC simulations, this model serves as a predictive tool to delineate a limit of applicability of the PIC method. Considering 3D simulations, it is shown that ACH provides an upper limit on the density scaled by the macroparticle weight squared, $nsw2$, for a fixed temperature. In addition, the consequences of violating the ACH condition ( $\Gamma sw<1$) are explored. In particular, the magnitude of ACH can be reduced if a large number of macroparticles are present per cell. Specifically, it is found that if $\u2273240$ particles per cell are present, then the temperature rise from ACH remains small even if $\Gamma sw>1$. However, this approach is not entirely satisfactory for momentum conserving algorithms because it implies not resolving the Debye length, and therefore, grid heating can occur. This is because $\Gamma sw$, the Debye length $\lambda Ds$, the average number of macroparticles per cell $Nc$, and cell size $\Delta x$ are related as $\Delta x/\lambda Ds=(4\pi Nc/3)1/33\Gamma sw$, see Sec. II B. Thus, if $\Gamma sw>1$ and $Nc>1$, then the Debye length is necessarily unresolved. Hence, to avoid ACH while reducing the computational cost implies unresolving the Debye length for certain conditions. This, in turn, indicates that energy conserving methods are more suitable for such conditions to control grid heating.^{19}

Furthermore, ACH is more relevant to 3D simulations, where it is common to have a relatively small number of macroparticles per cell, below the criteria found in this work.^{30–33} The reason ACH is less relevant to 1D or 2D simulations is that the process of projecting particles to a lower dimensional space effectively decreases the distance between them, mitigating ACH. ACH is most relevant to dense plasmas because they have physical values of the Coulomb coupling strength, $\Gamma $, that are large enough that the product $\Gamma w2/3$ can often exceed one. When the plasma is dense (and cold) enough, this can lead to situations where the macroparticle weight locally takes values less than unity if one simultaneously resolves the Debye length and has many macroparticles per cell. Unfortunately, many PIC papers neglect to give the macroparticle weights in their simulations; however, in one of the authors' prior publications,^{34} the macroparticle weights were $O(10\u22123)$. In order to better resolve massive spatial and temporal differences in physical densities, PIC codes often make use of a variable macroparticle weight and employ automatic merge and splitting algorithms to ensure that each cell contains many particles. This was the case in Ref. 34 where merge algorithms kept the number of macroparticle electrons per element to 100, and the domain had tetrahedral elements with a maximum size of $0.3\u2009\mu $m ( $Vc\u22481.48\xd710\u221223$ m^{3}) to resolve the Debye length and a plasma density of $1021$ m^{−3} resulting in macroparticle weights much less than one. Variable macroparticle weights can easily result in weights less than one without the numerical scientist's knowledge, and we have shown in previous work^{26} that a macroparticle weight less than unity is an invalid representation of physical reality, missing important physical effects such as DIH.

Finally, it is shown that ACH can induce a runaway heating process in simulations that include ionization of a neutral gas. As the electron temperature artificially increases due to ACH, it can trigger nonphysical ionization events, further increasing the electron density and perpetuating the cycle of heating and, ultimately, resulting in numerical instability. This positive feedback loop highlights the critical need for a better understanding and careful application of the PIC method, especially when simulating plasmas with high densities and large macroparticle weights.

The remainder of this paper is organized as follows. In Sec. II A, the numerical implementation of the PIC method utilized in this work is described. This is followed by an example of how ACH can be induced by increasing the macroparticle weight in Sec. II B. Then, a model for predicting ACH accounting for typical input parameters in a PIC simulation is described in Sec. III. Section IV B shows the PIC results and influence of the grid resolution on ACH as well as an extension of the model to include this effect. Section IV C shows the effect of shape function order on ACH. The limit of applicability of the PIC method, accounting for ACH, is described in Sec. IV G. Finally, the positive feedback loop of ACH when an electron impact ionization process is included is shown in Sec. IV D.

## II. PIC SIMULATIONS

### A. Simulation setup

An in-house developed electrostatic 3D3V PIC code was utilized, similar to what was described in previous work.^{26} The code adheres to the conventional method of interpolating particle charges onto a structured uniform grid using shape functions. Subsequently, field equations are solved, and the resulting fields are interpolated back to the particle positions to integrate the equations of motion using a momentum conserving scheme. This process is repeated throughout each iteration of the PIC simulation.^{17,26} Motivated by recent atmospheric pressure discharge simulations with high ionization fractions, which model electrons as particles with the PIC method and ions as a background fluid,^{35} we focus on studying ACH of the electron species. However, the results also pertain to ions if they are treated as particles in a PIC simulation, since the interaction between each species comes from the same potential, the Coulomb potential.

*V*is the volume of the simulation domain. The Poisson equation for the electrostatic potential was then solved using a 3D-FFT.

^{17}The electric field components were obtained by numerically differentiating the electric potential along each direction using finite central differences of order $O(2)$ and periodic boundary conditions. No external electric fields were included.

^{17}Shape functions of order 2 were used, unless otherwise specified. The shape functions used are described in Ref. 36, where the case of order two is

^{26}

The simulation setup is such that in the absence of ACH, density and temperature remain constant in time. So any deviation from the initial density and temperature can be prescribed to a numerical error. The number of macroparticles used in each simulation was given by $New=neV/w$, where $ne$ is the electron density, *V* is the volume of the simulation domain, and *w* is the macroparticle weight. The volume of the simulation domain was set by $V=L3$, where $L=(Nx\u22121)\Delta x$, $Nx$ is the number of nodes in each direction, and $\Delta x$ is the cell size. The simulation domain was cubic with periodic boundary conditions for fields and particles. The time step used for all the simulations was $\Delta t=10\u22122\omega pe\u22121$, where $\omega pe=(e2ne/\epsilon 0me)1/2$ is the electron plasma frequency. The positions of the macroparticles were initialized with a uniform random distribution and the velocities with a Maxwellian distribution at an initial temperature of 3 eV. No ion or neutral species were included since the objective of this part is to study artificial correlation heating of weakly coupled species, such as electrons, in a timescale given by $\omega pe\u22121$. Thus, ions and neutrals remain stationary at the timescale of interest and are modeled as a background uniform fluid. With this numerical setup, any increase in the electron temperature is attributable to a numerical error, particularly ACH or grid heating.

In Sec. IV D, an ionization module is included to study the possibility of a positive feedback loop between ACH and electron-impact ionization. This is done by including a Monte Carlo collision routine.^{17} In this part of the work, a background neutral uniform Xe gas is included at an initial density $nn=ng\u2212ne$, where $ng$ corresponds to the ideal gas density at atmospheric pressure and room temperature. Then, in each iteration and for each electron in the simulation, the probability of ionization is computed as $P=1\u2212exp(\u2212nn\sigma Xeve\Delta t)$, where $\sigma Xe$ is the corresponding electron impact ionization cross section for Xe^{37} and $ve$ is the electron speed. Then, this probability is compared to a random number *R* from a uniform distribution between 0 and 1. If $R<P$, then the ionization occurs, and a new electron is added at a random position in the simulation domain with a velocity sampled from a Maxwellian distribution at the initial electron temperature of 3 eV. After repeating this process for each electron, the background gas density is updated, subtracting from $nn$ the increase produced in $ne$, following the balance $nn=ng\u2212ne$. In the absence of ACH and within a weakly coupled plasma simulation using unity macroparticle weight, a steady state is achieved with this simulation setup. To maintain this steady state, particles must be removed at the same rate at which they are created, at conditions where ACH is absent. This process involves sampling an electron velocity, $ve$, from a Maxwellian distribution at an initial temperature of 3 eV for each particle, followed by employing the MCC scheme to remove the particle. Therefore, if $w=1$ (no ACH is present), the electron density remains constant during the simulation when the ionization module is active. In this case, the rates of ionization and removing of electrons are the same.

### B. Artificial correlation heating

To provide a concrete example, consider a plasma with an electron density of $2.5\xd71019\u2009m\u22123$ and an electron temperature of 3 eV. In this case, the electron Coulomb coupling parameter is approximately $\Gamma ee=2\xd710\u22123$, indicative of a weakly coupled regime. Using a time step of $10\u22122\omega pe\u22121$, grid resolution $\Delta x/\lambda De$ = 0.5, and unity macroparticle weight ( $w=1$), no numerical heating is observed, as expected since $\Gamma ew\u226a1$, see Fig. 1. However, if the macroparticle weight is increased enough that $\Gamma eew>1$, significant heating is observed. Since the heating is present, even though $\Delta x/\lambda De<1$, this is attributed to ACH rather than finite grid heating. For the particular case of $w=105$, the electron temperature increases to a peak of 9 eV within $1.5\u2009\omega pe\u22121$. This rise in temperature is the result of a large macroparticle weight causing a larger than unity macroparticle coupling parameter, $\Gamma eew=4.87$, as deduced from Eq. (1).

The magnitude of this numerical heating is not only solely attributable to the enhanced macroparticle coupling strength but also the limited number of macroparticles per cell $Nc$. This is illustrated in Fig. 2. Considering fixed plasma parameters and domain, a case with a smaller macroparticle weight will have more particles per cell on average (case a). This will lead to a largely uniform density from cell to cell, a small electric field force, and, therefore, little or no ACH. In contrast, a larger macroparticle weight leads to fewer particles in the domain and a larger variation in the initial macroparticle count from cell to cell (case b). This leads to an artificially enhanced electric field from cell to cell, which causes macroparticles to separate until the spatial configuration is such that the Coulomb potential energy is minimized. This process is similar to DIH,^{29} with the difference that the change in potential energy is artificially enhanced due to the large macroparticle weight, and that the grid resolution influences the amount of heating that can be observed. Since the total energy is conserved in a simulation, the decrease in the potential energy leads to a significant increase in temperature, resulting in artificial heating.

w
. | $Nc$ . | $\Gamma eew$ . | $\Delta x/aeew$ . |
---|---|---|---|

1 | $5.34\xd710$ | $2.26\xd710\u22123$ | 6.07 |

$103$ | $5.34\xd710\u22122$ | $2.26\xd710\u22121$ | $6.07\xd710\u22121$ |

$104$ | $5.34\xd710\u22123$ | 1.05 | $2.82\xd710\u22121$ |

$105$ | $5.34\xd710\u22124$ | 4.87 | $1.31\xd710\u22121$ |

w
. | $Nc$ . | $\Gamma eew$ . | $\Delta x/aeew$ . |
---|---|---|---|

1 | $5.34\xd710$ | $2.26\xd710\u22123$ | 6.07 |

$103$ | $5.34\xd710\u22122$ | $2.26\xd710\u22121$ | $6.07\xd710\u22121$ |

$104$ | $5.34\xd710\u22123$ | 1.05 | $2.82\xd710\u22121$ |

$105$ | $5.34\xd710\u22124$ | 4.87 | $1.31\xd710\u22121$ |

Simulations, where a large macroparticle weight is used, such that $\Gamma eew\u22731$, and the grid resolution resolves the Debye length, will inevitably lead to $Nc\u226a1$. As a consequence, the strong correlations will be resolved by the grid, and ACH will be significant.

## III. ACH MODEL

Artificial correlation heating arises from macroparticles moving from an initially randomly distributed configuration to the lowest potential energy configuration, as shown in Fig. 2. Here, a model for ACH is described, based on the macroparticle effect on the potential energy and conservation of energy arguments. First, the potential energy is calculated from statistical mechanics assuming infinitesimal grid resolution. Then, the change in potential energy from an initial uncorrelated state to an ordered state is calculated to predict the change in kinetic energy of macroparticles due to ACH. This model is then extended in Sec. IV B to any grid resolution $\Delta x/aeew$.

^{38}

*r*and thickness

*dr*centered on a chosen particle. Thus, $g(r)$ is the density of electrons surrounding an electron at $r=0$, normalized by the background density. The change in potential energy from the initial to final spatial configuration is then

^{38}

*T*is the temperature, $h(r)=g(r)\u22121$, and $h\u0302(k)$ denotes the Fourier transform of $h(r)$. The hypernetted chain approximation has been shown to be an accurate closure to model the radial distribution function for a one component plasma across different coupling regimes.

^{38–40}

*w*, the equation for the electron temperature increase due to ACH in a PIC simulation is

*w*, and $Te0$. It is important to underscore that the temperature $Temax$ predicted with this model neglects the effect of the grid resolution; thus, it assumes that the average distance between macroparticles is resolved by the grid, thus $Nc\u226a1$. Modifications to account for the grid resolution are provided in Sec. IV B.

## IV. DISCUSSION

### A. Limit of fine grid resolution

Figure 4 illustrates the maximum electron temperature obtained from PIC simulations across different electron densities and macroparticle weights. Each simulation was conducted with a grid spacing of approximately $\Delta x/aeew\u22480.1$, resulting in an average of $Nc=2.35\xd710\u22124$ macroparticles per cell. In this limit, the grid resolves the average distance between macroparticles. The specific parameters employed in these simulations are detailed in Table II. It is clear that as the density and macroparticle weight increases, the influence of ACH on the electron temperature becomes more significant. This agrees with the fact that increasing both $ne$ and *w* increases the macroparticle coupling parameter $\Gamma eew$. The theoretical model demonstrates good agreement with the PIC simulation results, showing its validity.

$ne$ (m^{−3})
. | w
. | $\Gamma eew$ . | $\Delta x/\lambda De$ . |
---|---|---|---|

$1021$ | 1 | $7.74\xd710\u22123$ | $1.52\xd710\u22122$ |

10^{2} | $1.67\xd710\u22121$ | $7.07\xd710\u22122$ | |

10^{4} | 3.60 | $3.28\xd710\u22121$ | |

$1022$ | 1 | $1.67\xd710\u22122$ | $2.24\xd710\u22122$ |

10^{2} | $3.60\xd710\u22121$ | $1.04\xd710\u22121$ | |

10^{4} | 7.74 | $4.82\xd710\u22121$ | |

$1023$ | 1 | $3.60\xd710\u22122$ | $3.28\xd710\u22122$ |

10^{2} | $7.74\xd710\u22121$ | $1.52\xd710\u22121$ | |

10^{4} | $1.67\xd710$ | $7.07\xd710\u22121$ |

$ne$ (m^{−3})
. | w
. | $\Gamma eew$ . | $\Delta x/\lambda De$ . |
---|---|---|---|

$1021$ | 1 | $7.74\xd710\u22123$ | $1.52\xd710\u22122$ |

10^{2} | $1.67\xd710\u22121$ | $7.07\xd710\u22122$ | |

10^{4} | 3.60 | $3.28\xd710\u22121$ | |

$1022$ | 1 | $1.67\xd710\u22122$ | $2.24\xd710\u22122$ |

10^{2} | $3.60\xd710\u22121$ | $1.04\xd710\u22121$ | |

10^{4} | 7.74 | $4.82\xd710\u22121$ | |

$1023$ | 1 | $3.60\xd710\u22122$ | $3.28\xd710\u22122$ |

10^{2} | $7.74\xd710\u22121$ | $1.52\xd710\u22121$ | |

10^{4} | $1.67\xd710$ | $7.07\xd710\u22121$ |

### B. Effect of grid resolution

Previous work indicated that the grid resolution has an impact on the temperature increase due to ACH.^{26} To further investigate this effect, a series of simulations were executed varying the grid resolution alongside different electron densities and macroparticle weights. In this set of simulations, for each combination of parameters, 20 independent simulations with different initial random seeds were conducted, and the maximum temperatures were averaged. Table III shows the parameters used. The outcomes, presented in Fig. 5(a), reveal that an increase in cell size correlates with a reduction in ACH, which is attributed to an increase in the average number of macroparticles per cell $Nc=(3/4\pi )(\Delta x/aeew)3$, as illustrated in Fig. 2. This increase in the cell size leads to a failure to resolve the artificially enhanced correlations, decreasing ACH.

n (m_{e}^{−3}) (2.5×)
. | w
. | $\Gamma ee$ . | $\Gamma eew$ . | $\Delta x/aeew$ . | N (2.37×)
. _{c} | $\Delta x/\lambda De$ . |
---|---|---|---|---|---|---|

$1020$ | $105$ | $4.87\xd710\u22123$ | 10.50 | $10\u22121$ | $10\u22124$ | 0.561 |

1 | $10\u22121$ | 5.61 | ||||

10 | $102$ | 56.1 | ||||

$1023$ | $103$ | $4.87\xd710\u22122$ | 4.87 | $10\u22121$ | $10\u22124$ | 0.382 |

1 | $10\u22121$ | 3.82 | ||||

10 | $102$ | 38.2 | ||||

$1024$ | $102$ | $1.05\xd710\u22121$ | 2.26 | $10\u22121$ | $10\u22124$ | 0.261 |

1 | $10\u22121$ | 2.61 | ||||

10 | $102$ | 26.1 |

n (m_{e}^{−3}) (2.5×)
. | w
. | $\Gamma ee$ . | $\Gamma eew$ . | $\Delta x/aeew$ . | N (2.37×)
. _{c} | $\Delta x/\lambda De$ . |
---|---|---|---|---|---|---|

$1020$ | $105$ | $4.87\xd710\u22123$ | 10.50 | $10\u22121$ | $10\u22124$ | 0.561 |

1 | $10\u22121$ | 5.61 | ||||

10 | $102$ | 56.1 | ||||

$1023$ | $103$ | $4.87\xd710\u22122$ | 4.87 | $10\u22121$ | $10\u22124$ | 0.382 |

1 | $10\u22121$ | 3.82 | ||||

10 | $102$ | 38.2 | ||||

$1024$ | $102$ | $1.05\xd710\u22121$ | 2.26 | $10\u22121$ | $10\u22124$ | 0.261 |

1 | $10\u22121$ | 2.61 | ||||

10 | $102$ | 26.1 |

This phenomenon is better understood by analyzing the radial distribution function obtained from PIC for a specific case of $ne=1022\u2009m\u22123$ and $w=1000$, as depicted in Fig. 6. An increase in cell size leads to an increase in the $g(r)$ at short distances. This is a consequence of the underestimation of the Coulomb potential within a cell, inherent to the PIC method. At a threshold where $\Delta x/aeew\u22737$, $g(r)$ converges to that of an uncorrelated state, and, therefore, ACH is effectively turned off. The numerical grid used in PIC simulations makes the potential around a charged particle asymmetric. Hence, the HNC approximation no longer serves as an accurate tool to determine the $g(r)$, necessitating an alternative approach to model grid resolution effects.

*w*. As shown in Fig. 5(a), the electron temperatures predicted by Eq. (13) closely match those observed in the PIC simulations, making the model given by Eqs. (10) and (13) a useful tool to predict ACH at any conditions and grid resolution.

It is important to underscore that avoiding ACH when $\Gamma eew\u22731$ requires that the grid resolution be $\Delta x/aeew\u227310$. This is equivalent to a number of particles per cell of $Nc\u2273240$, a value larger than what is typically used in 3D PIC simulations, particularly at large plasma densities such as in streamer discharges at atmospheric pressure, vacuum arc discharges, and laser-driven plasmas.^{30–33} According to Eq. (4), this requirement leads to under-resolving the Debye length when $\Gamma eew\u22651$. Thus, avoiding ACH by reducing the grid resolution is not a viable option for momentum conserving algorithms, since not resolving the Debye length will lead to traditional grid heating on timescales much longer than the plasma period.^{17} This indicates that energy conserving methods should be used at such conditions. Therefore, ACH imposes an additional constraint on PIC simulations, necessitating careful consideration when modeling plasmas with large macroparticle weights.

### C. Effect of shape function order

Previous work has shown that the order of shape functions used in a PIC simulation can affect how well correlations are resolved.^{26} As the order of shape function is increased, correlations are smoothed out, effectively decreasing the artificially enhanced electric field. As a result, ACH would be expected to decrease with order of the shape function. To test how this influences ACH, the set of simulations with $ne=2.5\xd71023$ m^{−3} from Table III was repeated but varying the order of shape functions from 2 to 6.^{36} Figure 7 shows the obtained maximum electron temperature for each simulation as a function of the grid resolution. It is observed that as the order of the shape function is increased, ACH does decrease. However, this effect does not make a significant enough change to justify the use of higher order shape functions to control ACH. On the contrary, as stated in the grid resolution effect analysis, increasing the cell size such that a minimum number of macroparticles per cell of $\u2273$240 is reached can avoid observing ACH without an increase in computational cost. However, this last approach leads to a decrease in the grid resolution in terms of the Debye length, which can trigger grid heating on a longer timescale.

### D. Numerical instability produced by ACH and ionization

Artificial correlation heating is potentially unstable in PIC simulations when ionization of a neutral gas is included. This can happen when the bulk of the electron energy distribution function is below the peak ionization cross section energy, and, thus, increasing the temperature will increase the ionization rate. For high temperature plasmas, $O$(100 eV), the ionization cross section decreases with increasing energy so there should be no instability. However, for low-temperature plasmas, $O$(1 eV), the temperature rise due to ACH will significantly increase the ionization rate through the steep electron energy dependence of the electron–neutral impact ionization cross section in the few eV energy range, see, for example, Ref. 37. The increase in the ionization rate following ACH corresponds to an artificial increase in the electron density, which then leads to more ACH, further increasing the electron temperature and potentially inducing a runaway heating process. To study the possibility of this runaway process, a set of simulations were run at an initial electron density of $ne=2.5\xd71019$ m^{−3}, temperature $Te=3$ eV, and initial grid resolution $\Delta x/\lambda De=0.5$, varying the macroparticle weight from 1 to $105$ as indicated in Table IV. Each simulation was repeated 20 times varying the random seed number, and the results displayed are the average of all the simulations.

n (m_{e}^{3}) $(2.5\xd7)$
. | $\Gamma ee$ $(2.3\xd7)$ . | $\Delta x\lambda De$ . | w
. | $\Gamma eew$ . | $\Delta x/aeew$ . | N $(4.9\xd7)$
. _{c} |
---|---|---|---|---|---|---|

$1019$ | $10\u22123$ | 0.50 | 1.0 | $2.3\xd710\u22123$ | 6.1 | $101$ |

10^{3} | $2.2\xd710\u22121$ | 0.61 | $10\u22122$ | |||

10^{4} | 1.1 | 0.28 | $10\u22123$ | |||

10^{5} | 4.9 | 0.13 | $10\u22124$ |

n (m_{e}^{3}) $(2.5\xd7)$
. | $\Gamma ee$ $(2.3\xd7)$ . | $\Delta x\lambda De$ . | w
. | $\Gamma eew$ . | $\Delta x/aeew$ . | N $(4.9\xd7)$
. _{c} |
---|---|---|---|---|---|---|

$1019$ | $10\u22123$ | 0.50 | 1.0 | $2.3\xd710\u22123$ | 6.1 | $101$ |

10^{3} | $2.2\xd710\u22121$ | 0.61 | $10\u22122$ | |||

10^{4} | 1.1 | 0.28 | $10\u22123$ | |||

10^{5} | 4.9 | 0.13 | $10\u22124$ |

Figure 8(a) shows the evolution of the electron temperature for each simulation over approximately 8 ps. It is clear that if the macroparticle weight is large enough that $\Gamma eew>1$, ACH significantly increases the electron temperature, and that this considerably increases the ionization rate, as shown in Fig. 8(c). Here, the ionization rate is calculated from the Xe electron-impact ionization cross section^{37} and a Maxwellian distribution at the corresponding electron temperature. Comparing the curves in Fig. 8(c) with macroparticle weights from $103$ to $105$, it is noticeable how the ionization rate can significantly increase by one or more orders of magnitude. This is reflected as an increase in the electron density with time at a faster than an exponential rate, as shown in Fig. 8(b).

Since the simulation domain volume as well as the cell size $\Delta x$ remains constant, increases in density and temperature affect the grid resolution. Figures 8(d) and 8(e) show how the grid resolution changes in terms of the interparticle distance and the electron Debye length, respectively. It is clear from Fig. 8(d) that fixing the initial grid resolution $\Delta x/\lambda De=0.5$ but increasing the macroparticle weight increases the grid resolution in comparison with the average interparticle distance. As was demonstrated in Sec. IV B, a grid that resolves $aeew$ is necessary for ACH to occur. It is also observed that for the case $w=105$, the significant change in the electron density increases $\Delta x/aeew$ and $\Delta x/\lambda De$ over time, potentially limiting ACH after the density has built up sufficiently. It is noteworthy that this implies that if an adaptive mesh refinement method was used to maintain $\Delta x/\lambda De=0.5$, ACH would be larger than what was observed here; further increasing the temperature and, thus, electron density to values larger than what it is shown in Figs. 8(a) and 8(b). A conclusion is that ACH combined with ionization can lead to a runaway heating mechanism.

### E. Artificial correlation heating in reduced dimensions

Running PIC simulations in reduced spatial dimensions is a common approach for systems with some kind of spatial symmetry. When using reduced spatial dimensions, the positions of all particles within the physical volume are projected onto an element in the corresponding reduced dimension.^{4} This method reduces the distance between numerical particles. Hence, if the physical Debye length is resolved, and a lower dimensionality of a problem is simulated, the average number of macroparticles per cell is significantly increased. This considerably reduces ACH, making it relevant only in 3D PIC simulations.

### F. Significance and relevance

The ACH criteria described earlier are mostly likely to be an important consideration when the plasma is dense or cool, so that the physical coupling parameter $\Gamma $ is not extremely small. In this case, a modest macroparticle weight can cause the macroparticle coupling parameter to exceed unity $\Gamma w>1$. Here, we mention a few examples where this is likely to be an important consideration. We note that it may be difficult to find explicit examples of violating the ACH criterion in the literature for two reasons: (1) ACH causes the plasma to heat to the point where the effect turns off ( $\Gamma w\u22481$), and this happens very fast (on a timescale of $\omega p\u22121$). So one would not expect to observe plasma conditions that violate the criterion. Rather, ACH may cause an unphysical, but unnoticed, heating effect. (2) Literature rarely reports the macroparticle weight, and many codes use a variable macroparticle weight. This makes it impossible to evaluate $\Gamma w$.

An illustrative example of where ACH could play a significant role is streamer discharges at atmospheric pressure that exhibit a high ionization degree at the streamer tip. For example, in Ref. 30, 3D PIC simulations were employed to investigate streamer discharges under atmospheric conditions. The authors imposed an upper limit on the electron density, set at $2.5\xd71021$ m^{−3}, which corresponds to an ionization degree of $10\u22124$ at atmospheric pressure. A target of 32 macroparticles per cell was maintained, with the smallest cell size being 1.25 $\mu $m. Assuming an electron temperature of 3 eV, the parameters stipulated result in a macroparticle weight of $w=132$ at the streamer's tip. These specifications yield a macroparticle Coulomb coupling parameter, $\Gamma eew$, of approximately 0.3, aligning closely with the limiting criteria established in the present study. This work acknowledged the necessity of capping the electron density to mitigate computational demands.^{30} However, our model suggests that, had a higher cap been set, ACH would likely have been a critical factor under the conditions explored. This underscores the potential of ACH as a notable source of numerical heating in scenarios with elevated ionization fractions, emphasizing its relevance to the research community.

Another example is vacuum arc discharges,^{31} with electron densities $\u22481026$ $m\u22123$. At these conditions, not only disorder-induced heating, studied in our previous work,^{26} is important, but also ACH can be triggered if a macroparticle weight larger than unity is used. Furthermore, 3D PIC simulations of laser-driven plasmas often use between 1 and 120 macroparticles per cell,^{32,33} values that are below the criteria to avoid ACH that we found in our work.

### G. ACH stability criterion

From the results obtained in this work, it is clear that ACH limits how high the electron density and macroparticle weight can be at a given electron temperature and grid resolution, and that it can be avoided if the $\Gamma w<1$ condition is met. However, it also shows that if the $\Gamma w<1$ condition is violated, the magnitude of the resulting ACH depends on the grid resolution and number of particles per cell. It is possible that a small amount of ACH is tolerable. To understand this limit of applicability of PIC simulations, the model given by Eqs. (10) and (13) was utilized to find the limiting curve defined by a 5% of electron heating due to ACH for a range of initial electron temperatures. This is, at each possible $Te$ and at a fixed $\Delta x/\lambda De$, the model is solved for increasing values of $new2$ until a temperature increase of 5% is predicted. This defines a limiting curve for each grid resolution $\Delta x/\lambda De$ at which only lower values of $new2$ can be simulated in order to avoid ACH. Figure 9 shows the limiting curves obtained for grid resolutions $\Delta x/\lambda De$ from 0.1 to 5. As the grid spacing increases, ACH becomes less important, and the limiting curve shifts to higher $new2$, increasing the allowed parameter space. However, ideally, the electron Debye length should be resolved in order to avoid grid heating, and this implies that there is a range of $new2$ values at each temperature that cannot be simulated using the PIC method.

Sometimes, it is not imperative to resolve the Debye length, and the consequent grid heating can be controlled by using energy conserving integration schemes as well as higher order shape functions and filtering methods.^{19–25} In such cases, instead of avoiding ACH by reducing the macroparticle weight until $\Gamma eew\u22641$, it can be controlled by increasing the cell size as described in Sec. IV B. This can also be observed from Fig. 9, where as $\Delta x/\lambda De$ is increased, the ACH limiting curve shifts to higher $new2$ values.

Figure 10 provides a simple procedure for ensuring that ACH remains below a tolerable level in PIC simulations. Once the input parameters $(ne,Te,w)$ are chosen, the macroparticle coupling parameter should be calculated using Eq. (1). If $\Gamma eew\u22641$, then the standard PIC criteria can be used. However, if $\Gamma eew\u22651$, the remaining question to ask is if it is necessary to resolve the electron Debye length.^{19–25} If the cell size $\Delta x$ can be larger than the Debye length, then to avoid or control ACH $\Delta x$ should be increased until $Nc\u2248240$. However, if the Debye length must be resolved, the temperature increase due to ACH can be calculated using Eqs. (10) and (13). If the expected temperature increase is significant for a particular application, then one must reduce the macroparticle weight *w* until ACH is reduced to a tolerable level.

## V. CONCLUSIONS

This study describes a novel numerical heating mechanism observed in 3D PIC simulations, artificial correlation heating. A comprehensive model for ACH is developed that incorporates variables such as electron density, temperature, macroparticle weight, and grid resolution. This model, validated against PIC simulations, maps out a new limit of the PIC method's applicability and provides a predictive tool for anticipating the influence of ACH on the electron temperature. The results can also be applied to ACH of ions in an analogous way. Its accuracy in predicting the onset and magnitude of ACH provides an instrument for a more reliable application of the PIC method in various plasma applications.

One of the key insights is understanding how ACH occurs. ACH is similar to the physical effect of DIH that arises in strongly coupled plasmas,^{29} but here is an artificial numerical effect that is introduced by a large macroparticle weight. Like DIH, it is caused by the Coulomb repulsion between particles, which causes them to move to the lowest potential configuration. During this spatial reordering, there is a conversion of potential to kinetic energy that arises at strong coupling, with the difference that here the strong coupling is induced artificially by a large macroparticle weight. Plasma discharges with large ionization degrees such as streamers at atmospheric pressure and vacuum arcs are particularly sensitive to ACH. It is important to note that the results presented in this work are independent of the unit system chosen for a PIC simulation since the criteria to avoid ACH depend on the dimensionless macroparticle coupling parameter and average number of macroparticles per cell. However, the results are expected to be primarily applicable to 3D PIC simulations. This is because at reduced dimensions, the number of macroparticles per cell, at a same weight and grid spacing, is significantly increased when compared to 3D simulations making ACH less important.

The findings highlight that careful selection of parameters like grid resolution and macroparticle weight is crucial to avoid ACH. Moreover, this study brings to light the identification of a runaway heating process induced by ACH, particularly in simulations that include ionization. This produces a numerical instability capable of increasing the plasma density by several orders of magnitude on a short timescale, exposing the importance of identifying and avoiding ACH. The main takeaway of this study is that it describes an upper limit of applicability of the PIC method in terms of the macroparticle weight and grid resolution. A general procedure is provided to guide on how to properly choose PIC input parameters in order to avoid or control ACH.

## ACKNOWLEDGMENTS

The authors thank Dr. Matthew Hopkins and Dr. Lucas Beving for helpful discussions. This work was supported in part by the U.S. Department of Energy under Award No. DE-SC0022201, and in part by Sandia National Laboratories. Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-NA0003525. This article describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**M. D. Acciarri:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **C. Moore:** Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal). **S. D. Baalrud:** Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

## REFERENCES

*Plasma Physics via Computer Simulation*

*Computer Simulation Using Particles*

*Theory of Simple Liquids*