# Benchmark for two-dimensional large scale coherent structures in partially magnetized E x 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, et al.

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

## 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 E x B plasmas-community collaboration & lessons learned. Plasma Sources Science and Technology, 2026, 35 (2), pp.025002. 10.1088/1361-6595/ae3985. hal-05502126

HAL Id: hal-05502126

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

Submitted on 10 Feb 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - 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³, 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¹⁰, Denis Eremin¹¹, Pablo Fajardo², Farbod Faraji¹², Gwenael Fubiani⁵, Laurent Garrigues⁵, Kentaro Hara⁸, Peter Hartmann⁹, Matthew Hopkins⁴, Igor D. Kaganovich¹, Aaron Knoll¹², Giovanni Lapenta¹³, Thierry Magin¹³, Alberto Marín Cebrián², Mario Merino Martínez², Pierpaolo Minelli⁷, Pietro Parodi¹³, Federico Petronio³, Maryam Reza¹², Andrei I. Smolyakov¹⁴, Dmytro Sydorenko¹,¹⁴,¹⁵, Francesco Taccogna⁷, Miles M. Turner¹⁶, Olivier Vermorel⁶, Gabriel Vigot⁶, Willca Villafana¹, Liang Xu¹¹, and Mina P. Zadeh¹⁴

¹Princeton Plasma Physics Laboratory, Princeton, New Jersey 08640, 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

⁵LAPLACE, Université de Toulouse, CNRS, INPT, UPS, 118 Route de Narbonne, 31062 Toulouse, France

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

⁷Institute for Plasma Science and Technology (ISTP), CNR, Bari, 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

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

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

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

¹⁶Dublin City University, Dublin, Ireland

June 5, 2025

**Abstract**

Low-temperature plasmas are essential to both fundamental scientific research and critical in-
dustrial applications. As in many areas of science, numerical simulations have become a vital
tool for uncovering new physical phenomena and guiding technological development. Code bench-
marking remains crucial for verifying implementations and evaluating performance. This work
continues the *Landmark Benchmark* initiative, a series specifically designed to support the veri-
fication of kinetic low-temperature plasma codes. In this study, sixteen simulation codes from a
collaborative community of fifteen international institutions modeled a partially magnetized E×B Penning discharge. The emergence of large scale coherent structures, or 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 required considerable effort, and we share several insights and best practices for conducting effective benchmarking campaigns.

# 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 of several electron-volts, are challenging to run and diagnose non-intrusively. Furthermore, the high-dimensional, non-linear, and multiscale nature of governing equations often prevents an analytical analysis. 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. This progress has been aided by the continuous improvement of computing performance, particularly with the widespread adoption of computational accelerators [8], as well as the integration of traditional modeling techniques with machine learning [9]. Nonetheless, computer science and algorithm design for LTPs is still a cutting edge area of research (see Section VII in Ref. [10]).

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 a choice of 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 model kinetic plasmas, including LTPs [11–13]. 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. Collisions between charged species and neutrals are implemented via Monte-Carlo algorithms [14]. PIC is also highly parallelizable, with particles evolved in a way that is nearly embarrassingly parallel, and field solvers can take advantage of distributed linear algebra packages [15–17]. Furthermore, the particle and collision routines are amenable to GPU acceleration, making the algorithm even more appealing for use on modern heterogeneous supercomputing architectures [8, 10, 18].

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 [19] 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 [20] as well as the two Landmark benchmarks [21, 22]. The aim of this work is to extend the Landmark series with a new benchmark configuration and through the involvement of a larger number of institutions.

Two previous benchmarks from the Landmark series, specifically 2a [21] and 2b [22] 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 confine electrons, whereas the heavier ions move freely under the influence of electric fields. Without collisions, and despite being confined by the applied field, electrons are eventually transported across the field lines from the cathode to the anode. The underlying mechanism for this “anomalous” transport is
still an active area of investigation, with numerous mechanisms being proposed. These include near-wall conductivity [23, 24] and a wide range of plasma waves and instabilities [25–33], 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 plasma “spoke” [34], a low-frequency, large-scale coherent structure observed in numerous E×B experiments [35–86] and simulations [87–117].

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) [57, 58, 87, 103–105, 108, 118]. 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 [119]) 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 [101, 103, 113]. 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 [34, 102, 104]. One theory also proposes that observations of the spoke are projections of a 3D helical structures [114]. 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 [34, 37, 39, 87, 88, 101, 102, 105], whether due to transport being enhanced by the wave itself, or the presence of a majority of the plasma current within the mode. 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.

<figure><img src="image_3.png" /><figcaption>Figure 1: Illustration of a typical Penning discharge experiment with a beamline source of energetic electrons. 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 sets up a radial ambipolar electric field (blue arrows) and corresponding azimuthal E×B drift (green arrow). The resulting azimuthal streaming of electrons against slowly moving ions leads to formation of the rotating spoke.</figcaption></figure>

The Penning discharge is an E×B device which exhibits many of the same phenomena observed in Hall thrusters, albeit with a simplified geometry [57–59, 64, 66]. This includes the presence of anomalous transport and rotating spokes. The Penning discharge 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 be the magnetic field. A radial ambipolar field emerges which drives E×B processes and associated instabilities, including the spoke. The axial dynamics take place at a relatively short timescale when compared to the radial and azimuthal motion of the plasma, and it is theoretically reasonable to decouple this axis from the radial-azimuthal plane. The numerical benchmark presented in this work will model a two-dimensional cross-section 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. 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 [61,66,67], 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. [120–122]). 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 presents the results of the benchmark and provides some important lessons learned through the lengthy process of reaching code agreement. Finally Section 4 provides 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 critical. 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 [87], typical Argon and Xenon gases would have led to excessive runtimes, dictating the choice of Helium-4. The simulation domain is a two-dimensional square {$x, y$}, centered about the origin, with edge lengths $L_x = L_y = 2.5$ cm (see Fig. 2). The square boundary of the simulation domain is grounded with $\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$ Gauss. 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_4.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).</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 each cell is perfectly square ($\Delta x = \Delta y$). The cell size is small enough 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 current of electrons $I_e$ and ions $I_i$ are injected into the center of the domain, each with 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}} \right) \quad (1)$$

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

The unit-depth of the simulation (i.e. in the unresolved z-direction) is 1 m, therefore all currents are listed in units of A/m. The injection profile is a cylindrical Heaviside function with form,

$$I_{\alpha}(x, y) = \begin{cases} I_{\alpha,0}, & x^2 + y^2 \le R_{inj} \\ 0, & \text{otherwise} \end{cases} \quad (2)$$

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<sub>i</sub></td><td>7291.712</td><td>m<sub>e</sub></td></tr><tr><td>Time step size</td><td>Δt</td><td>40</td><td>ps</td></tr><tr><td>Number of time steps</td><td>N<sub>t</sub></td><td>12,500,000</td><td>-</td></tr><tr><td>Total simulation time</td><td>t<sub>sim</sub></td><td>500</td><td>μs</td></tr><tr><td>Steps between diagnostic print</td><td>N<sub>print</sub></td><td>12,500</td><td>-</td></tr><tr><td>Time between diagnostic print</td><td>Δt<sub>print</sub></td><td>0.5</td><td>μs</td></tr><tr><td>Domain length in x/y-direction</td><td>L<sub>x</sub>, L<sub>y</sub></td><td>5.0 × 10<sup>-2</sup></td><td>m</td></tr><tr><td>Cells in x/y-direction</td><td>N<sub>x</sub>, N<sub>y</sub></td><td>256</td><td>-</td></tr><tr><td>Cell length in x/y-direction</td><td>Δx, Δ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>Gauss</td></tr><tr><td>Injection radius</td><td>R<sub>inj</sub></td><td>5.0 × 10<sup>-3</sup></td><td>m</td></tr><tr><td>Electron injection current</td><td>I<sub>e,0</sub></td><td>-20.0 × 10<sup>-3</sup></td><td>A/m</td></tr><tr><td>Electron injection temperature</td><td>T<sub>e,0</sub></td><td>15.0</td><td>eV</td></tr><tr><td>Ion injection current</td><td>I<sub>i,0</sub></td><td>8.0 × 10<sup>-3</sup></td><td>A/m</td></tr><tr><td>Ion injection temperature</td><td>T<sub>i,0</sub></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<sub>Δt,e</sub></td><td>49.932</td><td>-</td></tr><tr><td>Ion macro-particles injected per time step</td><td>N<sub>Δt,i</sub></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 cross-sections of the above listed quantities are used as the primary comparison between codes. Cross-sections 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 corresponds to extracting a single slice of grid data at the y-node location $j = 128$. 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 in $x \in [0, L_x/2]$, providing a single “radial” 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. We denote the time average of a measured quantity $q$ as $\langle q \rangle_t$. 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 in direction $\beta \in \{x, y, z\}$ at cell node $\{i, j\}$ (where $i$ here is not to be confused with ions) is computed as,

$$[nV_{\beta}^{m}]_{i,j} = W \sum_{p=0}^{P} V_{p,\beta}^{m} S\left(\frac{X_p - x_i}{\Delta x}\right) S\left(\frac{Y_p - y_j}{\Delta y}\right) \quad (3)$$

Where $P$ is the total number of macro-particles, $X_p$ and $Y_p$ are the position coordinates of particle $p$ and $V_{p,\beta}$ is the velocity of particle $p$ in direction $\beta$. $S(x)$ is the standard linear shape function,
$$S(x) = \begin{cases} 1+x & -1 < x < 0 \\ 1-x & 0 \le x < 1. \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{(nV_{\beta}^2)}{n} - \left( \frac{(nV_{\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, and for that reason we consider two definitions in this work.

**Temperature Type I:** The standard approach, 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{(nV_{\beta}^2)}{n} - \left( \frac{(nV_{\beta})}{n} \right)^2 \right\rangle_t. \qquad (6)$$

**Temperature II:** A non-standard approach, choosing to average each velocity moment in time prior to calculating the temperature,

$$\langle T_{\beta} \rangle_t^{II} = m \left[ \frac{\langle nV_{\beta}^2 \rangle_t}{\langle n \rangle_t} - \left( \frac{\langle nV_{\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^X = \frac{1}{3} \sum_{\beta} \langle T_{\beta} \rangle_t^X$. With $X \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 = \frac{1}{2}L_x, y = 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 \ \mu s \le t < 500 \ \mu s$, and peaks are identified via a peak-finding algorithm [123] (care must be taken to ensure that false peaks are ignored). 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_5.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 [123].</figcaption></figure>

## 2.4 Code details

A record 16 codes were involved in this iteration of the *Landmark* benchmarking effort, each of them testing independently developed PIC codes. These included the Princeton Plasma Physics Laboratory - PPPL (LTP-PIC code [18]), Laplace laboratory at the University of Toulouse (EP-PIC2D [124] and LePIC2D codes []), Ruhr University Bochum - RUB (RUBPIC code []), Wigner Research Center for Physics (MZ/X-PIC code []), Dublin City University - DCU (... code []), the French National Aerospace Research Centre - Onera (RHEI code []), Carlos III University of Madrid - UC3M (PICASO code []), University of Saskatchewan - USask (the Open-Source EDIPIC-2D code [125] - also developed with PPPL and the University of Alberta), École Polytechnique in Paris - LPP (LPPic code []), Stanford University (XPIC2D code []), The von Karman Institute for Fluid Dynamics - VKI (PANTERA code []), the European Centre for Research and Advanced Training in Scientific Computing - CERFACS (AVIP-PIC code []), the Italian National Agency for New Technologies, Energy and Sustainable Economic Development in Bari - CNR-Bari (PICCOLO code []), Sandia National Laboratories (Aleph code []), and the Imperial College London - ICL (IPPL-2D code []).

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 implicit time advancement (the RUBPIC code) and unstructured meshes (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 include Hypre [16], PARDISO [126], NAG [127], PETSc [15, 128, 129] and Trilinos [17]. 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>Runtime</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>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>RUB</td><td>RUBPIC</td><td>C/C++</td><td>NAG geometric multigrid</td><td>CUDA</td><td>Particle</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 GTX-1080Ti GPU + 1x Intel i7 CPU core</td><td>10 days</td></tr><tr><td>DCU</td><td>DCU</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</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><tr><td>USask</td><td>EDIPIC-2D</td><td>Fortran</td><td>PETSc</td><td>MPI</td><td>Domain</td><td>256x Intel E5-2683 v4 2.1 GHz CPU cores</td><td>42 days</td></tr></tbody></table>
<table><tr><td>LPP</td><td>LPPic</td><td>Fortran</td><td>Hypre</td><td>MPI</td><td>Domain</td><td>128x AMD Rome<br>2.6 GHz<br>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<br>7742 2.25 GHz<br>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<br>Skylake Xeon<br>Platinum 8174<br>CPU cores</td><td>72 hours</td></tr><tr><td>Cerfacs</td><td>AVIP-PIC</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>CNR-Bari</td><td>PICCOLO</td><td>Fortran</td><td>PETSc +<br>hypre</td><td>MPI+OpenMP</td><td>Hybrid</td><td>42x AMD<br>EPYC 7402<br>CPU cores</td><td>30 days</td></tr><tr><td>Sandia</td><td>Aleph</td><td>C/C++</td><td>Trilinos</td><td>MPI</td><td>Hybrid</td><td>640x Intel<br>Sandy Bridge<br>2.5 GHz<br>CPU cores</td><td>14 days</td></tr><tr><td>ICL</td><td>IPPL-2D</td><td>Julia</td><td>Julia direct<br>solver</td><td>MPI</td><td>Domain</td><td>-</td><td>-</td></tr></table>
# 3 Results & discussion

Before comparing simulations across participating groups, we review a few of the key results for a single simulation. Figure 4 shows contour plots of ion density over a single period of the rotating spoke, with an animation available on YouTube. Using the diagnostic described in Sec. 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.4 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 Sec. 3.2 we will discuss how this can cause issues when comparing measurements, such as temperature, between codes.

<figure><img src="image_6.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 Sec. 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 is a saturation of the collisionless Simon-Hoh instability (a type of gradient drift instability). Under this assumption, Eq. 4 in Ref. [87] 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 also demonstrates 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 [31, 113, 130, 131]. 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 collected 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 mean frequency of 43.4 kHz and an uncertainty, defined by the standard deviation, of 3.0 kHz. 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>Total time</th><th>Spoke frequency</th><th>Frequency discrepancy</th></tr></thead><tbody><tr><td>PPPL</td><td>4.0 ms</td><td>43.4 ± 3.0 kHz</td><td>-</td></tr><tr><td>Laplace<sup>I</sup></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>1.6 ms</td><td>42.6 kHz</td><td>0.8 kHz</td></tr></tbody></table>
Figure 5 compares time averaged “radial” discharge cross-sections 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 value 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.

<figure><img src="image_7.png" /><figcaption>Figure 5: Comparison of time-averaged cross-sections for long-time simulations, (a) ion-density, (b) plasma potential, (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 discharge cross-sections for (a) ion-density, (b) electric potential, and (c) electron temperature (Type I). Note that the results from Sec. 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_8.png" /><figcaption>Figure 6: Comparison of time-averaged cross-sections between all codes, (a) ion-density, (b) plasma potential, (c) electron temperature (Type I).</figcaption></figure>

Understanding the cause of this discrepancy in electron temperature across simulations resulted in significant further analysis and 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 large 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 are obviously in Section 3.1.

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 measurements within a given cell which can easily skew the time averaged results. Furthermore, it is unclear whether all codes take into consideration Bessel's correction, which becomes important in poorly sampled cells.

Further issues considered 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. Sensitivity tests with
the PPPL and 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 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
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 correct definition of time-averaged temperature, it is a
valid way to compare numerical methods if done consistently. Figure 7 shows plots of Temperature
II computed via Eq. 7, demonstrating agreement between all previously outlier codes within the
uncertainty bounds.

<figure><img src="image_9.png" /><figcaption>Figure 7: Comparison of time-averaged cross-sections of electron temperature (Type II) for select codes. NOTE: Cerfacs.</figcaption></figure>

Finally, Fig. 8 compares the spoke rotation frequency measured by all groups, showing agreement
with the measured frequency of 43.4 ± 3.0 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 all codes across a
variety of measurements.
<figure><img src="image_10.png" /><figcaption>Figure 8: Comparison of the spoke rotation frequency, showing agreement between all codes. The reference frequency of 43.4 kHz is produced from long time simulations with the PPPL code with uncertainty bounds of 3.0 kHz defined as the standard deviation and 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.

**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. Choice of pseudo-random number generator was also a question raised at code author meetings, although it was ultimately demonstrated this had little influence on the results. Nonetheless, we emphasize that preparing a suitable benchmark description is paramount to success in any such effort.

**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 prohibitive 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 low-temperature plasma community, where the phenomena under investigation are often highly multi-scale, especially when collective effects and if 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 [132] 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 important physics, but not longer, and multi-scale and chaotic behavior should be avoided unless it is the express purpose of the benchmark to test the emergence of this physics.
<figure><img src="image_11.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>

**Lesson 3: Avoid systems which have regions of poor statistic resolution.** This lesson pertains particularly to Monte-Carlo methods, such as the Particle-in-Cell algorithm. The discussion in Sec. 3.2 highlights that poor statistical sampling within the vacuum region 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 [21] and 2b [22] benchmark configurations provide a more appropriate setup which avoids these pitfalls.

**Lesson 4: More data is better.** Even with a well prepared benchmark description it is impossible for the researcher to predict all possible ways in which disagreement may arise between codes. To prepare for this eventually 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 collect uncertainty bounds for each simulation, which can further simplify comparison.

Despite the significant time and effort that went into resolving issues with this benchmark, the authors want to emphasize that the process is still worthwhile, especially when such a large number of groups are involved. Coupled with validation and simpler forms of verification, benchmarking is critical to ensuring that a complex code can accurately reproduce results which are relevant to the broader community. Ultimately, agreement between all codes has been demonstrated, and we hope that the above lessons can serve as a cautionary tale for future efforts to benchmark plasma simulation codes, especially if they lead to saved time and effort.
## 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 which are worth discussing.

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, particle decomposition where the list of particles in memory is split (usually evenly) between MPI tasks which each own 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 an `MPI_Allreduce` call at each time step to share field information. Therefore most codes implement a domain decomposition approach, whereby each task owns a certain area (or volume in 3D) of the total grid domain and manages particles which are only present in that area. This avoids the need for an expensive 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.

Despite this trend, the value of domain decompensation cannot be understated. For large simulations, particle decomposition becomes increasingly expensive and eventually intractable due to 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 [133–135] where the domain owned by each process is adjusted on the fly based on the spatial density of the simulation particles.

A second observation from Table 2, and alluded to often above, is simply the enormous computational time required for each penning benchmark 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 the low-frequency spoke instability. Although in the previous section we recommended avoiding such long time scale and chaotic phenomena as a future benchmark, 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 good high-performance computing practices are essential to reduce the compute time per simulation step, and therefore total simulation time. Details on how to achieve these have been documented extensively [8, 10, 18, 136]. However most broadly, this comes down to understanding how to efficiently map the fundamental PIC algorithms on to relevant (and constantly changing) computational architecture for maximum performance. With a trend towards increasingly heterogeneity in compute systems, this requirement will only continue to grow in importance over the coming years.

Code developers could also consider alternative algorithms which can step over the usually strict resolution requirements of PIC, these include an energy-conserving approach for larger cell sizes [137–139], and implicit [140–142] (often also energy-conserving) or semi-implicit [143–145] techniques for time-stepping. Advanced approaches including sparse-grid PIC [146–148] and reduced-order PIC [149, 150] which have also proven accurate at modeling $E \times B$ discharges at reduced computational cost. Of course, for accurate simulations 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.
# 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 phenomena. 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. Sixteen codes, a record for this community, were put to the test on this benchmark, with agreement found for 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 authors hope that describing these lessons will help this, and other simulation communities, design improved benchmarks for their respective software. The effort also highlighted important computational considerations for particle-in-cell code design, including flexibility of parallelization strategy and adherence to high-performance computing best practices.

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

# 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. 37, p. 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. 3, 2020.

[3] A. Anders, “Tutorial: Reactive high power impulse magnetron sputtering (r-hipims),” *Journal of Applied Physics*, vol. 121, no. 17, 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. 9, 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, no. 8, pp. 2312–2319, 1993.

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

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

[8] 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, p. 107913, 2021.

[9] 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. 10, p. 102002, 2025.
[10] 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. 12, 2020.

[11] R. W. Hockney and J. W. Eastwood, *Computer simulation using particles*. CRC Press, 2021.

[12] C. K. Birdsall and A. B. Langdon, *Plasma physics via computer simulation*. CRC press, 2018.

[13] T. Tajima, *Computational plasma physics: with applications to fusion and astrophysics*. CRC press, 2018.

[14] 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.

[15] 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.

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

[17] T. Trilinos Project Team, *The Trilinos Project Website*.

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

[19] 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, p. 110751, 2022.

[20] 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. 1, 2013.

[21] 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. 10, p. 105010, 2019.

[22] 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. 7, p. 075002, 2021.

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

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

[25] E. Y. Choueiri, "Plasma oscillations in hall thrusters," *Physics of Plasmas*, vol. 8, no. 4, pp. 1411-1426, 2001.

[26] 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, no. 4, pp. 582-591, 2001.

[27] T. Lafleur, S. Baalrud, and P. Chabert, "Theory for the anomalous electron transport in hall effect thrusters. i. insights from particle-in-cell simulations," *Physics of Plasmas*, vol. 23, no. 5, 2016.
[28] 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. 5, 2016.

[29] 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. 1, p. 015003, 2017.

[30] 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. 16, 2006.

[31] 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. 6, 2020.

[32] 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. 6, p. 065204, 2023.

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

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

[35] 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, p. 928, 1974.

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

[37] C. L. Ellison, Y. Raitses, and N. J. Fisch, “Direct measurement of spoke-induced, cross-field electron transport in a cylindrical hall thruster,” 2011.

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

[39] 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. 1, 2012.

[40] 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. 5, 2012.

[41] 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*, p. 4179, 2012.

[42] 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, no. 10, pp. 2656–2657, 2014.

[43] 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, no. 1, pp. 72–85, 2014.

[44] M. Cappelli, Y. Raitses, and I. Kaganovich, “Coherent structures in plasmas relevant to electric propulsion,” tech. rep., Technical Report, 2016.

[45] N. Hyatt, “Hall effect thruster characterization through potential, magnetic, and optical measurements,” in *55th AIAA Aerospace Sciences Meeting*, p. 0169, 2017.
[46] 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. 5, p. 054002, 2019.

[47] 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. 5, 2012.

[48] 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. 8, p. 084007, 2013.

[49] 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. 4, p. 045005, 2015.

[50] 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. 16, 2015.

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

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

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

[54] 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. 9, p. 095204, 2018.

[55] 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. 11, 2012.

[56] 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. 8, p. 084005, 2013.

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

[58] 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, no. 6, pp. 1681-1694, 1993.

[59] 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*, 2015.

[60] 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. 6, 2018.

[61] 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. 6, p. 065004, 2004.

[62] 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.
[63] 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. 2, 2009.

[64] 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. 22, p. 225006, 2011.

[65] A. Escarguel, “Optical diagnostics of a low frequency instability rotating around a magnetized plasma column,” *The European Physical Journal D*, vol. 56, no. 2, pp. 209–214, 2010.

[66] 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. 2, p. 025011, 2021.

[67] 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. 5, 2019.

[68] 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. 7, p. 07LT01, 2023.

[69] 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. 5, p. L053202, 2022.

[70] 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. 10, 2023.

[71] R. C. Przybocki and M. A. Cappelli, “Comparison of spoke dynamics in an obstructed direct-current magnetron operating on a range of noble gases,” *Plasma Sources Science and Technology*, 2025.

[72] 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, p. 130242, 2025.

[73] 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. 12, 2018.

[74] 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. 1, p. 01LT01, 2023.

[75] 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, no. 10, pp. 1238–1249, 2024.

[76] 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.

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

[78] 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. 8, p. 085013, 2022.
[79] 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. 6, p. 065010, 2015.

[80] 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, no. 4, pp. 903-917, 2016.

[81] 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. 11, 2022.

[82] 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. 25, 2014.

[83] 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. 10, p. 102003, 2014.

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

[85] A. Anders, "Self-organization and self-limitation in high power impulse magnetron sputtering," *Applied Physics Letters*, vol. 100, no. 22, 2012.

[86] 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. 3, p. 035001, 2016.

[87] 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. 7, 2018.

[88] 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. 6, 2018.

[89] J.-P. Boeuf and B. Chaudhury, "Rotating instability in low-temperature magnetized plasmas," *Physical review letters*, vol. 111, no. 15, p. 155005, 2013.

[90] J.-P. Boeuf, "Rotating structures in low temperature magnetized plasmas—insight from particle simulations," *Frontiers in Physics*, vol. 2, p. 74, 2014.

[91] 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*, p. 3512, 2014.

[92] 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 exb devices," in *2015 IEEE International Conference on Plasma Sciences (ICOPS)*, pp. 1–1, IEEE, 2015.

[93] D. Escobar and E. Ahedo, "Global stability analysis of azimuthal oscillations in hall thrusters," *IEEE Transactions on Plasma Science*, vol. 43, no. 1, pp. 149–157, 2014.

[94] 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. 10, 2015.

[95] K. Matyash, R. Schneider, O. Kalentev, Y. Raitses, and N. Fisch, "Numerical investigations of a cylindrical hall thruster," *arXiv preprint arXiv:1208.4151*, 2012.

[96] F. Taccogna and P. Minelli, "Three-dimensional fully kinetic particle-in-cell model of hall-effect thruster," in *International Electric Propulsion Conference*, no. 088, 2011.
[97] F. Taccogna, P. Minelli, M. Capitelli, and S. Longo, "Physics of hall-effect thruster by particle model," in *AIP conference proceedings*, vol. 1501, pp. 1390-1399, American Institute of Physics, 2012.

[98] 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*, pp. IEPC-paper, 2017.

[99] 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. 3, p. 035010, 2018.

[100] 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. 5, 2023.

[101] 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. 18, p. 185001, 2019.

[102] 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. 3, 2023.

[103] 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. 7, p. 075013, 2021.

[104] 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. 2, 2023.

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

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

[107] 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. 9, p. 095007, 2023.

[108] 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. 10, p. 105012, 2023.

[109] 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. 10, 2024.

[110] A. Rokhmanenkov and S. Kuratov, "Numerical simulation of penning gas discharge in 2d/3d setting," in *Journal of Physics: Conference Series*, vol. 1250, p. 012036, IOP Publishing, 2019.

[111] 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. 22, 2021.

[112] 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. 2, 2020.

[113] 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. 1, 2025.
[114] J. Chen, A. T. Powis, I. D. Kaganovich, and Z. Wang, "Three-dimensional helical-rotating plasma structures in beam-generated partially magnetized e × b plasmas," *arXiv preprint arXiv:2501.01161*, 2025.

[115] D. Poli, P. Fajardo, E. Ahedo, and G. Hagelaar, "Comparison of two fluid codes for modelling 2d instabilities in exb plasmas," 2024.

[116] 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. 18, p. 185005, 2020.

[117] 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. 8, 2020.

[118] 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. 5, 2013.

[119] 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.

[120] 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. 1, p. 014003, 2016.

[121] 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. 1, p. 014010, 2018.

[122] 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. 4, p. 044007, 2023.

[123] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, "SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python," *Nature Methods*, vol. 17, pp. 261-272, 2020.

[124] 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. 120, no. 21, p. 213303, 2016.

[125] D. Sydorenko, *Particle-in-cell simulations of electron dynamics in low pressure discharges with magnetic fields*. PhD thesis, University of Saskatchewan, 2006.

[126] O. Schenk and K. Gärtner, "Solving unsymmetric sparse systems of linear equations with pardiso," *Future Generation Computer Systems*, vol. 20, no. 3, pp. 475-487, 2004.

[127] NAG FORTRAN Library, *Numerical Algorithms Group Oxford*.

[128] 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.

[129] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Ejjkhout, 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.
[130] D. Forslund, R. Morse, and C. Nielson, “Electron cyclotron drift instability,” *Physical Review Letters*, vol. 25, no. 18, p. 1266, 1970.

[131] 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. 8, 2018.

[132] S. B. Pope, “Turbulent flows,” *Measurement Science and Technology*, vol. 12, no. 11, pp. 2020–2021, 2001.

[133] 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.

[134] 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, p. 107633, 2021.

[135] 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*, pp. 90–99, 2009.

[136] 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. 15, 2023.

[137] H. R. Lewis, “Energy-conserving numerical approximations for vlasov plasmas,” *Journal of Computational Physics*, vol. 6, no. 1, pp. 136–141, 1970.

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

[139] 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. 2, 2024.

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

[141] 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, no. 18, pp. 7018–7036, 2011.

[142] 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, p. 110934, 2022.

[143] A. B. Langdon, B. I. Cohen, and A. Friedman, “Direct implicit large time-step particle simulation of plasmas,” *Journal of Computational Physics*, vol. 51, no. 1, pp. 107–138, 1983.

[144] 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, no. 2, pp. 231–247, 1995.

[145] 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. 10, 2023.

[146] 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, no. 5, pp. 1809–1841, 2022.
[147] 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. 15, 2021.

[148] 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. 7, 2024.

[149] 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. 19, 2022.

[150] 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. 17, p. 175201, 2023.