# Impact of the radial electric field on plasma instabilities in a penning discharge: Insights from a 3D particle-in-cell model

G Fubiani, Yann Camenen, F Gaboriau, L Garrigues, M Papahn Zadeh, A Smolyakov, M Tyushev

► To cite this version:

G Fubiani, Yann Camenen, F Gaboriau, L Garrigues, M Papahn Zadeh, et al.. Impact of the radial electric field on plasma instabilities in a penning discharge: Insights from a 3D particle-in-cell model. Physics of Plasmas, 2025, 32, pp.080704. 10.1063/5.0271375. hal-05373368

HAL Id: hal-05373368

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

Submitted on 19 Nov 2025

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 Attribution 4.0 International License
# Impact of the radial electric field on plasma instabilities in a penning discharge: Insights from a 3D particle-in-cell model

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

G. Fubiani 🧮 ; Y. Camenen 🧮 ; F. Gaboriau 🧮 ; L. Garrigues 🧮 ; M. Papahn Zadeh 🧮 ; A. Smolyakov 🧮 ;
M. Tyushev 🧮

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

Check for updates

Phys. Plasmas 32, 080704 (2025)

https://doi.org/10.1063/5.0271375

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

View
Online

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

Export
Citation

## Articles You May Be Interested In

* Mode transitions and spoke structures in $E \times B$ Penning discharge

* *Phys. Plasmas* (January 2025)

* Full-scale modeling and experimental study of a gas neutron tube with a Penning ion source

* *J. Appl. Phys.* (October 2024)

* Azimuthal structures and turbulent transport in Penning discharge

* *Phys. Plasmas* (March 2023)

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

<figure><img src="image_7.png" /></figure>
# Impact of the radial electric field on plasma instabilities in a penning discharge: Insights from a 3D particle-in-cell model

F

Cite as: Phys. Plasmas **32**, 080704 (2025); doi: 10.1063/5.0271375

Submitted: 17 March 2025 · Accepted: 7 August 2025 ·

Published Online: 26 August 2025

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

View Online

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

Export Citation

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

CrossMark

G. Fubiani,<fnref n="1" />,ᵃ⁾ **ID** Y. Camenen,<fnref n="2" /> **ID** F. Gaboriau,<fnref n="1" /> **ID** L. Garrigues,<fnref n="1" /> **ID** M. Papahn Zadeh,<fnref n="3" /> **ID** A. Smolyakov,<fnref n="3" /> **ID** and M. Tyushev<fnref n="3" /> **ID**

## AFFILIATIONS

## ABSTRACT

This Letter describes the main features of the azimuthal and axial structures in a numerical model of a 3D low-pressure magnetized plasma column (a so-called Penning discharge). The simulations are performed with a 3D/3V explicit momentum conserving Particle-In-Cell (PIC) algorithm with Monte Carlo Collisions (MCC). The discharge is seeded by an electron beam injected on one end of the cylinder, ionizing helium atoms for gas pressures of mTorr-scale. We study how the magnetic field strength and a bias voltage on the column's end plates affect the radial electric field profile and the formation of rotating plasma instabilities. We observe kHz-scale rigid rotating structures localized azimuthally (so-called spokes) in the non-linear saturated regime when $E \cdot \nabla n > 0$ and large m-modes azimuthal waves in the $\sim 1$ MHz range when $E \cdot \nabla n < 0$. Axial modes in the direction parallel to the magnetic field lines are also observed, which are weakly correlated with azimuthal modes but related to the non-linear plasma-beam instabilities.

© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0271375

Penning-type linear devices use cylindrical configurations with the magnetic field parallel to their main axis.<fnref n="1" />⁻⁵ The external electron source may be provided by a beam of energetic electrons injected on one end of the cylinder with a given energy or low-energy thermal electrons emitted from an electrode. In both cases, these electrons are accelerated by the DC plasma sheath (from an externally applied bias voltage), thus providing energy to sustain the discharge through ionizing collisions with the background gas. Typically, the beam of primary electrons is narrow compared to the radius of the linear device. A typical regime of interest is when the electrons are strongly magnetized, $\omega_{ce} \gg \nu_{en}$ while depending on plasma parameters, ions may be only weakly magnetized, $\omega_{ci} \le \nu_{in}$, where $\omega_{ce}$ ($\omega_{ci}$) is the electron (ion) cyclotron frequency and $\nu_{en}$ ($\nu_{in}$) the electron-neutral (ion-neutral) collision frequency, respectively.

An axial magnetic field improves the electron confinement, increasing the electron residence time, thus enhancing ionization and leading to higher plasma densities. Electron confinement is further improved by the bias potential on the column's end plates, which traps

the electrons axially. Such configurations are relatively simple but efficient sources of dense plasmas used for numerous applications, e.g., material processing,⁶ nuclear detectors,⁷ mass separation devices,⁸ and many others.

A variety of different types of oscillations typically occur in Penning sources.<fnref n="2" />⁻⁵ These oscillations are generally associated with drift-type waves due to radial gradients of plasma parameters.⁹⁻¹² Relative streaming of ions and electrons is one common mechanism of such instabilities; in particular, azimuthal modes can be generated by the difference in drift velocities between electrons and ions. Another potential source of fluctuations is the presence of the electron beams and other inhomogeneities in phase space, such as flows of ions accelerated by the electric field, e.g., in the sheath region. Simon-Hoh Instability (SHI) is one example of perturbations observed in partially magnetized plasmas⁹⁻¹¹ and in the long wavelength limit $k_{\perp} \rho_e \ll 1$ when electron inertia may be neglected, where $k_{\perp}$ is the wavevector in the plane perpendicular to the magnetic field lines, $\rho_e = v_{the}/\omega_{ce}$ is the electron gyroradius, and $v_{the}$ is the thermal velocity, respectively. It

<footnote n="1">¹LAPLACE, Université de Toulouse, CNRS, Toulouse, France</footnote>

<footnote n="2">²Aix-Marseille Université, CNRS, PIIM UMR 7345, Marseille, France</footnote>

<footnote n="3">³Department of Physics and Engineering Physics, University of Saskatchewan, Saskatoon SK S7N 5E2, Canada</footnote>

<footnote>ᵃ) Author to whom correspondence should be addressed: gwenael.fubiani@cnrs.fr</footnote>
occurs when the electric field is pointing in the same direction as the density gradient, so that $E_r \partial n / \partial r > 0$, corresponding to a typical situation where $\partial n / \partial r < 0$ and the plasma potential confines the ions radially to enforce quasi-neutrality (due to the mobility of unmagnetized ions which exceeds that of the magnetized electrons). Additional effects such as electron inertia and collisions result in more general gradient-drift instabilities mechanisms, e.g., lower hybrid modes. The linear regime of gradient-drift instabilities may be studied in the local approximation¹² where the length scale of the equilibrium solution is assumed much larger than the considered instability wavelength. This is the so-called Boussinesq approximation. Note that the local criterion is often modified when including the global profiles into the calculation, so a more complex nonlocal analysis¹³ is required, while, in addition, the observed fluctuations in the nonlinear saturated regime could be significantly different from predictions of the linear theory. Such differences are especially apparent for large-scale azimuthally rotating structures (so-called spokes) which have been observed both numerically and experimentally²⁻¹⁴⁻¹⁸ in various working conditions of Penning discharges where typically the small scale have the largest growth rates compared to the low *m* mode. A nonlinear energy transfer from small to large scales (inverse cascade) may be involved in the formation of the dominant low *m* modes, but this is still an open problem. Therefore, the study of complex nonlinear physics in realistic geometries and taking into account the self-consistent ionization effects and boundary conditions, requires nonlinear numerical simulations.

The scope of this article is to model numerically in a 3D geometry the low-pressure (mTorr-scale) helium Penning discharge self-consistently supported by the ionization from the electron beam and taking into account realistic boundary conditions in the radial and axial directions. The focus is on the investigation of the azimuthal and axial wave patterns as a function of the magnetic field, electron beam energy, and applied bias voltage.

The two-dimensional numerical studies of Penning discharges have demonstrated the excitations of large-scale $m = 1$ spoke modes and smaller scale $m > 1$ azimuthal structures.⁷,⁸,¹⁴,¹⁵,¹⁹⁻²² There are only a few preliminary studies of spoke structures in full three-dimensional geometries at various conditions,⁶,⁷,²³,²⁴ which are much more demanding and require more computer resources and more efficient numerical codes. In this work, we present one of the first comprehensive studies of instabilities in a full 3D geometry for a Penning discharge, exploring a range of plasma conditions. Our emphasis is on the investigation of the transition between the large-scale (low *m*) spoke modes and smaller scale *m* ≫ 1 spiral arm structures. We also look for the manifestations of the axial instabilities that may be excited due to the energetic electron beam injection and by the sheath potential at the biased ends. The effects of the externally applied bias voltage and the resulting electric field on the character of the excited structures are of particular interest. Two-dimensional calculations¹⁴,¹⁵ have shown that while *m* = 1 spokes appear when the radial electric field is negative so that Simon-Hoh criterion is satisfied, small-scale modes are excited instead by a positive electric field, and hence the destabilization mechanism must be different in that case. Investigations of this trend in full 3D geometry are one of the main results of this paper. In practice, a transition from a radially inward to a radially outward electric field is obtained by either changing the strength of the magnetic field or the bias voltage on the end plates of the cylinder. In addition,

two-dimensional studies are limited to the radial-azimuthal plane and do not provide information on the axial mode structures and their role on the azimuthal modes. Preliminary 3D studies²⁴ have shown that both axially shearless modes and helical structures may be excited. We demonstrate that there is little correlation between the axial and azimuthal structures so the dominant low *m* azimuthal modes are of the flute type, i.e., have no axial shear.

Figure 1 shows the simulation domain for the linear device configuration studied in this work. The diameter of the cylinder is $D = 10$ cm and length $L_z = 53.3$ cm. The electron beam is injected at the bottom and is modeled as a shifted half-Maxwellian flux along (Oz) with a kinetic energy $E_b$ and temperature $T_b = 5$ eV. The top and bottom end plates are biased at a potential $V_{ep}$, and the side of the cylinder is maintained at ground potential (0 V). The end plates are biased to control electrostatically the axial plasma confinement. We perform parameter scans from ground potential to -40 V and we adjust the initial beam energy such that $E_b - V_{ep} = 55$ eV to approximately account for the fact that the primary beam of electrons is accelerated by the plasma sheath. The radial distribution of the injected beam is assumed uniform with a

<figure><img src="image_11.png" /><figcaption>FIG. 1. Schematic representation of the simulation domain. The electron beam is injected at the bottom (shifted half-Maxwellian flux) with a radius $R_b = 1.5$ cm. The end-plates at both ends of the column are biased at a potential $V_{ep}$. The side of the cylinder is at ground potential. The dimensions of the cylinder are a radius $R = 5$ cm and length $L_z = 53.3$ cm. The coordinate of the column's center is at $x = y = R$ and $z = L_z/2$, respectively.</figcaption></figure>
radius $R_p = 1.5$ cm. The magnetic field is aligned with the axis of the cylinder in the positive z-direction and has a uniform profile, $\mathbf{B} = B_z \mathbf{e}_z$, where $\mathbf{e}_z$ is a unit vector and $B_z > 0$. Its strength is varied between 25 and 100 G in the calculations, enough to magnetize the electrons but not the ions so that $\omega_{ci} \le \nu_{in}$.

The numerical algorithm is 3D in space, 3-velocities (3 V), electrostatic explicit momentum conserving PIC/MCC parallelized with the OpenMP and MPI libraries.²⁵ The simulations are performed on a Cartesian grid with a circular boundary in the xy-plane (staircase approximation). For the model of Fig. 1, the numerical resolution is $192^2 \times 1024$ grid nodes, corresponding to a grid spacing $\Delta x/\lambda_{De} \le 0.5$ with $\Delta x = \Delta y = \Delta z$ and a (normalized) time step $\omega_{pe}\Delta t \approx 0.1$ where $\lambda_{De}$ is the electron Debye length and $\omega_{pe}$ the electron plasma frequency, respectively. The average plasma density is about $n_p \sim 10^{14} \text{m}^{-3}$ associated with a number of macroparticles per cell (ppc) between 10 and 20 (this corresponds to a minimum number of macroparticles within a box of volume one electron Debye cube $N_d = n_{em}\lambda_{De}^3 \sim 100$ in the simulations, where $n_{em}$ is the macro-electron density, respectively). One of our goals is to study the effect of the radial electric field on the excitation of the azimuthal modes, maintaining approximately the same plasma density. Thus, the electron beam current $I_b$ is adjusted accordingly for all studied cases. The current magnitude ranges between 5 mA for the optimum plasma confinement ($V_{ep} = -40$ V) and 40 mA for the lowest ($V_{ep} = 0$ V). The simulations are initiated without any background plasma. The working gas is helium at room temperature with a pressure of 3.5 mTorr, i.e., a neutral density $n_g = 10^{20} \text{m}^{-3}$. The physical-chemistry includes interactions between charged particles and neutrals: elastic, inelastic, and ionization collisions for electrons as well as elastic and backscattering reactions for He⁺ ions. Cross-sections are taken from the LXCAT database.²⁶,²⁷ The frequency for the collisions between electrons and He⁺ averaged over the entire simulation domain is found to be of the order of $\nu_{ei} \approx 5 \times 10^3 \text{s}^{-1}$, i.e., much smaller than the total collision frequency against neutrals ($\nu_{en} \approx 7.5 \times 10^6 \text{s}^{-1}$). Thus, Coulomb collisions are negligible compared to scattering against neutrals due to the low plasma density. We assume that macroparticles are isotropically distributed after a collision in this model and that the kinetic energy is shared equally between electrons after an ionizing event. This may overestimate the number of beam electrons trapped by the plasma potential. An improvement to the model would be to account for the small deflection angles of fast particles.²⁸,²⁹ Finally, primary electrons are modeled as a shifted half-Maxwellian flux distribution function moving in the positive z-direction,

$$f(v_z) = v_z \exp \left[ \frac{-(v_z - v_b)^2}{v_{th}^2} \right], \quad (1)$$

where $v_b = \sqrt{2eE_b/m_e}$ is the beam velocity, $m_e$ is the electron mass, $v_{th} = \sqrt{2eT_b/m_e}$ is the electron thermal velocity, and $e$ is the elementary charge; the beam kinetic energy and temperature are both expressed in electron volts ($e$ being the electron charge in Coulomb). Note that a typical case in these simulations was run using one computer node on two sockets, each with a 16-core 2.9 GHz Intel Xeon Gold 6326 processor, 32 threads in total (16 OpenMP times 2 MPI threads) and 256 Gb of Random Access Memory (RAM). One complete 3D PIC/MCC cycle per time step $\Delta t = 2 \times 10^{-10}$ s took 4.6 s, leading to a run-time of 6.4 h per micro-second of physical time, i.e.,

<figure><img src="image_12.png" /><figcaption>FIG. 2. Time-averaged plasma potential in the center of the plasma column and ion confinement time $\tau_i$ as a function of the bias voltage for $B_z = 100$ G, $E_b - V_{ep} = 55$ eV. The beam current is $I_b = 10$ mA, except at $V_{ep} = 0$ V ($I_b = 40$ mA) and $V_{ep} = -40$ V ($I_b = 5$ mA), to maintain a similar average plasma density.</figcaption></figure>

about 27 days to complete 100 µs, respectively. We used $1.2 \times 10^9$ macroparticles and about $3.8 \times 10^7$ grid nodes in the calculation.

Figure 2 shows the time-averaged plasma potential (for a time window greater than 100 µs) and ion confinement time in the center of the plasma column ($x = y = 5$ cm and $z = 26.6$ cm) vs the end plate bias voltage for $B_z = 100$ G. The beam energy is varied from $E_b = 55$ eV for $V_{ep} = 0$ to $E_b = 15$ eV for $V_{ep} = -40$ V keeping $E_b - V_{ep} = 55$ eV. The plasma confinement is improving with the lowering of the bias potential $V_{ep} < 0$. The ion confinement time $\tau_i$ is calculated from the volume-averaged ion continuity equation¹⁶

$$\tau_i = \frac{e\pi R^2 L_z \langle n_i \rangle}{I_{iw}}, \quad (2)$$

where $I_{iw}$ is the ion current collected on all walls, and $\langle n_i \rangle$ is the ion density averaged over the whole volume. Note that $n_i \approx n_e$ as the plasma is quasi-neutral (except in the sheath).

With a fixed magnetic field strength, applying a negative end plate bias voltage improves axial plasma confinement and decreases the plasma potential, eventually making it negative. As a result, the radial electric field becomes negative (inward), as shown by the time-averaged potential $\langle \phi_c \rangle$ in Fig. 2. Magnetized plasma electrons generated by ionization (and some fraction of beam electrons scattered by collisions and axial potential fluctuations) are trapped by the potential, oscillate back-and-forth along field lines, and diffuse outward slowly due to collisions with neutrals and anomalous transport by instabilities. A similar behavior is observed with a significant improvement in plasma confinement caused by a larger magnetic field, where $\tau_i$ doubles from 21 µs to about 42 µs while $\langle \phi_c \rangle$ evolves from 1.6 V to -1.4 V when $B_z$ is increased from 25 G to 100 G together with a fixed bias $V_{ep} = -30$ V. The transition to a negative potential and hence to an inward-pointing electric field results in a noticeable transformation of the azimuthal instability wave pattern.

This trend is shown in Fig. 3 for a scan of the magnetic field strength. At the low field value of $B = 25$ G and positive potential in the center, the plasma density and potential show relatively weak fluctuation amplitudes of large $m \approx 8-10$ azimuthal modes, Figs. 3(a) and 3(b). The wave pattern displays more complex structures with a higher magnetization, $B = 50$ G, where the plasma exhibits both low
<figure><img src="image_13.png" /><figcaption>FIG. 3. Formation of azimuthal structures as the sign of the potential in the center flips to negative with an increase in the magnetic field. Plasma density (left) and potential (right) are shown in the saturated regime for $V_{ep} = -30$ V and $I_b = 10$ mA; (a-b) – $B_z = 25$ G, $\langle\phi_c\rangle = 1.6$ V, $E_b = 25$ eV, $n_{max} = 4 \times 10^{14}$ m$^{-3}$; (b-c) – $B_z = 50$ G, $\langle\phi_c\rangle = -0.4$ V, $E_b = 12$ eV, $n_{max} = 5 \times 10^{14}$ m$^{-3}$; (e-f) – $B_z = 100$ G, $\langle\phi_c\rangle = -1.4$ V, $E_b = 25$ eV, $n_{max} = 10^{15}$ m$^{-3}$.</figcaption></figure>

and high *m* modes rotating counterclockwise, Figs. 3(c) and 3(d). In this case, the averaged central potential becomes slightly negative, $\langle\phi_c\rangle \simeq -0.4$ V, so that the conditions are favorable for the gradient-drift Simon-Hoh-type destabilization mechanism occurring when the radial electric field $E_r$ and the plasma density gradient are in the same direction,¹²,³⁰ $E_r\partial n/\partial r > 0$. Increasing the magnetic field to $B = 100$ G results in a more negative electric field and, as shown in Figs. 3(e) and 3(f), the generation of a large-scale $m = 2$ azimuthal structure (spoke) rotating with a frequency of 36 kHz (period of 27 µs) coexisting with a significantly shorter wavelength ($m \gg 1$). This behavior is consistent with the results of 2D (azimuthal-radial) simulations.¹⁴,¹⁵ It is important to note that for large $E_r\partial n/\partial r > 0$ values, standard linear fluid models with magnetized electrons and non-magnetized ions predict the most unstable mode to be at scales $m \gg 1$, comparable to that of the electron Larmor radius, where gyroviscous effects start to affect the stability¹² ($k_\perp\rho_e \sim 1$). In contrast, in the non-linear simulations, a low *m* rotating spoke is observed.¹⁵ An inverse cascade process from high *m* linearly unstable modes is a possible route for the formation of these low *m* structures.

The formation of low *m* modes with a radial electric field pointing inward has also been observed in our simulations when biasing negatively the end plates while keeping the magnetic field amplitude fixed at 100 G. Figure 4 plots the plasma density and ion current density profiles in the xy-plane, in the middle of the column, for bias

<figure><img src="image_14.png" /><figcaption>FIG. 4. Plasma density (left) and ion current density (right) in the saturated regime for $B_z = 100$ G; (a-b) – $E_b = 55$ eV, $V_{ep} = 0$ V, $\langle\phi_c\rangle = 12$ V, $I_b = 40$ mA, $n_{max} \simeq 7.5 \times 10^{14}$ m$^{-3}$, $\Gamma_{max} = j_{max}/e \simeq 2.5 \times 10^{18}$ m$^{-2}$s$^{-1}$; (c-d) – $E_b = 15$ eV, $V_{ep} = -40$ V, $\langle\phi_c\rangle = -1.4$ V, $I_b = 5$ mA, $n_{max} \simeq 6 \times 10^{14}$ m$^{-3}$ and $\Gamma_{max} \simeq 1.4 \times 10^{18}$ m$^{-2}$s$^{-1}$. The green streamlines show the directions of the ion flux.</figcaption></figure>

voltages of 0 V and -40 V with $E_b - V_{ep} = 55$ eV (recall that the beam current strength was adjusted to maintain similar plasma density values). We performed a benchmark which provided similar results using twice the number of macroparticles and a time step such that $v_{b,max}\Delta t/\Delta x \simeq 0.8$ where $v_{b,max}$ corresponds to a beam electron with a kinetic energy $E_k = E_b + 3T_b$. Note that the cyclotron frequency for the singly charged helium ions is $\omega_{ci} \simeq 2.4 \times 10^5$ s$^{-1}$. The total ion-neutral frequency is calculated as $\nu_{in} \simeq 2 \times 10^5$ s$^{-1}$ giving a Hall parameter $h_i = \omega_{ci}/\nu_{in} \simeq 1.2$, respectively, so that the ions are rather collisional and weakly affected by the magnetic field. For the

<figure><img src="image_15.png" /><figcaption>FIG. 5. (a) Normalized ion density profile along the plasma column, in the yz-plane, and measured at $x = 6.5$ cm for the conditions of Figs. 4(a) and 4(b). The wave pattern is somewhat a flute mode for $z \ge 15$ cm. (b) Normalized ion source term $S_i$ showing a footprint of the beam-plasma instability along the beam path. $n_{max} = 6.4 \times 10^{14}$ m$^{-3}$ and $S_{max} = 3 \times 10^{20}$ m$^{-3}$ s$^{-1}$.</figcaption></figure>
parameters of case (a-b) in Fig. 4, we observe an instability in the form of spiral arms. The plasma potential is positive with respect to the grounded cylinder walls in this case, hence the radial electric field is pointing outward (positive), and we are in the conditions where $E_r \partial n / \partial r < 0$. The plasma wave exhibits a large mode number $m = 10$ rotating in the direction of the electron diamagnetic drift (counterclockwise). Ions diffuse to the wall and are weakly dragged azimuthally in the direction of the $E \times B$ drift, forming a rather regular pattern for the spiral arms, which is clearly pronounced in the ion flux density, as seen in Fig. 4(b).

The criterion for the SHI is not satisfied in this case, and the observed fluctuations have an azimuthal rotation frequency of $f \simeq 1.3$ MHz, which can be identified as the lower-hybrid (LH) mode destabilized by the density gradient and collisions (the electron cyclotron frequency is $\omega_{ce} = eB_z/m_e \simeq 1.75 \times 10^9$ s$^{-1}$, so that the LH frequency $\omega_{LH} = \sqrt{\omega_{ce}\omega_{ci}} \simeq 2 \times 10^7$ s$^{-1}$, and hence $f_{LH} = 3.3$ MHz). The LH instability includes the effects of the electron inertia leading to small-scale perturbations.¹² The mode is mostly flute along the plasma column, so that the spiral arms display a rigid rotational pattern, i.e., the position of the wavefronts is somewhat following a given magnetic field line but with a slight twist as shown in Fig. 5(a), which plots the normalized plasma density in the yz-plane at $x = 6.5$ cm. The first maxima on the left-hand side (LHS), for instance, is offset by about 5 mm along (Oy) over the upper half of the plasma column.

Next, for the parameters of case (c-d) in Fig. 4, when $V_{ep} = -40$ V and $E_b = 15$ eV, the potential in the center becomes negative, the radial electric field is pointing inward, and hence both the $E \times B$ and diamagnetic drift velocities point in the same counterclockwise rotation direction. This is a situation similar to the case of Figs. 3(e) and 3(f), where the necessary condition for the onset of the SHI is verified ($E_r \partial n / \partial r > 0$) resulting here also in a low-frequency large-scale $m = 2$ (spoke) mode. The ion flow is strongly distorted and driven radially to the cylinder walls by the time-varying rotating spoke electric field; it represents a significant fraction of the total current $I_{iw}$ (about 38%) as shown in Table I. More insight on the mode structure in provided by measuring the time variation of the electron density in the mid-plane $z = 26.6$ cm ($L_z/2$), at $x = 1$ cm and $y = 0$, respectively, as shown in Fig. 6. The spoke is formed around $t \simeq 50$ µs (the simulation is initiated without any background plasma, just the electron beam) and begins with a mode $m = 1$ exhibiting a rigid rotation at a frequency $\sim 22$ kHz until about $t \simeq 260$ µs which shows a bifurcation toward $m = 2$ (at twice the frequency). Note that the mode of Figs. 3(e) and 3(f) is in a transition phase toward $m = 3$. A

<table><caption>TABLE I. Total ion current $I_{iw}$ collected on the Penning discharge surfaces, normalized ion currents on the cylinder wall $I_{ic}/I_{iw}$ and the two end-plates $I_{ep}/I_{iw}$ vs the end-plate bias voltage $V_{ep}$. $I_b$ is the beam current. The magnetic field strength is fixed at 100 G.</caption><thead><tr><th>V<sub>ep</sub> (V)</th><th>I<sub>b</sub> (mA)</th><th>I<sub>iw</sub> (mA)</th><th>I<sub>ic</sub>/I<sub>iw</sub></th><th>I<sub>ep</sub>/I<sub>iw</sub></th></tr></thead><tbody><tr><td>0</td><td>40</td><td>9.1</td><td>0.85</td><td>0.15</td></tr><tr><td>-10</td><td>10</td><td>2</td><td>0.62</td><td>0.38</td></tr><tr><td>-20</td><td>10</td><td>2.2</td><td>0.37</td><td>0.63</td></tr><tr><td>-30</td><td>10</td><td>2.75</td><td>0.38</td><td>0.62</td></tr><tr><td>-40</td><td>5</td><td>1.6</td><td>0.38</td><td>0.62</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>FIG. 6. Electron density vs time measured numerically at the mid-plane of the plasma column $z = L_z/2$, $x = 1$ cm, and $y = 0$ for the conditions of Figs. 4(c) and 4(d). The main peaks are associated with the spoke crossing the probe location. It shows a constant rotation frequency with a mode $m = 1$ for $50 \le t \le 260$ µs splitting to $m = 2$ beyond that time window. Small-scale structures are also visible on the signal. $B = 100$ G.</figcaption></figure>

similar behavior and intermittent bifurcations between the low $m$ modes were also observed in the experiments.¹⁷,¹⁸ Instabilities are typically understood to act so as to diminish their original driving sources. The observed bifurcation to a higher $m$ mode may be attributable to the evolution of the background equilibrium toward a state of increased but marginal stability.

One important result of our simulations is the observation of the axial mode, due to a beam-plasma instability propagating along the direction of the beam path, as shown in Fig. 7. The wave pattern exhibits strong potential variations up to $\Delta\phi \sim 20$ V for grounded walls [case study of Figs. 4(a) and 4(b)] leaving a footprint only on the electron density (the ion density is somewhat uniform) with $k_z \lambda_{De} \sim 0.1$. A similar pattern is also visible on the ion source term $S_i = n_i v_i$.

<figure><img src="image_17.png" /><figcaption>FIG. 7. Axial electric field along the z-direction vs time in the center of the plasma column and for grounded walls (including the end-plates). The beam-plasma instability results in an axially propagating wave with a phase velocity consistent with the beam velocity. Case study of Figs. 4(a) and 4(b).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIG. 8. Plasma potential profile along (Oz) at the center of the column $x = y = 5$ cm for grounded walls and a –40 V end-plate bias voltage, respectively. Test cases of Fig. 4.</figcaption></figure>

displayed in Fig. 5(b), where $v_i$ is the ionization frequency. For these parameters, $n_b = I_b/(ev_b\pi R_b^2) \simeq 8 \times 10^{13} m^{-3}$ (given at the pre-sheath entrance after injection, i.e., before being scattered away by elastic and inelastic collisions with neutrals) with $I_b = 40$ mA, $v_b = \sqrt{2eE_b/m_e} \simeq 4.8 \times 10^6$ m/s and $E_b \simeq 67$ eV on average in the plasma bulk. The spatial scale of the axial modes is consistent with the resonant condition for the beam-plasma instability $\omega_{pe} = k_z v_b$. We measure in the center of the column a plasma frequency $f_{pe} \simeq 163$ MHz and a wavelength $\lambda_{pe} \simeq 2.9$ cm, giving an estimate for the beam velocity $v_b^* = f_{pe}\lambda_{pe} \simeq 4.7 \times 10^6$ m/s and hence $v_b^* \simeq v_b$. Note that the axial fluctuation amplitude is significantly smaller for a lower beam current, as shown for the case of a –40 V end plate bias voltage in Fig. 8, where $I_b = 5$ mA and $E_b \simeq 54$ eV are in the plasma.

In conclusion, we modeled numerically in 3D the generation of plasma instabilities in a Penning discharge using a 3D/3V Cartesian explicit momentum conserving PIC/MCC algorithm. The amplitude of the plasma potential is modified by acting on the plasma confinement, i.e., either through changing the magnetic field strength or biasing negatively the end plates. We find that when the plasma potential in the plasma bulk is positive with respect to the walls of the device, then the wave pattern exhibits a high mode number $m \gg 1$ rotating azimuthally counterclockwise at frequencies of $\sim 1$ MHz comparable to the LH frequency. A spoke appears in the non-linear saturated regime, i.e., a large-scale rotating structure, when the plasma potential confines the ions (its amplitude is negative, below the value on the walls, which are grounded) and the SHI criterion is satisfied. The rotation frequency is about tens of kHz, and the mode number is of order unity ($m \le 3$ in our simulations). Small wavelength structures are superimposed on the spokes. It is expected that these small-scale, lower hybrid-type oscillations result in ion heating and may also play a role in anomalous electron transport. Assessing their actual contribution under the current conditions is left for future work. In all cases, the direction of rotation is that of the electron diamagnetic drift, and the modes are close to flute along the axis of the plasma column. Finally, a short-wavelength beam-plasma instability develops along the direction of the beam propagation.

## SUPPLEMENTARY MATERIALS

See the supplementary material for a video of a rotating spoke for G and an end plate potential of V in a helium gas discharge. This case study corresponds to Figs. 4(c) and 4(d) before the transition to an $m=2$ mode (see Fig. 6).

This work has received funding from the French Fédération de Recherche sur la Fusion Magnétique (FR-FCM).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Gwenael Fubiani:** Conceptualization (lead); Data curation (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (equal); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Yann Camenen:** Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Validation (supporting); Writing – review & editing (supporting). **Freddy Gaboriau:** Formal analysis (supporting); Investigation (supporting); Validation (supporting); Writing – review & editing (supporting).

**Laurent Garrigues:** Formal analysis (supporting); Investigation (supporting); Validation (supporting); Writing – review & editing (supporting). **Mina Papahn Zadeh:** Investigation (supporting). **Andrei I. Smolyakov:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal).

**Mikhail Sergeevich Tyushev:** Writing – review & editing (supporting).

## DATA AVAILABILITY

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

## REFERENCES

¹C. Rebont, N. Claire, T. Pierre, and F. Doveil, Phys. Rev. Lett. 106, 225006 (2011).

²S. Jaeger, T. Pierre, and C. Rebont, Phys. Plasmas 16, 022304 (2009).

³N. Claire, A. Escarguel, C. Rebont, and F. Doveil, Phys. Plasmas 25, 061203 (2018).

⁴E. Rodriguez, V. Skoutnev, Y. Raitses, A. Powis, I. Kaganovich, and A. Smolyakov, Phys. Plasmas 26, 053503 (2019).

⁵C. Cheon, J. Choi, J. B.-W. Koo, and J. Y. Kim, Plasma Sources Sci. Technol. 32, 07LT01 (2023).

⁶A. Likhanskii, A. I. Smolyakov, M. Papahn Zadeh, N. Chopra, and Y. Raitses, in *76th Annual Gaseous Electronics Conference*, Vol. 68, ER2.00004.

⁷A. S. Rokhmanenkov and S. E. Kuratov, J. Phys. Conf. Ser. 1250, 012036 (2019).

⁸V. S. Smirnov, S. A. Kislenko, A. V. Gavrikov, A. A. Klimkov, G. D. Liziakin, and V. P. Smirnov, Plasma Sources Sci. Technol. 32, 095007 (2023).

⁹A. Simon, Phys. Fluids 6, 382 (1963).

¹⁰F. C. Hoh, Phys. Fluids 6, 1184 (1963).

¹¹Y. Sakawa, C. Joshi, P. Kaw, F. Chen, and V. Jain, Phys. Fluids B: Plasma Phys. 5, 1681 (1993).

¹²A. I. Smolyakov, O. Chapurin, W. Frias, O. Koshkarov, I. Romadanov, T. Tang, M. Umansky, Y. Raitses, I. D. Kaganovich, and V. P. Lakhin, Plasma Phys. Controlled Fusion 59, 014041 (2017).

¹³S. Aggarwal, Y. Camenen, A. Escarguel, and A. Poyé, J. Plasma Phys. 89, 905890310 (2023).
<footnote n="14">¹⁴M. Tyushev, M. Papahn Zadeh, N. S. Chopra, Y. Raitses, I. Romadanov, A. Likhanskii, G. Fubiani, L. Garrigues, R. Groenewald, and A. Smolyakov, *Phys. Plasmas* **32**, 013511 (2025).</footnote>

<footnote n="15">¹⁵M. Tyushev, M. Papahn Zadeh, V. Sharma, M. Sengupta, Y. Raitses, J.-P. Boeuf, and A. Smolyakov, *Phys. Plasmas* **30**, 033506 (2023).</footnote>

<footnote n="16">¹⁶G. Fubiani, R. Agnello, I. Furno, L. Garrigues, P. Guittienne, G. Hagelaar, A. Howling, R. Jacquier, B. Reman, A. Simonin, and F. Taccogna, *Phys. Plasmas* **28**, 063503 (2021).</footnote>

<footnote n="17">¹⁷M. A. Valinurov, R. A. Usmanov, A. V. Gavrikov, A. P. Oiler, and V. P. Smirnov, *Phys. Lett. A* **534**, 130242 (2025).</footnote>

<footnote n="18">¹⁸M. Lee, K.-J. Chung, K. Hara, and J. Y. Kim, *Phys. Plasmas* **31**, 123508 (2024).</footnote>

<footnote n="19">¹⁹J. P. Boeuf and M. Takahashi, *Phys. Rev. Lett.* **124**, 185005 (2020).</footnote>

<footnote n="20">²⁰J.-P. Boeuf, *Front. Phys.* **2**, 74 (2014).</footnote>

<footnote n="21">²¹J. Carlsson, I. Kaganovich, A. Powis, Y. Raitses, I. Romadanov, and A. Smolyakov, *Phys. Plasmas* **25**, 061201 (2018).</footnote>

<footnote n="22">²²R. Lucken, A. Bourdon, M. A. Lieberman, and P. Chabert, *Phys. Plasmas* **26**, 070702 (2019).</footnote>

<footnote n="23">²³K. Matyash, R. Schneider, S. Mazouffre, S. Tsikata, and L. Grimaud, *Plasma Sources Sci. Technol.* **28**, 044002 (2019).</footnote>

<footnote n="24">²⁴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:2501.01161* (2025).</footnote>

<footnote n="25">²⁵G. Fubiani, L. Garrigues, G. Hagelaar, N. Kohen, and J. P. Boeuf, *New J. Phys.* **19**, 015002 (2017).</footnote>

<footnote n="26">²⁶Biagi-v7.1 database, http://www.lxcat.net, Retrieved November 18, 2024.</footnote>

<footnote n="27">²⁷Phelps database, http://www.lxcat.net, Retrieved November 18, 2024.</footnote>

<footnote n="28">²⁸K. Nanbu, *Plasma Sci., IEEE Trans.* **28**, 971 (2000).</footnote>

<footnote n="29">²⁹J. F. J. Janssen, L. C. Pitchford, G. J. M. Hagelaar, and J. van Dijk, *Plasma Sources Sci. Technol.* **25**, 055026 (2016).</footnote>

<footnote n="30">³⁰J.-P. Boeuf and A. Smolyakov, *Phys. Plasmas* **30**, 050901 (2023).</footnote>