# Benchmark for two-dimensional large scale coherent structures in partially magnetized ExB plasmas
## -Community collaboration & lessons learned

Andrew T Powis, Eduardo Ahedo, Alejandro Álvarez Laguna, Nicolas

Barléon, Enrique Bello-Benítez, Lucas Beving, Jean-Pierre Boeuf, Guillaume

Bogopolsky, Anne Bourdon, Filippo Cichocki, et al.

► To cite this version:

Andrew T Powis, Eduardo Ahedo, Alejandro Álvarez Laguna, Nicolas Barléon, Enrique Bello-Benítez, et al.. Benchmark for two-dimensional large scale coherent structures in partially magnetized ExB plasmas -Community collaboration & lessons learned. 2026. hal-05447477

HAL Id: hal-05447477

https://hal.science/hal-05447477v1

Preprint submitted on 7 Jan 2026

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY-SA 4.0 - Attribution - ShareAlike - International
License
Benchmark for two-dimensional large scale coherent structures
in partially magnetized E×B plasmas - Community
collaboration & lessons learned

Andrew T. Powis¹, Eduardo Ahedo², Alejandro Álvarez Laguna³, Nicolas Barléon⁶,
Enrique Bello-Benítez², Lucas Beving⁴, Jean-Pierre Boeuf⁵, Guillaume Bogopolsky⁶,
Anne Bourdon³, Filippo Cichocki⁷, Bénédicte Cuenot⁶, Andrew Denig⁸, Zoltán
Donkó⁹, Paul-Quentin Elias¹⁰, Miguel P. Encinar², Denis Eremin¹¹, Pablo Fajardo²,
Farbod Faraji¹², Gwenael Fubiani⁵, Laurent Garrigues⁵, Kentaro Hara⁸, Peter
Hartmann⁹, Matthew Hopkins⁴, Igor D. Kaganovich¹, Aaron Knoll¹², Giovanni
Lapenta¹⁴, Thierry E. Magin¹⁵, Alberto Marín-Cebrián², Mario Merino², Pierpaolo
Minelli¹⁶, Mina Papahn Zadeh¹⁷, Pietro Parodi¹³, Federico Petronio³, Maryam Reza¹²,
Andrei I. Smolyakov¹⁷, Dmytro Sydorenko¹⁸, Francesco Taccogna¹⁶, Miles M.
Turner¹⁹, Olivier Vermorel⁶, Willca Villafana¹, and Liang Xu¹¹

¹Princeton Plasma Physics Laboratory, Princeton, New Jersey 08540, USA

²Universidad Carlos III de Madrid, 28911 Leganés, Spain

³Laboratoire de Physique des Plasmas (LPP), CNRS, Sorbonne Université, École
Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

⁴Sandia National Laboratories, Albuquerque, New Mexico 87132, USA

⁵Université de Toulouse, Toulouse INP, CNRS, LAPLACE, Toulouse, France

⁶CERFACS - 42, avenue Gaspard Coriolis, F-31057 Toulouse, France

⁷Nuclear (NUC) Department, ENEA C.R. Frascati, Frascati, Rome, Italy

⁸Stanford University, 496 Lomita Mall, Stanford, California 94305, USA

⁹HUN-REN Wigner Research Centre for Physics, Konkoly Thege Miklós str. 29-33,
1121 Budapest, Hungary

¹⁰DPHY, ONERA, Université Paris-Saclay, Palaiseau, France

¹¹Ruhr University Bochum, Universitaetsstrasse 150, D-44801 Bochum, Germany

¹²Imperial College London, Exhibition Road, London SW7 2AZ, United Kingdom

¹³Aeronautics and Aerospace Department, von Karman Institute for Fluid Dynamics,
Waterloosesteenweg 72, 1640 Sint-Genesius-Rode, Belgium

¹⁴Department of Mathematics, KU Leuven, Celestijnenlaan 200B, 3001 Leuven,
Belgium

¹⁵Aero-Thermo-Mechanics Department, Université Libre de Bruxelles, Avenue F.D.
Roosevelt 50, 1050 Bruxelles, Belgium

¹⁶Institute for Plasma Science and Technology (ISTP), CNR, Bari, Italy

¹⁷University of Saskatchewan, Saskatoon, Saskatchewan S7N 5E2, Canada

¹⁸University of Alberta, Edmonton, Alberta T6G2E1, Canada

¹⁹Dublin City University, Dublin, Ireland

October 27, 2025
Abstract

Low-temperature plasmas are essential to both fundamental scientific research and critical industrial applications. As in many areas of science, numerical simulations have become a vital tool for uncovering new physical phenomena and guiding technological development. Code benchmarking remains crucial for verifying implementations and evaluating performance. This work continues the *Landmark* benchmark initiative, a series specifically designed to support the verification of low-temperature plasma codes. In this study, seventeen simulation codes from a collaborative community of nineteen international institutions modeled a partially magnetized $E \times B$ Penning discharge. The emergence of large scale coherent structures, or rotating plasma spokes, endows this configuration with an enormous range of time scales, making it particularly challenging to simulate. The codes showed excellent agreement on the rotation frequency of the spoke as well as key plasma properties, including time-averaged ion density, plasma potential, and electron temperature profiles. Achieving this level of agreement came with challenges, and we share lessons learned on how to conduct future benchmarking campaigns. Comparing code implementations, computational hardware, and simulation runtimes also revealed interesting trends, which are summarized with the aim of guiding future plasma simulation software development.

# 1 Introduction

The study of low-temperature plasmas (LTPs) is motivated by scientific exploration as well as their importance to numerous critical industrial technologies [1]. LTPs exhibit complex, multiscale, and often chaotic behavior, similar to those observed in other regimes of plasma physics, with additional complexity introduced by chemical reactions and surface processes. On the applications side, LTPs are present in all of the most critical steps in microcircuit fabrication [2], with additional use cases including other forms of materials processing [3], power transmission technologies [4,5], and spacecraft electric propulsion [6,7]. Research into LTPs therefore provides enormous benefits to both scientific and engineering communities.

As in many domains of modern science, simulations have come to form a third foundational pillar of discovery for LTPs. Plasma experiments, even in the low-temperature regime, are challenging to run and diagnose non-intrusively. Furthermore, the high-dimensional, non-linear, and multiscale nature of the governing equations often prevents analytical solution. Simulations therefore play a critical role in solving these equations under realistic conditions, and assist in interpreting experimental results. Despite their high computational cost, LTP simulations are also playing an increasingly important role in computer aided engineering of related technologies [8]. This progress has been aided by the continuous improvement of computing performance, particularly with the widespread adoption of computational accelerators [9], as well as the integration of traditional modeling techniques with machine learning [10]. Nonetheless, computer science and algorithm design for LTPs is still a cutting edge area of research (for further discussion, see Section VII in Ref. [11]).

Many of the interesting LTP phenomena and applications referenced above operate in the kinetic regime, with the plasma governed by the Boltzmann equation in combination with electrostatic or electromagnetic field equations. Due to its intuitive nature and relatively simple implementation, the particle-in-cell (PIC) algorithm has become the most common approach for modeling kinetic plasmas, including LTPs [12–14]. The algorithm uses a mixed Eulerian/Lagrangian framework, with the fields solved on an up to three-dimensional Eulerian grid and particles evolved in a Lagrangian sense, in the up to six-dimensional phase space. Furthermore, collisions between charged species and neutrals are often implemented via Monte Carlo algorithms [15]. PIC is also highly parallelizable, with particles evolved in a way that is nearly embarrassingly parallel, and field solvers can take advantage of mature distributed linear algebra packages [16–18]. Furthermore, the particle and collision routines are amenable to GPU acceleration, making the algorithm even more appealing for use on modern heterogeneous supercomputing architectures [9, 11, 19, 20].

Despite the relatively simple and interpretable nature of PIC, the requirement for versatile, scalable and portable software can significantly complicate code design and implementation. Therefore, as in any other field, code validation and verification play an important role in developing PIC codes for LTPs. Verification can be pursued in several ways, including through unit tests, the method of manufactured solutions [21] and through benchmarking against other codes. In the past several years the LTP simulation community has introduced a series of benchmarks against which developers can test their software, these include the popular Turner benchmark of 2013 [22] as well as the *Landmark*
benchmarks [23–25]. The aim of this work is to extend the *Landmark* series with a new benchmark configuration which models low-frequency and large-scale structures emerging in partially magnetized E×B plasmas. This is pursued through two-dimensional PIC simulations of a Penning discharge cross-section, wherein a low-frequency rotating spoke is commonly observed.

Two previous benchmarks from the *Landmark* series, specifically 2a [24] and 2b [25] modeled collisionless partially magnetized E×B LTPs in the simplified two-dimensional geometry of a Hall thruster propulsion system. Within this regime, the magnetic field strength is sufficient to magnetically confine electrons, whereas the heavier ions are weakly-magnetized and can move freely under the influence of electric fields. Despite being confined by the applied magnetic field, it has been observed numerically, and experimentally, that electrons are eventually transported across the field lines from the cathode to the anode, even when ignoring collisions with other species. The underlying mechanism for this “anomalous” transport is still an active area of investigation, with numerous mechanisms being proposed. These include near-wall conductivity [26, 27] and a wide range of plasma waves and instabilities [28–39], with the true mechanism likely being a combination of effects and also dependent on the device under consideration. One important plasma mode which does not emerge within the Landmark 2 configurations is the so-called rotating plasma “spoke” [40], a low-frequency, large-scale coherent structure observed in numerous E×B experiments [41–92] and simulations [93–123]. Numerical modeling of this mode within the Penning discharge represents an order of magnitude increase in the range of time scales when compared to the physics modeled in the *Landmark 2* cases, providing a significant stress-test of code capabilities.

<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of a typical Penning discharge experiment with a beamline source of energetic electrons emitted from the cathode. The outer cylindrical pressure vessel is grounded, with dielectric end-plates. The approximately uniform axial magnetic field (red arrows) is sustained by electromagnets. Dynamics between the magnetized electrons and weakly magnetized ions, which are produced at the center and diffuse outwards, produces a radial ambipolar electric field (blue arrows) in the opposite direction to that of a typical ambipolar field found in un-magnetized plasmas. This results in azimuthal E×B drift (green arrow), with the streaming of electrons against slowly moving ions leading to formation of the rotating spoke.</figcaption></figure>

Over the past decade there has been a significant effort to develop an understanding of the spoke formation mechanism. It is generally agreed that the spoke is seeded by a gradient-drift instability (GDI), often referred to as a modified or collisionless Simon-Hoh instability (CSHI) [63, 64, 83, 93, 109–111, 114, 124]. Charge separation, resulting from the relative E×B drift of magnetized electrons to the non- or weakly-magnetized ions (or possibly due to diamagnetic drift [125]) fuel linear growth and rotation in the E×B direction. At long times, the spoke can coalesce into a mode or modes with amplitude comparable to the peak density within the device. Various theories have been proposed to
explain the non-linear and saturated behavior of the spoke, including formation from an inverse cascade of the GDI (or possibly higher frequency) waves [107, 109, 119]. It has also been proposed that when combined with a gradient in magnetic field the spoke can transform into an ionization wave, and even reverse direction [40, 108, 110]. One theory also proposes that observations of the spoke are projections of a 3D helical structure [120]. Just as for anomalous transport, the true answer is likely dependent on the specific system under investigation and may result from several of the above effects, as well as coupling to other waves and instabilities within the device. The structure is often strongly correlated with anomalous transport [40, 43, 45, 93, 94, 107, 108, 111], whether due to transport being enhanced by the wave itself, or the presence of a majority of the plasma current within the spoke arms. While there has been much light shed on the spoke formation, structure, and effects on device performance, this brief discussion highlights that the matter is by no means settled and that this complex phenomena is worthy of continued study.

The Penning discharge is an E×B device which exhibits similar behavior to that observed in Hall thrusters, albeit with a simplified geometry [63–65, 70, 72]. This includes the presence of anomalous transport and rotating spokes. The device geometry, shown in Fig. 1, consists of a central plasma source and uniform axial magnetic field. Plasma is confined axially by the dielectric end plates, and radially by the magnetic field. A radial ambipolar field emerges which drives E×B processes and associated instabilities, including the spoke. The axial dynamics, along the magnetic field lines take place at a relatively short timescale when compared to the radial and azimuthal motion of the plasma. It is therefore common to assume that this axis can be decoupled from the radial-azimuthal plane. The numerical benchmark presented in this work will model a two-dimensional slice of the Penning discharge, including E×B motion and formation of the rotating spoke.

As a numerical benchmark, the collisionless Penning discharge is compelling for several reasons. First, despite the complex physics taking place within the discharge, the configuration is relatively simple, and it is therefore less challenging to define a clear and concise benchmark description. Second, the discharge exhibits complex multiscale collective behavior, which should emerge from PIC code simulations. This includes microscale instabilities, such as the electron-cyclotron drift instability and magnetized two-stream instability observed in the Landmark 2a and 2b benchmarks [24, 25], as well as large scale phenomena such as the spoke. Third, although not the original purpose of this effort, the long time-scale and non-uniform nature of the discharge strains code performance with respect to runtime and parallelisation strategy. Fourth, there exist numerous natural extensions to this benchmark, most importantly the inclusion of collisional effects. Fifth, there are several experimental Penning discharges around the world [67, 72, 73], making this, or a future collisional version of this benchmark a potential pathway towards a code validation exercise, more of which are needed within the LTP community (see for example Refs. [126–129]). Finally, there are many fascinating and unresolved physics questions associated with partially magnetized E×B discharges and the hope is that engaging numerous groups in the modeling of this device may lead to further and deeper investigation into these unsolved problems.

The paper proceeds as follows, Section 2 outlines the setup for the benchmark and numerical diagnostics as well as details of the participating codes. Section 3 compares the results of the different codes as well as highlighting lessons learned during the benchmark and providing suggestions on high-performance computing (HPC) best practices for future code development. Finally, Section 4 offers some concluding thoughts on the effort.

# 2 Methodology

To ensure a successful benchmarking effort, a complete, clear and concise description of the simulation configuration is crucial. In this section, we outline the simulation setup (Section 2.1) and the diagnostics used for code comparisons (Sections 2.2 and 2.3). Section 2.4 provides details on each of the participating institutes and their associated PIC codes.

## 2.1 Penning benchmark configuration

A collisionless Penning discharge configuration with helium-4 ions is considered for this benchmark. Since the spoke rotation frequency scales with the inverse square root of ion mass [93], typical argon and xenon gases would have led to excessive runtimes, dictating the choice of helium-4. There exists
a wide range of algorithms for modeling curved boundaries in PIC simulations, each of which can lead to slightly different simulation results. To control for the effect of such implementations on the benchmarking exercise, we relax the curved geometry of the Penning discharge experiment (see Fig. 1) to that of a grid-aligned square, an approach which has been pursued in numerous other numerical experiments [112,119–121,130]. The two-dimensional square domain {$x, y$} is centered about the origin, with edge lengths $L_x = L_y = 5$ cm (see Fig. 2). The boundary is grounded with electric potential $\phi = 0$ V, and perfectly conducting, such that all particles reaching the boundary are absorbed (deleted from the simulation). There is a uniform, constant, externally applied magnetic field in the $z$ direction with strength $B = 100$ G. With ion temperature $\approx 0.5$ eV, their cyclotron radius is comparable to the device scale, which while still making them weakly magnetized, can subtly influence the physics.

<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic of the Penning discharge benchmark simulation domain, including geometry for the conducting boundaries (thick black lines), injection region (blue circle), and location of the density probe (green dot). The light gray circular boundary is not part of the simulation, but rather indicates the equivalent geometry for a Penning discharge experiment, similar to that shown in Fig. 1.</figcaption></figure>

The simulation is run for a total of $N_t$ time steps each of length $\Delta t$ up to a final time $t_{sim} = 500$ µs. This is sufficient to resolve approximately a dozen rotations of the spoke, and $\Delta t$ is small enough to resolve the electron plasma frequency $\omega_{p,e}$ in the most dense part of the resulting quasi-steady state plasma. The grid is Cartesian and uniform with an equal number of cells in each direction ($N_x = N_y$) and with square cells ($\Delta x = \Delta y$). The cell size is sufficient to appropriately resolve the smallest electron plasma Debye length $\lambda_{D,e}$ at quasi-steady state. The values for all numerical and physical parameters are listed in Table 2.1.

At $t = 0$ s the domain is completely empty, with no macro-particles present. A constant source of electrons and ions is created in the center of the domain, sampled randomly from a Maxwellian velocity distribution function defined by,

$$f_{\alpha}(v_x, v_y, v_z) = f_{\alpha,0} \exp \left( \frac{v_x^2 + v_y^2 + v_z^2}{2T_{\alpha,0}/m_{\alpha}} \right), \quad (1)$$

where $\alpha \in \{e, i\}$ is the species index, $m_{\alpha}$ is the mass, $T_{\alpha,0}$ is the characteristic injection temperature in Joules, and $f_{0,\alpha}$ is a normalizing constant.

The injection profile $\dot{n}_{\alpha}$ is a cylindrical Heaviside function with form,

$$\dot{n}_{\alpha}(x, y) = \begin{cases} \dot{n}_{\alpha,0} & \text{for } x^2 + y^2 \le R_{\text{inj}}^2, \\ 0 & \text{otherwise,} \end{cases} \quad (2)$$
giving a total “injection” current for each species of $I_α = πq_α\dot{n}_{α,0}R_{inj}^2$. Where $q_α$ is the charge and $\dot{n}_{α,0}$ is the constant particle number injection rate for each species.

In these simulations, macro-particles are treated as constant weight, therefore when computing the number of macro-particles to inject per time-step, there is often a remaining fraction of a particle. This is handled by drawing a random number and testing it against the fractional remainder to determine if an additional particle should be injected per time step.

<table><caption>Table 1: Physical and numerical parameters for the Penning discharge benchmark, including details on averaging times for diagnostics.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>helium-4 ion mass</td><td>$m_i$</td><td>7291.712</td><td>$m_e$</td></tr><tr><td>Time step size</td><td>$\Delta t$</td><td>40</td><td>ps</td></tr><tr><td>Number of time steps</td><td>$N_t$</td><td>12,500,000</td><td>-</td></tr><tr><td>Total simulation time</td><td>$t_{sim}$</td><td>500</td><td>µs</td></tr><tr><td>Steps between diagnostic print</td><td>$N_{print}$</td><td>1,250</td><td>-</td></tr><tr><td>Time between diagnostic print</td><td>$\Delta t_{print}$</td><td>0.05</td><td>µs</td></tr><tr><td>Domain length in x/y-direction</td><td>$L_x, L_y$</td><td>5.0 × 10<sup>-2</sup></td><td>m</td></tr><tr><td>Cells in x/y-direction</td><td>$N_x, N_y$</td><td>256</td><td>-</td></tr><tr><td>Cell length in x/y-direction</td><td>$\Delta x, \Delta y$</td><td>1.953 × 10<sup>-4</sup></td><td>m</td></tr><tr><td>Applied magnetic field strength</td><td>$B$</td><td>100</td><td>G</td></tr><tr><td>Injection radius</td><td>$R_{inj}$</td><td>5.0 × 10<sup>-3</sup></td><td>m</td></tr><tr><td>Electron injection current</td><td>$I_{e,0}$</td><td>-20.0 × 10<sup>-3</sup></td><td>A/m</td></tr><tr><td>Electron injection temperature</td><td>$T_{e,0}$</td><td>15.0</td><td>eV</td></tr><tr><td>Ion injection current</td><td>$I_{i,0}$</td><td>8.0 × 10<sup>-3</sup></td><td>A/m</td></tr><tr><td>Ion injection temperature</td><td>$T_{i,0}$</td><td>0.025</td><td>eV</td></tr><tr><td>Macro-particle weight</td><td>$W$</td><td>1.0 × 10<sup>5</sup></td><td>-</td></tr><tr><td>Electron macro-particles injected per time step</td><td>$N_{\Delta t, e}$</td><td>49.932</td><td>-</td></tr><tr><td>Ion macro-particles injected per time step</td><td>$N_{\Delta t, i}$</td><td>19.973</td><td>-</td></tr></tbody></table>

## 2.2 Time averaged diagnostics

Instantaneous snapshots of parameters are printed at time-step intervals of $N_{print}$ which corresponds to $\Delta t_{print}$ time between outputs. The following grid-based quantities are saved at each output time:

- Ion density $n_i(x, y)$

- Electron temperature $T_e(x, y)$

- Plasma potential $\phi(x, y)$

Time-averaged slices of the above listed quantities are used as the primary comparison between codes. Data slices are taken along the x-direction with $x \in [-L_x/2, L_x/2]$ and $y = 0$ (i.e. along the x-axis). This single slice is then averaged over the left and right hand sides of the domain. Specifically, the data at $x \in [-L_x/2, 0]$ is flipped and averaged with the data at $x \in [0, L_x/2]$, providing a single half-profile from the center of the domain to the right hand edge.

Although the system becomes quasi-neutral shortly after initialization, it takes some time for the large scale spoke structure to start-up and for the system to reach quasi-steady state. Considering this, time averaging is performed over $200 \ \mu s \le t < 500 \ \mu s$, with each snapshot printed at interval $\Delta t_{print}$. In total, each quantity is averaged over 6,000 instantaneous snapshots. As will be discussed in Section 3.1, several codes were run to much longer times, averaging quantities over a larger series of data.

Since the plasma potential is a grid-based quantity, measurement at each time step is straightforward. Measuring the plasma species density and temperature requires the interpolation of particle velocity moments onto the grid prior to sampling. This is performed via the standard, bi-linear, cloud-in-cell approach. Dropping the species index, the $m^{th}$ velocity moment ($nU_\beta^m$) in direction $\beta \in \{x, y, z\}$ is computed as,
$$ (nU_{\beta}^{m}) = W \sum_{p=0}^{P} V_{p,\beta}^{m} S\left(\frac{X_{p}-x}{\Delta x}\right) S\left(\frac{Y_{p}-y}{\Delta y}\right), \qquad (3) $$

where P is the total number of macro-particles, $X_p$ and $Y_p$ are the position coordinates of particle p, $V_{p,\beta}$ is the velocity of particle p in direction $\beta$, and $S(x)$ is the standard linear shape function:

$$ S(x) = \begin{cases} 1+x & -1 < x < 0 \\ 1-x & 0 \le x < 1 \\ 0 & \text{otherwise.} \end{cases} \qquad (4) $$

The ion density is therefore computed from the 0th velocity moment, and temperature is defined as the standard deviation of the velocity distribution function about the mean velocity. Specifically, for a given species, cell location, and instantaneous snapshot, (not indexed here), the temperature in direction $\beta$ is defined as,

$$ T_{\beta} = m \left[ \frac{(nU_{\beta}^{2})}{n} - \left( \frac{(nU_{\beta})}{n} \right)^{2} \right], \qquad (5) $$

For reasons that will be discussed in Section 3.2, complications arise when computing the time-average of plasma species temperature for the Penning discharge configuration. Therefore, we consider two definitions in this work.

**Temperature Type I:** The standard definition, computing the temperature using Eq. (5) at each snapshot time and then averaging temporally across all times,

$$ \langle T_{\beta} \rangle_{t}^{I} = m \left\langle \frac{(nU_{\beta}^{2})}{n} - \left( \frac{(nU_{\beta})}{n} \right)^{2} \right\rangle_{t}, \qquad (6) $$

where $\langle \cdot \rangle_t$ denotes time averaging.

**Temperature Type II:** A non-standard definition, choosing to average each velocity moment in time prior to calculating the temperature [131],

$$ \langle T_{\beta} \rangle_{t}^{II} = m \left[ \frac{\langle nU_{\beta}^{2} \rangle_{t}}{\langle n \rangle_{t}} - \left( \frac{\langle nU_{\beta} \rangle_{t}}{\langle n \rangle_{t}} \right)^{2} \right]. \qquad (7) $$

The final temperature is then computed the same for both procedures as, $\langle T \rangle_t^\Theta = \frac{1}{3} \sum_{\beta} \langle T_{\beta} \rangle_t^\Theta$ with $\Theta \in \{I, II\}$.

## 2.3 Measuring spoke rotation frequency

Many techniques could be used to measure spoke rotation frequency, however the chosen approach is both simple and reproducible. We rely on a single ion density probe (see Fig. 2) located at $(x, y) = (\frac{1}{2}L_x, 0)$, on the x axis, halfway between the origin and the right hand boundary. $n_i(\frac{1}{2}L_x, 0)$ is measured at the $N_{print}$ sampling rate for 200 µs $\le t < 500$ µs, and peaks are identified via a peak-finding algorithm [132]. To avoid detection of false peaks due to the noisy data, a minimum peak height of $n_i = 10^{15}$ 1/m³ was set, as well as a minimum separation of 200 data points between peaks (10 µs). Figure 3 shows a plot of density at the probe location against time and identification of the relevant density peaks. The spoke period is taken as the mean of the distance between subsequent peaks and then inverted for the frequency.
<figure><img src="image_4.png" /><figcaption>Figure 3: Plot of ion density at probe location ($x = \frac{1}{2}L_x, y = 0$) against time. Peaks (red dots) are identified via a peak-finding algorithm from the `scipy.signal` package [132].</figcaption></figure>

## 2.4 Code details

A record forty one scientists from nineteen institutions testing seventeen codes were involved in this iteration of the *Landmark* benchmarking effort. These included the Princeton Plasma Physics Laboratory - PPPL (LTP-PIC code [19]), Laplace laboratory at the University of Toulouse (EP-PIC2D [133] and LePIC2D codes [134]), Ruhr University Bochum - RUB (RUBPIC code [25]), Wigner Research Center for Physics (MZ/X-PIC code), Dublin City University - DCU, ONERA (RHEI code [135]), Universidad Carlos III de Madrid - UC3M (PICASO [136, 137] and PICASO-G codes), University of Saskatchewan - USask (the open-source EDIPIC-2D code [138] - also developed with PPPL and the University of Alberta), the Laboratoire de Physique des Plasmas at École Polytechnique - LPP (LPPic code [139, 140]), Stanford University (XPIC2D code [141, 142]), the von Karman Institute for Fluid Dynamics - VKI (the open-source PANTERA code [143]), the European Centre for Research and Advanced Training in Scientific Computing - CERFACS (AVIP-PIC code [144, 145]), the Bari section of the Institute for Plasma Science and Technology of the Italian National Research Council (CNR-ISTP-Bari) with the Nuclear (NUC) Department of ENEA Frascati Research Center (PICCOLO code [146-148]), Sandia National Laboratories (Aleph code [149]) and the Imperial College London - ICL (IPPL-2D code [150]).

Salient features of these codes, the computing architecture used to model the Penning benchmark, and simulation runtimes are shown in Table 2. While most authors implemented the standard explicit momentum-conserving PIC algorithm on a structured grid, some authors considered an unstructured mesh (AVIP-PIC and Aleph codes). Codes were written in either C/C++ or Fortran, with the exception being the IPPL-2D code written with Julia. Many codes implemented external linear algebra packages for solving the Poisson equation, these packages included Hypre [17], PARDISO [151], NAG [152], PETSc [16, 153, 154] and Trilinos [18]. Some codes were parallelised on CPUs with MPI and/or OpenMP and GPUs with OpenACC or CUDA. Authors also considered different ways to decompose the problem into parallel components including particle decomposition, domain decomposition or a hybrid (mixed) approach. A deeper discussion on these approaches and their implications for performance is provided in Section 3.4.
<table><caption>Table 2: Names and details of codes participating in the Penning discharge benchmark.</caption><thead><tr><th>Institute</th><th>Code name</th><th>Language</th><th>Poisson solver</th><th>Parallelization</th><th>Decomposition</th><th>Architecture</th><th>Wall-clock-time</th></tr></thead><tbody><tr><td>PPPL</td><td>LTP-PIC</td><td>C/C++</td><td>Hypre</td><td>MPI+OpenMP w/ OpenACC</td><td>Hybrid (particle used here)</td><td>4x NVIDIA V100 GPUs + 4x IBM Power9 CPU cores</td><td>82 hours</td></tr><tr><td>Laplace<sup>I</sup></td><td>EP-PIC2D</td><td>Fortran</td><td>PARDISO</td><td>MPI+OpenMP</td><td>Particle</td><td>72x Intel Skylake 2.30 GHz CPU cores</td><td>71 hours</td></tr><tr><td>Laplace<sup>V</sup></td><td>LePIC2D</td><td>Fortran</td><td>in house multigrid</td><td>MPI+OpenMP</td><td>Particle</td><td>28x Intel Xeon CPU E5-2690 v4 2.60 GHz</td><td>128 hours</td></tr><tr><td>RUB</td><td>RUBPIC</td><td>C/C++</td><td>NAG geometric multigrid</td><td>CUDA</td><td>Hybrid</td><td>1x NVIDIA V100 GPU + 1x Intel Xeon Gold 5220 CPU core</td><td>80 hours</td></tr><tr><td>Wigner</td><td>MZ/X-PIC</td><td>C/C++</td><td>in house (spectral, direct)</td><td>CUDA</td><td>Particle</td><td>1x NVIDIA P100 GPU + 1x Intel i7 CPU core</td><td>10 days</td></tr><tr><td>DCU</td><td>DCU</td><td>C</td><td>in house multigrid</td><td>OpenMP</td><td>Domain</td><td>32x AMD EPYC 7302 CPU cores</td><td>10 days</td></tr><tr><td>ONERA</td><td>RHEI</td><td>Fortran</td><td>in house (geometric multigrid)</td><td>MPI+OpenMP</td><td>Hybrid</td><td>48x Intel Sapphire Rapids 8468 CPU cores</td><td>16 days</td></tr><tr><td>UC3M</td><td>PICASO</td><td>Fortran</td><td>PARDISO</td><td>OpenMP</td><td>Particle</td><td>40x Intel Xeon Silver 4316 2.30 GHz CPU cores</td><td>7 days 13 hrs</td></tr></tbody></table>
<table><tr><td>UC3M<sup>G</sup></td><td>PICASO-G</td><td>CUDA/Python</td><td>in house (spectral, direct)</td><td>CUDA</td><td>Particle</td><td>1x NVIDIA A100 GPU</td><td>10 hrs</td></tr><tr><td>USask</td><td>EDIPIC-2D</td><td>Fortran</td><td>PETSc +Hypre</td><td>MPI</td><td>Hybrid</td><td>256x Intel E5-2683 v4 2.1 GHz CPU cores</td><td>42 days</td></tr><tr><td>LPP</td><td>LPPic</td><td>Fortran</td><td>Hypre</td><td>MPI</td><td>Domain</td><td>128x AMD Rome 2.6 GHz CPU cores</td><td>55 days</td></tr><tr><td>Stanford</td><td>XPIC2D</td><td>C/C++</td><td>Hypre</td><td>MPI</td><td>Particle</td><td>32x AMD EPYC 7742 2.25 GHz CPU cores</td><td>37 days</td></tr><tr><td>VKI</td><td>PANTERA</td><td>Fortran</td><td>PETSc</td><td>MPI</td><td>Hybrid</td><td>768x Intel Skylake Xeon Platinum 8174 CPU cores</td><td>72 hours</td></tr><tr><td>CERFACS</td><td>AVIP-PIC</td><td>Fortran</td><td>PETSc (GAMG+CG)</td><td>MPI</td><td>Domain</td><td>192x AMD EPYC 9654 (Genoa) 2.4 GHz CPU cores</td><td>7 days 17 hrs</td></tr><tr><td>CNR-ENEA</td><td>PICCOLO</td><td>Fortran</td><td>PETSc + hypre</td><td>MPI+OpenMP</td><td>Hybrid</td><td>42x AMD EPYC 7402 CPU cores</td><td>30 days</td></tr><tr><td>Sandia</td><td>Aleph</td><td>C/C++</td><td>Trilinos (direct solver)</td><td>MPI</td><td>Hybrid</td><td>640x Intel Sandy Bridge 2.5 GHz CPU cores</td><td>14 days</td></tr><tr><td>ICL</td><td>IPPL-2D</td><td>Julia</td><td>Julia direct solver</td><td>MPI</td><td>Domain</td><td>32x AMD EPYC 2.25 GHz CPU cores</td><td>-</td></tr></table>
# 3 Results & discussion

Before comparing simulations across participating groups, we review a few of the key results from a single simulation. Figure 4 shows contour plots of ion density over a single period of the rotating spoke, with an animation of this figure available on YouTube (or full URL available in Ref. [155]). Using the diagnostic described in Section 2.3, the mean period over 13 rotations encompassed by the 300 µs averaging time is found to be 23.1 µs for a corresponding rotation frequency of 43.2 kHz. Another important observation is that outside of the central injection region, the discharge is mostly a vacuum, except when the spoke sweeps by. In Section 3.2 we will discuss how this can cause issues when comparing measurements, such as temperature, between codes.

<figure><img src="image_5.png" /><figcaption>Figure 4: Contour plots of ion density at four different phases of spoke rotation after the system has reached quasi-steady state. Times are given relative to the left most figure.</figcaption></figure>

As discussed in Section 1 spoke formation may be the result of various mechanisms; however, due to the uniform magnetic field and collisionless nature of these simulations, it is hypothesized that the spoke shown here is a saturation of the collisionless Simon-Hoh instability (a type of GDI). Under this assumption, Eq. 4 in Ref. [93] provides a formula for the spoke rotation frequency as a function of the discharge radius, ion mass, radial electric field, and gradient length scale. For this simulation we have a radial electric field of approximately 100 V/m and gradient length scale of 7.1 mm giving a predicted spoke frequency of 53 kHz, which is reasonable compared to the measured frequency.

Figure 4, and associated animation, also demonstrate the complex, chaotic and multiscale nature of the discharge. On the largest scales we clearly have the rotation of the relatively slowly evolving spoke; however, this rotation does not appear to be rigid, and the animation in particular shows the persistent shedding of density at various phases of spoke rotation. Although harder to observe, there are also finer scale structures which exist within the core of the discharge and along the spoke arm, which may be a form of ion-acoustic or electron-cyclotron-drift instabilities [37, 119, 156, 157]. Since this simulation is collisionless, radial transport of magnetized electrons can only be induced by these waves and instabilities and is highly correlated with the spoke rotation. A further discussion or analysis of this transport is beyond the scope of this benchmarking paper, however we refer the reader to many of the excellent references in Section 1 which explore this in depth.

## 3.1 Long time scale simulation comparison

The first simulation comparisons are made between codes which were run for extended times, with the goal of collecting uncertainty bounds for time averaged diagnostics. The list of codes, run times, and associated measured spoke frequencies are shown in Table 3. From 0.2 – 4.0 ms the probe in the PPPL code measures 163 passages of the spoke for a, previously mentioned, mean frequency of 43.2 kHz. The uncertainty bounds are defined by computing averaged frequencies over each group of 13 sequential spoke rotations and then considering the maximum and minimum frequencies from the set. This gives a min/max mean frequency of 41.1 kHz and 46.1 kHz, respectively. Both Laplace codes (EP-PIC2D and LePIC2D) agree with this value within the uncertainty bounds.
<table><caption>Table 3: Details for long time simulations including spoke frequency measurements.</caption><thead><tr><th>Institute</th><th>Code name</th><th>Total time</th><th>Spoke frequency</th><th>Frequency discrepancy</th></tr></thead><tbody><tr><td>PPPL</td><td>LTP-PIC</td><td>4.0 ms</td><td>43.2 (+2.9 -2.1) kHz</td><td>-</td></tr><tr><td>Laplace<sup>I</sup></td><td>EP-PIC2D</td><td>5.0 ms</td><td>43.3 kHz</td><td>0.1 kHz</td></tr><tr><td>Laplace<sup>V</sup></td><td>LePIC2D</td><td>1.6 ms</td><td>42.6 kHz</td><td>0.8 kHz</td></tr></tbody></table>

Figure 5 compares time averaged half-profile slices for (a) ion density, (b) electric potential, and (c) electron temperature (Type I). The uncertainty bounds (shown as gray in the figure) are obtained by considering statistics of a rolling average, over 300 µs, of the PPPL data only. Each average begins at the start of the subsequent snapshot, providing a total of 7,000 averages. The lower and upper bounds are computed as the minimum and maximum values within this set of averages at each grid location. This definition proved critical to assessing agreement within the larger comparison between all codes, as will be discussed in Section 3.2. Using this uncertainty bounds, good agreement is observed between all three codes considered here.

<figure><img src="image_6.png" /><figcaption>Figure 5: Comparison of time-averaged slices for long-time simulations: (a) ion-density, (b) plasma potential, and (c) electron temperature (Type I).</figcaption></figure>
## 3.2 Simulation comparison for all codes

Proceeding with a comparison between all codes, Fig. 6 shows plots of time averaged half-profile slices for (a) ion density, (b) electric potential, and (c) electron temperature (Type I). Note that the results from Section 3.1 for long time scale simulations are reproduced here for the relevant codes. Using the same uncertainty bounds defined in Section 3.1, all codes agree for plots of ion density and plasma potential. Turning to the electron temperature, we observe mixed results, with most codes falling within the tight uncertainty bounds, but with several lying outside of these bounds either at the center or edges of the discharge.

<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison of time-averaged slices between all codes: (a) ion density, (b) plasma potential, and (c) electron temperature (Type I). See supplementary material for the mean and uncertainty bounds data for each plot.</figcaption></figure>

Understanding the cause of this discrepancy in electron temperature across simulations resulted in significant further analysis and a large number of additional simulations. Fundamentally, there are two issues with the benchmark which resulted in this outcome. The first is that the system is highly multiscale, especially with the presence of the low-frequency rotating spoke. Therefore, the proposed 300 µs averaging time is likely too short to achieve convergence between simulations. With consideration to the simulation run times listed in Table 2, it was decided that asking all groups to simulate longer times to improve averaging would be prohibitive. Therefore, this issue was remediated by having a few select groups perform very long time simulations to; (1) check that convergence could be achieved at these time scales and (2) define the uncertainty bounds which could be used to encapsulate accuracy from all codes, the results of which were discussed in Section 3.1. This demonstrated excellent agreement between the three independently developed codes (from PPPL, and two from Laplace) across all tested diagnostics, suggesting that if all codes were run to this length better agreement would be observed.
The second issue arises due to the persistence of vacuum regions throughout the simulation, as shown in Fig. 4. When capturing snapshots of temperature data there can be ambiguity in how to handle cells with few or zero macro-particles. When there are zero macro-particles in a given cell, most groups set the temperature to zero, whereas in reality the temperature is undefined. When averaging over time, these points should be excluded from the procedure, however most codes include the data point as a zero. This would ordinarily not be a problem since most PIC simulations fill the area with particles, however the vacuum regions in this configuration exacerbate the issue. Things become more subtle when considering cells with few, or a fraction of a particle (which is possible for cloud-in-cell shape functions). In particular, if there are only a few particles, the statistics for sampling temperature are poor, leading to significant noise. This noise can lead to anomalously high temperature calculations within a given cell which can easily skew the time averaged results. Different codes also computed temperature in slightly different ways, with some using sample variance and others relying on population variance, which can be related to each other via Bessel's correction [158]. For a well-sampled distribution, these values would agree closely, however in the vacuum region, which are poorly sampled by particles, the discrepancy can be large.

Further issues included whether codes were defined with boundaries at cell-centers or cell-nodes, the type and convergence tolerance of the numerical methods used to solve the Poisson equation, and the effect of the Pseudo-Random-Number-Generator algorithm. However, it is worth noting that sensitivity tests with the PPPL and the two Laplace codes revealed that neither of these factors played a strong role in influencing the temperature profiles.

To remedy this discrepancy between codes, we asked several of the “outlier” code authors to rerun their simulations with the alternative temperature diagnostic described by Eq. (7) (*Temperature Type II*). Averaging velocity moments prior to calculating temperature has several advantages. This includes that velocity moments are truly zero when there are no particles in the cell, and time-averaged moments gather significant statistics to end up being non-zero everywhere (except at the boundaries). Furthermore, it reduces the chance of small values being present in the denominator of the temperature calculation, further improving statistical accuracy. While this is not strictly the standard definition of time-averaged temperature, it is a valid way to compare numerical methods if done consistently. Figure 7 shows plots of *Temperature Type II* computed via Eq. (7), demonstrating agreement between nearly all previously outlier codes within the uncertainty bounds.

<figure><img src="image_8.png" /><figcaption>Figure 7: Comparison of time-averaged slices of electron temperature (Type II) for select codes. See supplementary material for the mean and uncertainty bounds data of the plot.</figcaption></figure>

While investigating the causes of discrepancy between groups, numerous simulations were run with both the Sandia and CERFACS codes in an attempt to understand how the unstructured nature of their meshes could influence the results of the benchmark. Although no clear influence was identified, suspected impacts include that (1) the location of the triangular cell center may not align with those from a Cartesian grid, and (2) the smaller area of the triangular cell may lead to an increase in measured fluctuations in regions of low statistical resolution. The Aleph code resolved outstanding
discrepancies by printing complete particle phase space data at each time step and mapping this data onto a Cartesian grid in post-processing. Despite best attempts, the Cerfacs code still fell around 1.3% outside of the Temperature Type II measurement uncertainty bounds near the center of the simulation. We suspect this is due to the unstructured nature of the code since it has been well benchmarked in previous Landmark efforts [24, 25].

The newly developed UC3M<sup>G</sup> code also showed some minor discrepancies close to the uncertainty bounds for the density and temperature (Type I) measurements. After extensive investigation, this discrepancy was conclusively attributed to the group's use of a spectral field solver, which showed root-mean-square fluctuations 3-15% in the electric field larger than those from a finite-difference field solver.

Finally, Fig. 8 compares the spoke rotation frequency measured by all groups, showing agreement within the measured frequency bounds of [41.1, 46.1] kHz from long time simulations using the PPPL code. Despite the challenges associated with this benchmark, we believe that the results presented in this section constitute sufficient evidence to conclude that there is agreement between almost all codes across a variety of diagnostics.

<figure><img src="image_9.png" /><figcaption>Figure 8: Comparison of the spoke rotation frequency, showing agreement between all codes. The reference frequency of 43.2 kHz is produced from long time simulations with the PPPL code with uncertainty bounds of -2.1 kHz and +2.9 kHz delineated by the red dashed lines (see Section 3.1 for further details).</figcaption></figure>

## 3.3 Lessons learned for benchmark selection

Despite the best efforts of the authors, there were several significant challenges associated with completing this benchmarking comparison, some of which are discussed in Sec. 3.2. In this section, we aim to highlight some of the lessons learned from this process with the hope that future benchmarking efforts can avoid the same pitfalls.

**Benchmark lesson 1: Prepare a well defined benchmark description.** Based on the author's experience participating in previous PIC benchmarks, a best attempt was made to prepare a clear, concise, and sufficient description of the Penning discharge setup as well as the diagnostics used to compare codes. Despite this effort, several important details were missing, which in hindsight would have proved useful. These included an appropriate definition of how to handle temperature in vacuum regions and the specific equations for calculating radial injection locations on the Cartesian grid. The choice of the pseudo-random number generator was also a question raised at code author meetings, although it was ultimately demonstrated that this had little influence on the results. Nonetheless, we emphasize that preparing a suitable benchmark description is paramount to success in any such effort.

**Benchmark lesson 2: Pick a low cost benchmark.** A major challenge associated with obtaining convergence between all codes was the large simulation cost, both in computing resources and time, of modeling the Penning discharge (see Table 2). This made it challenging for many groups to
rerun their simulations or run for longer times, which proved essential for demonstrating agreement. Of course, a benchmark should run for sufficient time such that it tests the capabilities of each code to reproduce important behavior relevant to the domain of physics, but ideally it should not be longer than that. This is particularly challenging in the LTP community, where the phenomena under investigation are often highly multi-scale, especially when collective effects and collisions are included. Figure 9 shows how the spectrum of the total number of ions in the simulation follows a power law, indicative of chaotic behavior [159] and the highly multiscale nature of the discharge, making it less ideal as a benchmark configuration. In conclusion, a benchmark should be selected that is long enough to capture the physics being tested, but ideally not longer.

<figure><img src="image_10.png" /><figcaption>Figure 9: (a) Plot of total number of ions within the simulation vs time, and (b) frequency spectrum of the data from plot (a), demonstrating the chaotic nature of the Penning discharge benchmark configuration.</figcaption></figure>

**Benchmark lesson 3: Avoid systems which have regions of poor statistic resolution.**

This lesson pertains particularly to particle-based methods, such as the PIC algorithm, which is an approximate solution to the Vlasov equation. The discussion in Section 3.2 highlights that poor statistical sampling within the vacuum regions of the Penning discharge led to significant discrepancies in temperature measurements between the codes. In general, poor statistical sampling can lead to high, and unpredictable levels of noise which can easily skew results. It is therefore recommended to avoid simulations which consist of these poorly resolved regions. As an example, the *Landmark* 2a [24] and 2b [25] benchmark configurations provide a more appropriate setup which avoids these pitfalls.

**Benchmark lesson 4: More data is better.** Even with a well prepared benchmark description, it is impossible for the researchers to predict all possible ways in which disagreement may arise between codes. To prepare for this eventuality, it is wise to request that participants collect more than the minimum data required to compare results. Case in point, if participants in this benchmarking effort had been requested to record velocity moments, rather than temperature, it would have allowed for a rapid comparison between the Type I and Type II temperature diagnostics, eliminating the need for further simulations. Additionally, collecting more complete data can make it easier to determine
uncertainty bounds for each simulation, which can further simplify comparison.

Despite the significant time that went into resolving issues with this benchmark, the authors want to emphasize that such efforts will continue to be critical for ensuring that complex codes can accurately reproduce results which are relevant to the broader community. More complex benchmark configurations, such as those with collisions, are also required, and the ultimate goal of the community should be to perform rigorous validation studies against experimental data. The lessons learned here will be essential to achieving that goal.

## 3.4 High performance computing lessons learned

Although the original intent of this benchmarking effort was not to assess or even stress test the computational performance of the participating codes, the results shown in Table 2 (comparing code implementations and performance) reveal some interesting trends that pertain to code implementation and high-performance computing (HPC) strategies.

**HPC Lesson 1: Parallel decomposition strategy matters.** One trend which became apparent throughout the benchmark is the total runtime for a given code compared to the parallel decomposition strategy. For PIC simulations, there are generally two approaches. The simplest, particle decomposition, is where the list of particles in memory is split (usually evenly) between MPI tasks, each of which owns a copy of the entire simulation grid. Particle decomposition is simple to implement and easy to load balance, however the downside of each task owning an entire copy of the grid is the requirement for at least two MPI Allreduce calls at each time step to share grid information. Therefore, most codes implement a domain decomposition approach, whereby each task owns a subset of the global mesh and all particles which reside on that mesh. This avoids the need for an expensive global MPI call, and improves scaling performance, however introduces code complexity due to particle communication across domain boundaries.

The Penning discharge configuration of this benchmark also highlights the second challenge with domain decomposition, that of load balancing. As shown in Fig. 4, the spoke, and associated simulation particles, comprise only a narrow angular region of the domain outside of the central injection region. Therefore, at any given time, a uniformly domain decomposed grid will have extremely poor load balancing, effectively restricting code performance to that of a few cores handling a majority of the particles. By comparing the “Decomposition” and “Runtime” columns of Table 2, it is clear that codes which implement particle decomposition performed, on average, better than those with domain decomposition. We note that the situation may improve for collisional simulations, where scattering can better distribute plasma density throughout the domain.

Despite this trend, the value of domain decomposition cannot be understated. For large simulations, particle decomposition becomes increasingly expensive and eventually intractable due to communication bottlenecks and memory limitations. However, the lesson learned here is that a hybrid approach to parallel decomposition can improve versatility of code performance across a wide range of applications. Such an approach has two levels of task parallelism, with the outer (higher) level being domain decomposition and then the inner (lower) level being particle decomposition within each domain. This flexibility allows a code to better handle load balancing when a system has large plasma gradients, yet still achieve reasonable parallel scaling. Alternatively, one could consider dynamic load balancing [160–162] where the domain owned by each process is adjusted on the fly based on the spatial density of the simulation particles.

**HPC Lesson 2: GPUs are effective at accelerating PIC codes.** A second trend which emerges from Table 2 relates to the hardware targeted by each code. In general, codes which are able to take advantage of GPU acceleration perform faster than those using CPUs. In the last decade, general purpose graphical processing units have provided an affordable, massively parallel alternative to multi-core CPUs. At the time of writing, these chips have a core count on the order of 10,000, and thread counts over 1 million for each GPU device, making them ideal for acting on large chunks of contiguous memory, such as particle phase-space data. Furthermore, the implementation of efficient atomic operations in languages such as CUDA (which can be targeted by abstraction layers such as OpenMP and OpenACC) leads to very fast interpolation routines. GPUs also benefit from large caches of high speed memory, which can store the entire local grid chunk, further accelerating interpolation and grid based routines. Details on writing low-temperature plasma PIC codes for GPUs can be found in Refs [9, 19, 163].
**HPC Lesson 3: Adhering to HPC best practices is essential.** Even when optimizing over parallel strategies and chip architectures, it is clear from Table 2 that an enormous amount of time is required for each Penning discharge simulation. This results from a familiar problem to all computational plasma physicists, the vast range of time scales exhibited by plasma phenomena, which are particularly present in simulations with low-frequency spoke oscillations. Although in the previous section, we recommended avoiding such long time scales in future benchmarks, this obviously does not eliminate the need to model such waves and instabilities for future scientific and engineering purposes. If groups are to continue relying on the standard explicit PIC method, then high-performance computing best practices are essential to reduce the computing time per simulation step, and therefore total simulation time. Details on how to achieve these have been documented extensively [9, 11, 19, 164].

This lesson also extends to selecting external packages which will be integrated into a code. For LTP simulations, this is usually in the form of the linear algebra package chosen (or written in-house) to solve the Poisson equation. The choice of solver can play an important role, with direct solvers being best for small domains and Krylov and/or multigrid methods being best for larger domains [165].

Broadly, this lesson comes down to understanding how to efficiently map the fundamental PIC algorithms onto relevant (yet constantly changing) computational architecture for maximum performance. With a trend towards increasing heterogeneity in computer systems, this requirement will only continue to grow in importance over the coming years.

**HPC Lesson 4: Consider alternative algorithms.** Code developers could also consider alternative algorithms which can step over the usually strict numerical resolution requirements of PIC, these include energy-conserving PIC, which allow for larger cell sizes [166–168], as well as implicit [169–171] or semi-implicit [172–174] time-stepping techniques, allowing for larger time steps. Advanced PIC algorithms, which may offer further performance improvements, include sparse-grid PIC [175–179] or reduced-order PIC [180–182], both of which have proven accurate for modeling E×B discharges. Of course, for accurate simulation results, care must always be taken to resolve the appropriate time and length scales of the relevant physical phenomena which give rise to the spoke and associated anomalous transport.

We also wish to highlight some of the expected performance challenges associated with the unstructured Sandia and CERFACS codes. In general, the handling of particle data in memory, and transfer of data between adjacent cells can incur an increased overhead. Furthermore, when particle phase data is stored contiguously there is an increased cost to identifying the parent cell, which is otherwise a trivial process on a Cartesian grid. Finally, unstructured meshes are unable to take advantage of the highest performing linear algebra preconditioners such as geometric multigrid. Although an unstructured grid can lead to increased simulation cost, the benefit is significantly increased flexibility, especially when modeling complex engineering systems.

# 4 Conclusion

The computational study of low-temperature plasmas is essential to future progress in plasma science and engineering. A critical part of this effort is the verification, or benchmarking, of codes used to simulate this phenomenon. The two-dimensional Penning discharge configuration detailed in this article has many attributes which make it an ideal stress test of the required physics modeling and computational performance of low-temperature plasma kinetic simulation software. Unlike the *Landmark 2* series of benchmarks, this new configuration spans an order of magnitude larger range of time scales, and captures the emergence of the large scale rotating spoke structure. Seventeen codes, a record for this community, were put to the test on this benchmark, with agreement found for nearly all across four important measurements.

Several challenges were encountered throughout this process, resulting from a lack of detail in the benchmark description, high computational cost of the configuration, and poor statistical resolution in regions of the discharge. The effort also highlighted important computational considerations for particle-in-cell code design, including flexibility of parallelization strategy, adherence to high-performance computing best practices, and taking advantage of modern computing hardware. The authors hope that describing these lessons will help this, and other simulation communities, design improved benchmarks for their respective software.

Despite a growing number of benchmarks within this community, new configurations are required to test more complex low-temperature plasma physics, including collisions, chemical reactions, and
surface processes. Furthermore, validation against experimental data is severely lacking, and essential
to the future use of these codes for modeling complex engineering systems. We hope that this effort
serves as a small step along the road towards developing robust, efficient, and useful plasma simulation
software.

Acknowledgments

This research was supported by the U.S. Department of Energy through the PPPL CRADA agreement with Applied Materials, the Princeton Collaborative Research Facility (PCRF) under contract No. DE-AC02-09CH11466, and the Laboratory Directed Research and Development (LDRD) program. The United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purpose(s).

The work of Federico Petronio was partially supported by Agence de l'Innovation de Defense – AID - via Centre Interdisciplinaire d'Etudes pour la Défense et la Sécurité – CIEDS - (project 2023 – validHETion). This project was provided with computer and storage resources by GENCI at TGCC thanks to the grant 2024-A0160510439 on the supercomputer Joliot Curie's Irene-ROME partition.

For what concerns the CNR-ENEA study, simulations with the PICCOLO code were carried out
with Eni HPC4 Supercomputing Cluster, under the Joint Research Agreement ENI-CNR on 'Fusione
a Confinamento Magnetico'.

The Wigner Research Center for Physics would like to thank the Hungarian Office for Research, Development and Innovations for their support through grant number 152843.

This work was supported by DOE grant No. DE-SC0022201. This work used the capabilities of the SNL Plasma Research Facility, supported by DOE SC FES. This article has been authored by an employee of National Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee owns all right, title and interest in and to the article and is solely responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan https://www.energy.gov/downloads/doi-public-access-plan.

The work of Pietro Parodi is funded by an FWO Strategic Basic PhD fellowship (reference 1S24022N). The authors of the VKI group gratefully acknowledge the Gauss Centre for Supercomputing for providing computing time on SuperMUC-NG at Leibniz Supercomputing Centre through the project "Heat flux regulation by collisionless processes in heliospheric plasmas – ARIEL".

This work conducted by the Ruhr University Bochum group was funded by the German Research Foundation within the framework of the project "Modellierung und Simulation von Hochleistungsmagnetonentladungen" (No. 550860775).

The contribution of E. Ahedo, M.P. Encinar, P. Fajardo, A. Marín-Cebrián has been supported by the R&D project PID2022-140035OB-I00 (HEEP) funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. The contribution of E. Bello-Benítez, M. Merino has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Starting Grant ZARATHUSTRA, grant agreement No 950466).

This work conducted at Laplace was granted access to the HPC resources of CALMIP supercomputing center under Allocation No. 2013-P1125.

G. Bogopolsky acknowledges financial support from Safran Spacecraft Propulsion and the Association Nationale de la Recherche et de la Technologie (ANRT) as part of a CIFRE convention.

References

[1] I. Adamovich, S. Agarwal, E. Ahedo, L. L. Alves, S. Baalrud, N. Babaeva, A. Bogaerts, A. Bourdon, P. J. Bruggeman, C. Canal, et al., "The 2022 plasma roadmap: Low temperature plasma science and technology," *Journal of Physics D: Applied Physics*, vol. 55, no. 373001, 2022.
[2] K. J. Kanarik, "Inside the mysterious world of plasma: A process engineer's perspective," *Journal of Vacuum Science & Technology A*, vol. 38, no. 031004, 2020.

[3] A. Anders, "Tutorial: Reactive high power impulse magnetron sputtering (R-HIPIMS)," *Journal of Applied Physics*, vol. 121, no. 171101, 2017.

[4] W. Villafana, A. T. Powis, S. Sharma, I. D. Kaganovich, and A. V. Khrabrov, "Establishing criteria for the transition from kinetic to fluid modeling in hollow cathode analysis," *Physics of Plasmas*, vol. 31, no. 093504, 2024.

[5] D. M. Goebel, R. Poeschel, and R. Schumacher, "Low voltage drop plasma switch for inverter and modulator applications," *Review of scientific instruments*, vol. 64, pp. 2312-2319, 1993.

[6] C. Charles, "Plasmas for spacecraft propulsion," *Journal of Physics D: Applied Physics*, vol. 42, no. 163001, 2009.

[7] E. Ahedo, "Plasmas for space propulsion," *Plasma Physics and Controlled Fusion*, vol. 53, no. 124037, 2011.

[8] W. T. Osowiecki, M. J. Coogans, S. Sriraman, R. Ranjan, Y. Lu, and D. M. Fried, "Achieving sustainability in the semiconductor industry: The impact of simulation and AI," *IEEE Transactions on Semiconductor Manufacturing*, vol. 37, pp. 464-474, 2024.

[9] Z. Juhasz, J. Ďurian, A. Derzsi, Š. Matejčík, Z. Donkó, and P. Hartmann, "Efficient GPU implementation of the Particle-in-cell/Monte-Carlo collisions method for 1D simulation of low-pressure capacitively coupled plasmas," *Computer Physics Communications*, vol. 263, no. 107913, 2021.

[10] F. Faraji and M. Reza, "Machine learning applications to computational plasma physics and reduced-order plasma modeling: a perspective," *Journal of Physics D: Applied Physics*, vol. 58, no. 102002, 2025.

[11] I. D. Kaganovich, A. Smolyakov, Y. Raitses, E. Ahedo, I. G. Mikellides, B. Jorns, F. Taccogna, R. Gueroult, S. Tsikata, A. Bourdon, et al., "Physics of E×B discharges relevant to plasma propulsion and similar technologies," *Physics of Plasmas*, vol. 27, no. 120601, 2020.

[12] R. W. Hockney and J. W. Eastwood, *Computer Simulation Using Particles*. New York, NY: CRC Press, 2021.

[13] C. K. Birdsall and A. B. Langdon, *Plasma physics via computer simulation*. New York, NY: CRC press, 2018.

[14] T. Tajima, *Computational plasma physics: with applications to fusion and astrophysics*. New York, NY: CRC press, 2004.

[15] V. Vahedi and M. Surendra, "A Monte Carlo collision model for the particle-in-cell method: applications to argon and oxygen discharges," *Computer Physics Communications*, vol. 87, no. 1-2, pp. 179-198, 1995.

[16] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Ejkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, and H. Zhang, "PETSc Web page." http://www.mcs.anl.gov/petsc, 2015.

[17] "hypre: High performance preconditioners." https://llnl.gov/casc/hypre, https://github.com/hypre-space/hypre.

[18] T. Trilinos Project Team, *The Trilinos Project Website*. https://trilinos.github.io.

[19] A. T. Powis, *Particle Methods for Modeling Magnetospheric Diagnostics and Low-Temperature Plasma Physics*. PhD thesis, Princeton University, 2021.

[20] P. Mertmann, D. Eremin, T. Mussenbrock, R. P. Brinkmann, and P. Awakowicz, "Fine-sorting one-dimensional particle-in-cell algorithm with Monte-Carlo collisions on a graphics processing unit," *Computer Physics Communications*, vol. 182, pp. 2161-2167, 2011.
[21] P. Tranquilli, L. Ricketson, and L. Chacón, “A deterministic verification strategy for electrostatic particle-in-cell algorithms in arbitrary spatial dimensions using the method of manufactured solutions,” *Journal of Computational Physics*, vol. 448, no. 110751, 2022.

[22] M. M. Turner, A. Derzsi, Z. Donko, D. Eremin, S. J. Kelly, T. Lafleur, and T. Mussenbrock, “Simulation benchmarks for low-pressure plasmas: Capacitive discharges,” *Physics of Plasmas*, vol. 20, no. 013507, 2013.

[23] A. Smolyakov, T. Zintel, L. Couedel, D. Sydorenko, A. Umnov, E. Sorokina, and N. Marusov, “Anomalous electron transport in one-dimensional electron cyclotron drift turbulence,” *Plasma Physics Reports*, vol. 46, pp. 496–505, 2020.

[24] T. Charoy, J.-P. Boeuf, A. Bourdon, J. A. Carlsson, P. Chabert, B. Cuenot, D. Eremin, L. Garrigues, K. Hara, I. D. Kaganovich, et al., “2D axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas,” *Plasma Sources Science and Technology*, vol. 28, no. 105010, 2019.

[25] W. Villafana, F. Petronio, A. Denig, M. Jimenez, D. Eremin, L. Garrigues, F. Taccogna, A. Alvarez-Laguna, J.-P. Boeuf, A. Bourdon, et al., “2D radial-azimuthal particle-in-cell benchmark for E×B discharges,” *Plasma Sources Science and Technology*, vol. 30, no. 075002, 2021.

[26] A. Morozov, “Effect of near-wall conductivity in a strongly magnetized plasma,” *Zh. Prikl. Mekh. Tekh. Fiz.*, vol. 3, p. 19, 1968.

[27] A. Morozov and V. Savel'Ev, “Theory of the near-wall conductivity,” *Plasma Physics Reports*, vol. 27, pp. 570–575, 2001.

[28] E. Y. Choueiri, “Plasma oscillations in Hall thrusters,” *Physics of Plasmas*, vol. 8, pp. 1411–1426, 2001.

[29] E. Chesta, C. M. Lam, N. B. Meezan, D. P. Schmidt, and M. A. Cappelli, “A characterization of plasma fluctuations within a Hall discharge,” *IEEE transactions on plasma science*, vol. 29, pp. 582–591, 2001.

[30] J. Adam, A. Héron, and G. Laval, “Study of stationary plasma thrusters using two-dimensional fully kinetic simulations,” *Physics of Plasmas*, vol. 11, pp. 295–305, 2004.

[31] F. Taccogna, S. Longo, M. Capitelli, and R. Schneider, “Anomalous transport induced by sheath instability in Hall effect thrusters,” *Applied Physics Letters*, vol. 94, no. 251502, 2009.

[32] A. Héron and J. Adam, “Anomalous conductivity in Hall thrusters: Effects of the non-linear coupling of the electron-cyclotron drift instability with secondary electron emission of the walls,” *Physics of Plasmas*, vol. 20, no. 082313, 2013.

[33] T. Lafleur, S. Baalrud, and P. Chabert, “Theory for the anomalous electron transport in Hall effect thrusters. II. Kinetic model,” *Physics of Plasmas*, vol. 23, no. 053503, 2016.

[34] T. Lafleur and P. Chabert, “The role of instability-enhanced friction on ‘anomalous’ electron and ion transport in Hall-effect thrusters,” *Plasma Sources Science and Technology*, vol. 27, no. 015003, 2017.

[35] C. Boniface, L. Garrigues, G. Hagelaar, J. Boeuf, D. Gawron, and S. Mazouffre, “Anomalous cross field electron transport in a Hall effect thruster,” *Applied Physics Letters*, vol. 89, no. 161503, 2006.

[36] F. Taccogna, P. Minelli, Z. Asadi, and G. Bogopolsky, “Numerical studies of the E×B electron drift instability in Hall thrusters,” *Plasma Sources Science and Technology*, vol. 28, no. 064002, 2019.

[37] T. Charoy, T. Lafleur, A. Tavant, P. Chabert, and A. Bourdon, “A comparison between kinetic theory and particle-in-cell simulations of anomalous electron transport in E×B plasma discharges,” *Physics of Plasmas*, vol. 27, no. 063510, 2020.
[38] Z. A. Brown and B. A. Jorns, “Anomalous cross-field transport in a Hall thruster inferred from direct measurement of instability growth rates,” *Physical Review E*, vol. 108, no. 065204, 2023.

[39] J.-P. Boeuf, “Tutorial: Physics and modeling of Hall thrusters,” *Journal of Applied Physics*, vol. 121, no. 011101, 2017.

[40] G. Janes and R. Lowder, “Anomalous electron diffusion and ion acceleration in a low-density plasma,” *The Physics of Fluids*, vol. 9, pp. 1115–1123, 1966.

[41] Y. B. Esipchuk, A. Morozov, G. Tilinin, and A. Trofimov, “Plasma oscillations in closed-drift accelerators with an extended acceleration zone,” *Soviet Physics Technical Physics*, vol. 18, no. 928, 1974.

[42] J. Parker, Y. Raitses, and N. Fisch, “Transition in electron transport in a cylindrical Hall thruster,” *Applied Physics Letters*, vol. 97, no. 091501, 2010.

[43] C. L. Ellison, Y. Raitses, and N. J. Fisch, “Direct measurement of spoke-induced, cross-field electron transport in a cylindrical Hall thruster,” in *32nd International Electric Propulsion Conference, Weisbaden, Germany*, 2011.

[44] M. S. McDonald and A. D. Gallimore, “Rotating spoke instabilities in Hall thrusters,” *IEEE Transactions on Plasma Science*, vol. 39, pp. 2952–2953, 2011.

[45] C. Ellison, Y. Raitses, and N. J. Fisch, “Cross-field electron transport induced by a rotating spoke in a cylindrical Hall thruster,” *Physics of Plasmas*, vol. 19, no. 013503, 2012.

[46] M. E. Griswold, C. Ellison, Y. Raitses, and N. Fisch, “Feedback control of an azimuthal oscillation in the E×B discharge of Hall thrusters,” *Physics of Plasmas*, vol. 19, no. 053506, 2012.

[47] Y. Raitses, M. Griswold, L. Ellison, J. Parker, and N. Fisch, “Studies of rotating spoke oscillations in cylindrical Hall thrusters,” in *48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, Atlanta, Georgia*, 2012.

[48] D. Liu, R. E. Huffman, R. D. Branam, and W. A. Hargus, “Ultrahigh speed images of Hall thruster azimuthal instabilities,” *IEEE Transactions on Plasma Science*, vol. 42, pp. 2656–2657, 2014.

[49] M. J. Sekerak, B. W. Longmier, A. D. Gallimore, D. L. Brown, R. R. Hofer, and J. E. Polk, “Azimuthal spoke propagation in Hall effect thrusters,” *IEEE Transactions on Plasma Science*, vol. 43, pp. 72–85, 2014.

[50] M. Cappelli, Y. Raitses, and I. Kaganovich, “Coherent structures in plasmas relevant to electric propulsion,” Tech. Rep. VA-TR-2016-0229, AFRL-AFOSR, 2016.

[51] N. Hyatt, “Hall effect thruster characterization through potential, magnetic, and optical measurements,” in *55th AIAA Aerospace Sciences Meeting, Grapevine, Texas*, 2017.

[52] S. Mazouffre, L. Grimaud, S. Tsikata, K. Matyash, and R. Schneider, “Rotating spoke instabilities in a wall-less Hall thruster: Experiments,” *Plasma Sources Science and Technology*, vol. 28, no. 054002, 2019.

[53] A. Anders, P. Ni, and A. Rauch, “Drifting localization of ionization runaway: Unraveling the nature of anomalous transport in high power impulse magnetron sputtering,” *Journal of Applied Physics*, vol. 111, no. 053304, 2012.

[54] J. Winter, A. Hecimovic, T. de los Arcos, M. Böke, and V. Schulz-von der Gathen, “Instabilities in high-power impulse magnetron plasmas: from stochasticity to periodicity,” *Journal of Physics D: Applied Physics*, vol. 46, no. 084007, 2013.

[55] A. Hecimovic, V. Schulz-von Der Gathen, M. Böke, A. Von Keudell, and J. Winter, “Spoke transitions in HiPIMS discharges,” *Plasma Sources Science and Technology*, vol. 24, no. 045005, 2015.
[56] P. Poolcharuansin, F. L. Estrin, and J. W. Bradley, "The use of segmented cathodes to determine the spoke current density distribution in high power impulse magnetron sputtering plasmas," *Journal of Applied Physics*, vol. 117, no. 163304, 2015.

[57] A. Anders and Y. Yang, "Direct observation of spoke evolution in magnetron sputtering," *Applied Physics Letters*, vol. 111, no. 064103, 2017.

[58] A. Hecimovic, N. Britun, S. Konstantinidis, and R. Snyders, "Sputtering process in the presence of plasma self-organization," *Applied Physics Letters*, vol. 110, no. 014103, 2017.

[59] M. Panjan and A. Anders, "Plasma potential of a moving ionization zone in DC magnetron sputtering," *Journal of Applied Physics*, vol. 121, no. 063302, 2017.

[60] J. Hnilica, P. Klein, M. Šlapanská, M. Fekete, and P. Vašina, "Effect of magnetic field on spoke behaviour in HiPIMS plasma," *Journal of Physics D: Applied Physics*, vol. 51, no. 095204, 2018.

[61] A. P. Ehiasarian, A. Hecimovic, T. de los Arcos, R. New, V. Schulz-von der Gathen, M. Böke, and J. Winter, "High power impulse magnetron sputtering discharges: Instabilities and plasma self-organization," *Applied Physics Letters*, vol. 100, no. 114101, 2012.

[62] N. Brenning, D. Lundin, T. Minea, C. Costin, and C. Vitelaru, "Spokes and charged particle transport in HiPIMS magnetrons," *Journal of Physics D: Applied Physics*, vol. 46, no. 084005, 2013.

[63] K. Thomassen, "Turbulent diffusion in a Penning-type discharge," *Physics of Fluids*, vol. 9, pp. 1836-1842, 1966.

[64] Y. Sakawa, C. Joshi, P. Kaw, F. Chen, and V. Jain, "Excitation of the modified Simon-Hoh instability in an electron beam produced plasma," *Physics of Fluids B: Plasma Physics*, vol. 5, pp. 1681-1694, 1993.

[65] Y. Raitses, I. Kaganovich, and A. Smolyakov, "Effects of the gas pressure on low frequency oscillations in E×B discharges," in *Proceedings of the Joint Conference of 30th ISTS, 34th IEPC, and 6th NSAT, Kobe, Japan*, 2015.

[66] N. Claire, A. Escarguel, C. Rebont, and F. Doveil, "Ion velocity analysis of rotating structures in a magnetic linear plasma device," *Physics of Plasmas*, vol. 25, no. 061203, 2018.

[67] T. Pierre, A. Escarguel, D. Guyomarc'h, R. Barni, and C. Riccardi, "Radial convection of plasma structures in a turbulent rotating magnetized-plasma column," *Physical Review Letters*, vol. 92, no. 065004, 2004.

[68] M. Matsukuma, T. Pierre, A. Escarguel, D. Guyomarc'h, G. Leclert, F. Brochard, E. Gravier, and Y. Kawai, "Spatiotemporal structure of low frequency waves in a magnetized plasma device," *Physics Letters A*, vol. 314, no. 1-2, pp. 163-167, 2003.

[69] S. Jaeger, T. Pierre, and C. Rebont, "Direct observation of a cross-field current-carrying plasma rotating around an unstable magnetized plasma column," *Physics of plasmas*, vol. 16, no. 022304, 2009.

[70] C. Rebont, N. Claire, T. Pierre, and F. Doveil, "Ion velocity distribution function investigated inside an unstable magnetized plasma exhibiting a rotating nonlinear structure," *Physical Review Letters*, vol. 106, no. 225006, 2011.

[71] A. Escarguel, "Optical diagnostics of a low frequency instability rotating around a magnetized plasma column," *The European Physical Journal D*, vol. 56, pp. 209-214, 2010.

[72] J. Y. Kim, J. Y. Jang, J. Choi, J.-i. Wang, W. I. Jeong, M. Elgarhy, G. Go, K.-J. Chung, and Y. Hwang, "Magnetic confinement and instability in partially magnetized plasma," *Plasma Sources Science and Technology*, vol. 30, no. 025011, 2021.
[73] E. Rodríguez, V. Skoutnev, Y. Raitses, A. Powis, I. Kaganovich, and A. Smolyakov, "Boundary-induced effect on the spoke-like activity in E×B plasma," *Physics of Plasmas*, vol. 26, no. 053503, 2019.

[74] C. Cheon, J. Choi, J. B.-W. Koo, and J. Y. Kim, "Experimental evidence of various mode numbers of azimuthal waves in an E × B Penning source for semiconductor processing," *Plasma Sources Science and Technology*, vol. 32, no. 07LT01, 2023.

[75] J. Y. Kim, J. Choi, Y. Hwang, and K.-J. Chung, "High-magnetic-confinement mode in partially magnetized E×B plasmas," *Physical Review E*, vol. 105, no. L053202, 2022.

[76] M. Lee, C. Cheon, J. Choi, H. J. Lee, Y. Hwang, K.-J. Chung, and J. Y. Kim, "Azimuthal and axial dispersion relations of low-frequency oscillations in partially magnetized plasma with crossed electric and magnetic fields," *Physics of Plasmas*, vol. 30, no. 103502, 2023.

[77] R. Przybocki and M. Cappelli, "Comparison of spoke dynamics in an obstructed direct-current magnetron operating on a range of noble gases," *Plasma Sources Science and Technology*, vol. 34, no. 015008, 2025.

[78] M. Valinurov, R. Usmanov, A. Gavrikov, A. Oiler, and V. Smirnov, "Study of rotating "spoke" structure in reflex discharge with thermionic cathode," *Physics Letters A*, vol. 534, no. 130242, 2025.

[79] V. Skoutnev, P. Dourbal, E. Rodriguez, and Y. Raitses, "Fast sweeping probe system for characterization of spokes in E×B discharges," *Review of Scientific Instruments*, vol. 89, no. 123501, 2018.

[80] J. Choi, C. Cheon, Y. Hwang, K.-J. Chung, and J. Y. Kim, "Rotating spoke frequency in E×B Penning source with two-ion-species plasma," *Plasma Sources Science and Technology*, vol. 32, no. 01LT01, 2023.

[81] A. Pushkarev, X. Zhu, S. Polisadov, P. Tang, Z. Yang, and M. Lei, "Spoke-type structures in an ion diode with magnetic insulation of electrons," *Plasma Physics Reports*, vol. 50, pp. 1238-1249, 2024.

[82] Y. Lv, L.-Q. Wei, W.-B. Li, W.-C. Gu, F. Wu, H. Li, Y.-J. Ding, D.-R. Yu, S. Wang, and N. Guo, "Rotating spoke analysis of Hall thruster performance during the extinguishing process," *Acta Astronautica*, vol. 173, pp. 378-382, 2020.

[83] T. Ito, C. V. Young, and M. A. Cappelli, "Self-organization in planar magnetron microdischarge plasmas," *Applied Physics Letters*, vol. 106, no. 254104, 2015.

[84] J. Held, M. George, and A. von Keudell, "Spoke-resolved electron density, temperature and potential in direct current magnetron sputtering and HiPIMS discharges," *Plasma Sources Science and Technology*, vol. 31, no. 085013, 2022.

[85] M. Panjan, S. Loquai, J. E. Klemberg-Sapieha, and L. Martinu, "Non-uniform plasma distribution in DC magnetron sputtering: origin, shape and structuring of spokes," *Plasma Sources Science and Technology*, vol. 24, no. 065010, 2015.

[86] M. J. Sekerak, A. D. Gallimore, D. L. Brown, R. R. Hofer, and J. E. Polk, "Mode transitions in Hall-effect thrusters induced by variable magnetic field strength," *Journal of Propulsion and Power*, vol. 32, pp. 903-917, 2016.

[87] A. Guglielmi, F. Gaboriau, and J. Boeuf, "Simultaneous measurements of axial motion and azimuthal rotation of non-uniformities ("spokes") in a Hall thruster," *Physics of Plasmas*, vol. 29, no. 112108, 2022.

[88] Y. Yang, J. Liu, L. Liu, and A. Anders, "Propagation direction reversal of ionization zones in the transition between high and low current magnetron sputtering," *Applied Physics Letters*, vol. 105, no. 254101, 2014.
[89] A. Hecimovic, M. Böke, and J. Winter, “The characteristic shape of emission profiles of plasma spokes in HiPIMS: the role of secondary electrons,” *Journal of Physics D: Applied Physics*, vol. 47, no. 102003, 2014.

[90] A. Hecimovic and A. Von Keudell, “Spokes in high power impulse magnetron sputtering plasmas,” *Journal of Physics D: Applied Physics*, vol. 51, no. 453001, 2018.

[91] A. Anders, “Self-organization and self-limitation in high power impulse magnetron sputtering,” *Applied Physics Letters*, vol. 100, no. 224104, 2012.

[92] A. Hecimovic, C. Maszl, V. Schulz-Von Der Gathen, M. Böke, and A. Von Keudell, “Spoke rotation reversal in magnetron discharges of aluminium, chromium and titanium,” *Plasma Sources Science and Technology*, vol. 25, no. 035001, 2016.

[93] A. T. Powis, J. A. Carlsson, I. D. Kaganovich, Y. Raitses, and A. Smolyakov, “Scaling of spoke rotation frequency within a Penning discharge,” *Physics of Plasmas*, vol. 25, no. 072110, 2018.

[94] J. Carlsson, I. Kaganovich, A. Powis, Y. Raitses, I. Romadanov, and A. Smolyakov, “Particle-in-cell simulations of anomalous transport in a Penning discharge,” *Physics of Plasmas*, vol. 25, no. 061201, 2018.

[95] J.-P. Boeuf and B. Chaudhury, “Rotating instability in low-temperature magnetized plasmas,” *Physical Review Letters*, vol. 111, no. 155005, 2013.

[96] J.-P. Boeuf, “Rotating structures in low temperature magnetized plasmas—insight from particle simulations,” *Frontiers in Physics*, vol. 2, no. 74, 2014.

[97] D. Escobar Anton and E. Ahedo, “Analysing the azimuthal spoke oscillation of Hall thrusters via numerical simulation,” in *50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, Cleveland, OH*, 2014.

[98] J. A. Carlsson, I. D. Kaganovich, A. V. Khrabrov, A. Smolyakov, D. Sydorenko, and Y. Raitses, “Multi-dimensional kinetic simulations of instabilities and transport in E×B devices,” in *2015 IEEE International Conference on Plasma Sciences (ICOPS), Antalya, Turkey*, IEEE, 2015.

[99] D. Escobar and E. Ahedo, “Global stability analysis of azimuthal oscillations in Hall thrusters,” *IEEE Transactions on Plasma Science*, vol. 43, pp. 149–157, 2014.

[100] D. Escobar and E. Ahedo, “Low frequency azimuthal stability of the ionization region of the Hall thruster discharge. II. Global analysis,” *Physics of Plasmas*, vol. 22, no. 102114, 2015.

[101] K. Matyash, R. Schneider, O. Kalentev, Y. Raitses, and N. Fisch, “Numerical investigations of a cylindrical Hall thruster,” *arXiv preprint arXiv:1208.4151*, 2012.

[102] F. Taccogna and P. Minelli, “Three-dimensional particle-in-cell model of Hall thruster: The discharge channel,” *Physics of Plasmas*, vol. 25, no. 061208, 2018.

[103] F. Taccogna, P. Minelli, M. Capitelli, and S. Longo, “Physics of Hall-effect thruster by particle model,” in *28th International Symposium on Rarefied Gas Dynamics, Zaragoza, Spain*, vol. 1501, pp. 1390–1399, American Institute of Physics, 2012.

[104] K. Matyash, R. Schneider, S. Mazouffre, S. Tsikata, and L. Grimaud, “Investigation of rotating spoke instabilities in a wallless Hall thruster. part II: Simulation,” in *35th International Electric Propulsion Conference, Atlanta, Georgia*, 2017.

[105] R. Kawashima, K. Hara, and K. Komurasaki, “Numerical analysis of azimuthal rotating spokes in a crossed-field discharge plasma,” *Plasma Sources Science and Technology*, vol. 27, no. 035010, 2018.

[106] J.-P. Boeuf and A. Smolyakov, “Physics and instabilities of low-temperature E×B plasmas for spacecraft propulsion and other applications,” *Physics of Plasmas*, vol. 30, no. 050901, 2023.
[107] O. Koshkarov, A. Smolyakov, Y. Raitses, and I. Kaganovich, "Self-organization, structures, and anomalous transport in turbulent partially magnetized plasmas with crossed electric and magnetic fields," *Physical Review Letters*, vol. 122, no. 185001, 2019.

[108] M. Tyushev, M. Papahn Zadeh, V. Sharma, M. Sengupta, Y. Raitses, J.-P. Boeuf, and A. Smolyakov, "Azimuthal structures and turbulent transport in Penning discharge," *Physics of Plasmas*, vol. 30, no. 033506, 2023.

[109] L. Xu, D. Eremin, and R. P. Brinkmann, "Direct evidence of gradient drift instability being the origin of a rotating spoke in a crossed field plasma," *Plasma Sources Science and Technology*, vol. 30, no. 075013, 2021.

[110] J.-P. Boeuf, "Spoke formation in low temperature E×B plasmas: Transition from gradient-drift instability to ionization wave," *Physics of Plasmas*, vol. 30, no. 022112, 2023.

[111] J. Boeuf, "Micro instabilities and rotating spokes in the near-anode region of partially magnetized plasmas," *Physics of Plasmas*, vol. 26, no. 072113, 2019.

[112] A. Mansour and K. Hara, "Full fluid moment modeling of rotating spokes in Penning-type configuration," *Plasma Sources Science and Technology*, vol. 31, no. 055012, 2022.

[113] V. Smirnov, S. Kislenko, A. Gavrikov, A. Klimkov, G. Liziakin, and V. Smirnov, "Particle-in-cell simulation of a cylindrical plasma mass separator based on Penning discharge with thermionic cathodes," *Plasma Sources Science and Technology*, vol. 32, no. 095007, 2023.

[114] L. Xu, H. Sun, D. Eremin, S. Ganta, I. Kaganovich, K. Bera, S. Rauf, and X. Wu, "Rotating spokes, potential hump and modulated ionization in radio frequency magnetron discharges," *Plasma Sources Science and Technology*, vol. 32, no. 105012, 2023.

[115] S. S. Ganta, K. Bera, S. Rauf, I. Kaganovich, A. Khrabrov, A. T. Powis, D. Sydorenko, and L. Xu, "Investigating instabilities in magnetized low-pressure capacitively coupled RF plasma using particle-in-cell (pic) simulations," *Physics of Plasmas*, vol. 31, no. 102107, 2024.

[116] A. S. Rokhmanenkov and S. E. Kuratov, "Numerical simulation of Penning gas discharge in 2D/3D setting," *Journal of Physics: Conference Series*, vol. 1250, p. 012036, jun 2019.

[117] M. Sengupta, A. Smolyakov, and Y. Raitses, "Restructuring of rotating spokes in response to changes in the radial electric field and the neutral pressure of a cylindrical magnetron plasma," *Journal of Applied Physics*, vol. 129, no. 223302, 2021.

[118] M. Sengupta and A. Smolyakov, "Mode transitions in nonlinear evolution of the electron drift instability in a 2D annular E×B system," *Physics of Plasmas*, vol. 27, no. 022309, 2020.

[119] M. Tyushev, M. Papahn Zadeh, N. Chopra, Y. Raitses, I. Romadanov, A. Likhanskii, G. Fubiani, L. Garrigues, R. Groenewald, and A. Smolyakov, "Mode transitions and spoke structures in E×B Penning discharge," *Physics of Plasmas*, vol. 32, no. 013511, 2025.

[120] J. Chen, A. T. Powis, I. D. Kaganovich, Z. Wang, and Y. Yu, "Three-dimensional helical-rotating plasma structures in beam-generated partially magnetized plasmas," *Physical Review Letters*, vol. 135, no. 045301, 2025.

[121] D. Poli, P. Fajardo, E. Ahedo, and G. Hagelaar, "Comparison of two fluid codes for modelling 2D instabilities in E×B plasmas," in *38th International Electric Propulsion Conference, Toulouse, France*, 2024.

[122] J.-P. Boeuf and M. Takahashi, "Rotating spokes, ionization instability, and electron vortices in partially magnetized E×B plasmas," *Physical Review Letters*, vol. 124, no. 185005, 2020.

[123] J.-P. Boeuf and M. Takahashi, "New insights into the physics of rotating spokes in partially magnetized E×B plasmas," *Physics of Plasmas*, vol. 27, no. 083520, 2020.
[124] W. Frias, A. I. Smolyakov, I. D. Kaganovich, and Y. Raitses, "Long wavelength gradient drift instability in Hall plasma devices. II. Applications," *Physics of Plasmas*, vol. 20, no. 052108, 2013.

[125] K. Hara, A. R. Mansour, and S. Tsikata, "Theory of gradient drift instabilities in low-temperature, partially magnetised plasmas," *Journal of Plasma Physics*, vol. 88, no. 905880408, 2022.

[126] J. Carlsson, A. Khrabrov, I. Kaganovich, T. Sommerer, and D. Keating, "Validation and benchmarking of two particle-in-cell codes for a glow discharge," *Plasma Sources Science and Technology*, vol. 26, no. 014003, 2016.

[127] Z. Donkó, A. Derzsi, I. Korolov, P. Hartmann, S. Brandt, J. Schulze, B. Berger, M. Koepke, B. Bruneau, E. Johnson, et al., "Experimental benchmark of kinetic simulations of capacitively coupled plasmas in molecular gases," *Plasma Physics and Controlled Fusion*, vol. 60, no. 014010, 2018.

[128] D. Eremin, E. Kemaneci, M. Matsukuma, T. Mussenbrock, and R. P. Brinkmann, "Modeling of very high frequency large-electrode capacitively coupled plasmas with a fully electromagnetic particle-in-cell code," *Plasma Sources Science and Technology*, vol. 32, no. 044007, 2023.

[129] D. Eremin, Y. Sharova, L. C. Heijmans, A. M. Yakunin, M. van de Kerkhof, R. P. Brinkmann, and E. H. Kemaneci, "Electromagnetic particle-in-cell modeling of an electron cyclotron resonance plasma discharge in hydrogen," *Plasma Sources Science and Technology*, vol. 34, no. 065001, 2024.

[130] R. Lucken, A. Tavant, A. Bourdon, M. Lieberman, and P. Chabert, "Saturation of the magnetic confinement in weakly ionized plasma," *Plasma Sources Science and Technology*, vol. 29, no. 065014, 2020.

[131] A. Domínguez-Vázquez, F. Taccogna, and E. Ahedo, "Particle modeling of radial electron dynamics in a controlled discharge of a Hall thruster," *Plasma Sources Science and Technology*, vol. 27, no. 064006, 2018.

[132] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al., "SciPy 1.0: Fundamental algorithms for scientific computing in python," *Nature methods*, vol. 17, pp. 261-272, 2020.

[133] L. Garrigues, G. Fubiani, and J. P. Boeuf, "Appropriate use of the particle-in-cell method in low temperature plasmas: Application to the simulation of negative ion extraction," *Journal of Applied Physics*, vol. 21, no. 213303, 2016.

[134] G. Fubiani, L. Garrigues, G. Hagelaar, N. Kohen, and J.-P. Boeuf, "Modeling of plasma transport and negative ion extraction in a magnetized radio-frequency plasma source," *New Journal of Physics*, vol. 19, no. 015002, 2017.

[135] J. Porto, P. Q. Elias, and A. Ciardi, "Anisotropic electron heating in an electron cyclotron resonance thruster with magnetic nozzle," *Physics of Plasmas*, vol. 30, no. 023506, 2023.

[136] A. Marín-Cebrián, E. Bello-Benítez, A. Domínguez-Vázquez, and E. Ahedo, "Non-maxwellian electron effects on the macroscopic response of a Hall thruster discharge from an axial-radial kinetic model," *Plasma Sources Science and Technology*, vol. 33, no. 025008, 2024.

[137] E. Bello-Benítez, A. Marín-Cebrián, and E. Ahedo, "Effect of injection conditions on the nonlinear behavior of the ECDI and related turbulent transport," *Physics of Plasmas*, vol. 32, no. 012307, 2025.

[138] D. Sydorenko, A. Khrabrov, W. Villaphana, A. Powis, and S. Ethier, "EDIPIC-2DWeb page," https://github.com/PrincetonUniversity/EDIPIC-2D, 2021.

[139] V. Croes, T. Lafleur, Z. Bonaventura, A. Bourdon, and P. Chabert, "2D particle-in-cell simulations of the electron drift instability and associated anomalous electron transport in Hall-effect thrusters," *Plasma Sources Science and Technology*, vol. 26, no. 034001, 2017.
[140] P. Parodi and F. Petronio, “Step-by-step verification of particle-in-cell Monte Carlo collision codes,” *Physics of Plasmas*, vol. 32, no. 013902, 2025.

[141] K. Hara and S. Tsikata, “Cross-field electron diffusion due to the coupling of drift-driven microinstabilities,” *Physical Review E*, vol. 102, no. 023202, 2020.

[142] K. Hara, T. Robertson, J. Kenney, and S. Rauf, “Effects of macroparticle weighting in axisymmetric particle-in-cell Monte Carlo collision simulations,” *Plasma Sources Science and Technology*, vol. 32, no. 015008, 2023.

[143] P. Parodi, S. Boccelli, F. Bariselli, and T. Magin, “Pantera: A PIC-MCC-DSMC software for the simulation of rarefied gases and plasmas,” *SoftwareX*, vol. 31, no. 102244, 2025.

[144] W. Villafana, B. Cuenot, and O. Vermorel, “3D particle-in-cell study of the electron drift instability in a Hall Thruster using unstructured grids,” *Physics of Plasmas*, vol. 30, no. 033503, 2023.

[145] G. Bogopolsky, *Exploration of methods for the numerical simulation of Hall thrusters*. PhD Thesis, Université de Toulouse, 2025.

[146] F. Cichocki, V. Sciortino, P. Minelli, and F. Taccogna, “PIC modeling challenges in diverse low temperature plasma scenarios,” in *26th Europhysics Conference on Atomic and Molecular Physics of Ionized Gases*, Brno, Czech Republic, 2024.

[147] F. Cichocki, F. Napoli, D. Iannarelli, and F. Taccogna, “Numerical study of rotating magnetic fields for enhanced plasma propulsion,” in *39th International Electric Propulsion Conference*, London, United Kingdom, 2025.

[148] F. Taccogna, M. Rutigliano, P. Minelli, F. Paganucci, S. Camarri, G. Becatti, and A. Micciani, “Particle-in-cell model of a compact diamond-based disk emitter plasma cathode,” in *39th International Electric Propulsion Conference*, London, United Kingdom, 2025.

[149] H. Timko, P. Crozier, M. Hopkins, K. Matyash, and R. Schneider, “Why perform code-to-code comparisons: A vacuum arc discharge simulation case study,” *Contributions to Plasma Physics*, vol. 52, pp. 295-308, 2012.

[150] F. Faraji, M. Reza, and A. Knoll, “Verification of the generalized reduced-order particle-in-cell scheme in a radial–azimuthal E×B plasma configuration,” *AIP Advances*, vol. 13, no. 025315, 2023.

[151] O. Schenk and K. Gärtner, “Solving unsymmetric sparse systems of linear equations with PARDISO,” *Future Generation Computer Systems*, vol. 20, pp. 475–487, 2004.

[152] NAG FORTRAN Library, *Numerical Algorithms Group Oxford*. www.nag.com/numeric/fl/FLdescription.asp.

[153] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, “Efficient management of parallelism in object oriented numerical software libraries,” in *Modern Software Tools in Scientific Computing* (E. Arge, A. M. Bruaset, and H. P. Langtangen, eds.), pp. 163–202, Birkhäuser Press, 1997.

[154] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Ejkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, and H. Zhang, “PETSc users manual,” Tech. Rep. ANL-95/11 - Revision 3.6, Argonne National Laboratory, 2015.

[155] “Animation of contour plots of ion density of the Penning discharge benchmark..” https://www.youtube.com/watch?v=XWr-3vAUuvA&feature=youtu.be.

[156] D. Forslund, R. Morse, and C. Nielson, “Electron cyclotron drift instability,” *Physical Review Letters*, vol. 25, no. 1266, 1970.
[157] S. Janhunen, A. Smolyakov, D. Sydorenko, M. Jimenez, I. Kaganovich, and Y. Raitses, "Evolution of the electron cyclotron drift instability in two-dimensions," *Physics of Plasmas*, vol. 25, no. 082308, 2018.

[158] J. Reichmann, W. *Use and Abuse of Statistics*. New York: Oxford University Press, 1962. Nontechnical exposition.

[159] S. B. Pope, *Turbulent Flows*. Cambridge, United Kingdom: Cambridge University Press, 1st ed., 2000. Comprehensive text on turbulence theory and modeling, covering statistical methods, Reynolds-averaged and LES approaches.

[160] K. Germaschewski, W. Fox, S. Abbott, N. Ahmadi, K. Maynard, L. Wang, H. Ruhl, and A. Bhat-tacharjee, "The plasma simulation code: A modern particle-in-cell code with patch-based load-balancing," *Journal of Computational Physics*, vol. 318, pp. 305-326, 2016.

[161] K. G. Miller, R. P. Lee, A. Tableman, A. Helm, R. A. Fonseca, V. K. Decyk, and W. B. Mori, "Dynamic load balancing with enhanced shared-memory parallelism for particle-in-cell codes," *Computer Physics Communications*, vol. 259, no. 107633, 2021.

[162] H. Nakashima, Y. Miyake, H. Usui, and Y. Omura, "OhHelp: a scalable domain-decomposing dynamic load balancing for particle-in-cell simulations," in *Proceedings of the 23rd international conference on Supercomputing, Yorktown Heights, NY*, pp. 90-99, 2009.

[163] R. Jambunathan and D. A. Levin, "Chaos: An octree-based PIC-DSMC code for modeling of electron kinetic properties in a plasma plume using mpi-cuda parallelization," *Journal of computational physics*, vol. 373, pp. 571-604, 2018.

[164] F. Taccogna, F. Cichocki, D. Eremin, G. Fubiani, and L. Garrigues, "Plasma propulsion modeling with particle-based algorithms," *Journal of Applied Physics*, vol. 134, no. 150901, 2023.

[165] Y. Saad, *Iterative Methods for Sparse Linear Systems*. Philadelphia, PA: Society for Industrial and Applied Mathematics, 2nd ed., 2003.

[166] H. R. Lewis, "Energy-conserving numerical approximations for Vlasov plasmas," *Journal of Computational Physics*, vol. 6, pp. 136-141, 1970.

[167] D. C. Barnes and L. Chacón, "Finite spatial-grid effects in energy-conserving particle-in-cell algorithms," *Computer Physics Communications*, vol. 258, no. 107560, 2021.

[168] A. T. Powis and I. D. Kaganovich, "Accuracy of the explicit energy-conserving particle-in-cell method for under-resolved simulations of capacitively coupled plasma discharges," *Physics of Plasmas*, vol. 31, no. 023901, 2024.

[169] G. Lapenta, "Exactly energy conserving semi-implicit particle in cell formulation," *Journal of Computational Physics*, vol. 334, pp. 349-366, 2017.

[170] G. Chen, L. Chacón, and D. C. Barnes, "An energy-and charge-conserving, implicit, electrostatic particle-in-cell algorithm," *Journal of Computational Physics*, vol. 230, pp. 7018-7036, 2011.

[171] D. Eremin, "An energy-and charge-conserving electrostatic implicit particle-in-cell algorithm for simulations of collisional bounded plasmas," *Journal of Computational Physics*, vol. 452, no. 110934, 2022.

[172] A. B. Langdon, B. I. Cohen, and A. Friedman, "Direct implicit large time-step particle simulation of plasmas," *Journal of Computational Physics*, vol. 51, pp. 107-138, 1983.

[173] M. R. Gibbons and D. W. Hewett, "The Darwin direct implicit particle-in-cell (DADIPIC) method for simulation of low frequency plasma phenomena," *Journal of Computational Physics*, vol. 120, pp. 231-247, 1995.

[174] H. Sun, S. Banerjee, S. Sharma, A. T. Powis, A. V. Khrabrov, D. Sydorenko, J. Chen, and I. D. Kaganovich, "Direct implicit and explicit energy-conserving particle-in-cell methods for modeling of capacitively coupled plasma devices," *Physics of Plasmas*, vol. 30, no. 103509, 2023.
[175] F. Deluzet, G. Fubiani, L. Garrigues, C. Guillet, and J. Narski, “Sparse grid reconstructions for particle-in-cell methods,” *ESAIM: Mathematical Modelling and Numerical Analysis*, vol. 56, pp. 1809–1841, 2022.

[176] L. Garrigues, B. Tezenas du Montcel, G. Fubiani, F. Bertomeu, F. Deluzet, and J. Narski, “Application of sparse grid combination techniques to low temperature plasmas particle-in-cell simulations. I. Capacitively coupled radio frequency discharges,” *Journal of Applied Physics*, vol. 129, no. 153303, 2021.

[177] L. Garrigues, B. Tezenas du Montcel, G. Fubiani, and B. Reman, “Application of sparse grid combination techniques to low temperature plasmas particle-in-cell simulations. II. Electron drift instability in a Hall thruster,” *Journal of Applied Physics*, vol. 129, no. 153304, 2021.

[178] L. Garrigues, M. Chung-To-Sang, G. Fubiani, C. Guillet, F. Deluzet, and J. Narski, “Acceleration of particle-in-cell simulations using sparse grid algorithms. I. Application to dual frequency capacitive discharges,” *Physics of Plasmas*, vol. 31, no. 073907, 2024.

[179] L. Garrigues, M. Chung-To-Sang, G. Fubiani, C. Guillet, F. Deluzet, and J. Narski, “Acceleration of particle-in-cell simulations using sparse grid algorithms. II. Application to partially magnetized low temperature plasmas,” *Physics of Plasmas*, vol. 31, no. 073908, 2024.

[180] F. Faraji, M. Reza, and A. Knoll, “Enhancing one-dimensional particle-in-cell simulations to self-consistently resolve instability-induced electron transport in Hall thrusters,” *Journal of Applied Physics*, vol. 131, no. 193302, 2022.

[181] M. Reza, F. Faraji, and A. Knoll, “Concept of the generalized reduced-order particle-in-cell scheme and verification in an axial-azimuthal Hall thruster configuration,” *Journal of Physics D: Applied Physics*, vol. 56, no. 175201, 2023.

[182] M. Reza, F. Faraji, and A. Knoll, “Latest verifications of the reduced-order particle-in-cell scheme: Penning discharge and axial-radial Hall thruster case,” in *AIAA SCITECH 2024 Forum*, 2024.