Scaling of spoke rotation frequency within a Penning discharge

Hall plasmas, consisting of magnetized electrons and unmagnetized ions, exhibit a wide range of plasma instabilities.^1,2 In some systems, these instabilities may result in the formation of a long wavelength, low frequency fluctuation in plasma density, propagating in the E × B direction. This rotating structure is commonly referred to as a plasma “spoke.”

The spoke has been well characterized within a number of experimental devices, including Hall thrusters,^3–15 planar magnetrons,^16–23 cylindrical magnetrons,^24,25 and Penning discharges.^26–28 These devices feature cylindrical geometry, with electron drift and spoke propagation occurring in the azimuthal direction.

The plasma spoke may play an important role in the anomalous transport of electrons across the applied magnetic field, a source of reduced efficiency.²⁹ Using separated probes within a Hall thruster, Janes and Lowder³ measured an azimuthal electric field correlated with the passage of the spoke and suggested that the resulting E_θ × B_r drift was enhancing electron transport towards the thruster anode. More recently, the use of segmented anodes within a cylindrical Hall thruster³⁰ demonstrated that half of the discharge current was being carried through the spoke, evidence for enhanced transport through the structure.^6,8

Due to its highly non-linear, turbulent, and global nature, the spoke has continued to evade a clear theoretical understanding. Proposed mechanisms include a type of ionization wave,³ whereby the azimuthal field at the front of the spoke provides sufficient energy to ionize neutrals and propagate the plasma density perturbation. The ionization wave would propagate at the Critical Ionization Velocity (CIV)³¹ defined as $vciv=2Uion/mn$⁠, where U_ion and m_n are the ionization energy and mass of the neutral species, respectively.

Alternatively, the spoke may be the result of collective effects, with the Collisionless Simon-Hoh Instability (CSHI) being a likely candidate.^26,27 Charge separation between drifting electrons and unmagnetized (non-drifting) ions generates an azimuthal electric field. If the background electric field and density gradients within the system are aligned E₀ ⋅ ∇n₀ > 0, then the resulting azimuthal field will act to enhance perturbations in plasma density, driving the instability.

Investigation of the plasma spoke within a Penning discharge offers several advantages, most immediately being improved access for diagnostics. The applied magnetic field is ideally uniform and aligned with the axis of the device, allowing the electron dynamics parallel to the device axis to be decoupled from those in the transverse direction. For purposes of simulation and theory, the system can therefore be treated as approximately two-dimensional. The lack of a magnetic field gradient reduces the number of energy sources for instabilities, eliminating possible driving mechanisms for spoke formation. In many systems, the plasma is also weakly collisional, further simplifying theoretical considerations. This makes it an ideal system within which to study spoke formation and anomalous transport.

Simulations have played an increasingly important role in understanding the formation of the spoke and its connection to anomalous transport. Kinetic techniques are required to self-consistently capture transport effects, with the Particle-in-Cell Monte-Carlo Collision (PIC-MCC) method commonly being used. The challenge with PIC-MCC comes from the numerical requirement to resolve the smallest time and length scales associated with the plasma, the electron plasma frequency and Debye length, respectively. These constraints result in a significant computational overhead required to capture the vastly multi-scale physics observed within devices such as the Penning discharge.

This challenge is usually dealt with in one of two ways. The first approach is to scale the system in some way, by reducing device size (to reduce the grid size), decreasing the plasma density, or increasing the relative permittivity (to increase the cell size). With improving computational capabilities and employing scaling, it has finally become possible to simulate a 2D profiles of a device^32–38 and scaled 3D devices.^39–42 Low mode number, rotating structures have been observed in a number of these simulations; in particular, Matyash found close agreement between the frequency of the rotating structure observed within a model of a wall-less Hall thruster⁴³ and experiments at the French National Center for Scientific Research (CNRS).¹⁵ 

The second technique adopts a hybrid approach, modeling electrons as a fluid, such that the associated length and time scales need not be resolved. The primary disadvantage of this technique is that kinetic electron transport effects are not captured self-consistently and must be incorporated via a model informed from experimental evidence. None-the-less hybrid models have had success in reproducing the spoke. Most recently, a two-dimensional axial-azimuthal hybrid code captured the motion of the spoke within a simulated thruster channel with SPT-100 like parameters.⁴⁴ The spoke velocity was on the order of what has previously been observed in experiments, and suggested to be the result of an instability driven by the magnetic field and density gradients.^4,45

Despite these capabilities, significant work is required to produce unscaled, fully kinetic simulations of entire devices. Perhaps equally as challenging will be developing tools for analysing the enormous amount of data produced by such simulations.

The present paper is a follow up of the work commenced in Ref. 38 whereby large scale coherent structures were observed within two-dimensional PIC-MCC simulations of a Penning discharge cross-section. These simulations are modified to allow quasi-steady-state operation of the discharge, revealing formation of a single mode, azimuthally rotating structure. By comparing this structure with those observed within experiments, and studying how it scales with discharge properties, an attempt is made to determine the fundamental mechanism responsible for its formation.

Simulations were run using the Large-Scale Plasma code (LSP).⁴⁶ LSP is a multi-purpose and versatile PIC code, widely benchmarked, and validated within the community.^46–49 It was assumed that perturbations to the magnetic field resulted in negligible Lorentz forces when compared to those generated by electrostatic fields. Therefore, the system could be treated electrostatically, and Poisson's equation inverted to obtain the electric potential from charge density. PPPL modifications to the code⁵⁰ include incorporation of the latest version of the Portable Extensible Toolkit for Scientific Computation (PETSc)^51–53 for improved performance and scalability. Poisson's equation was inverted via the SuperLU,^54,55 LU factorization package accessed via the PETSc interface.

Simulation were modeled off the Penning discharge experiment at the PPPL Hall Thruster Experiment laboratory.²⁸ Within this device, an RF plasma cathode injects electrons along the axial magnetic field, ionizing a low pressure gas of either argon or xenon. Electron motion in the radial direction is inhibited by an axially applied uniform magnetic field. Ions are weakly magnetized and therefore significantly more mobile, giving rise to an ambipolar radial electric field. The applied magnetic field and ambipolar electric field result in electrons undergoing E_r× B_z drifts in the azimuthal direction. The plasma is surrounded by a grounded cylindrical metal anode with 10 cm diameter. The ends of the cylinder are dielectric, preventing the short-circuit effect. The device geometry and simulation geometry are shown in Fig. 1.

A slice of the device cross-section was modeled on a uniform Cartesian grid, with particles evolving in 2D-3V phase space. While a two-dimensional model appears to be a good approximation for this system, these simulations will not capture axial modes or the effects of the dielectric end plates. A Cartesian grid was chosen over a cylindrical grid to avoid numerical instabilities associated with the grid singularity at zero radius.

A uniform magnetic field is applied perpendicular to the simulation domain. To shed light on the formation mechanism of the spoke, both collisionless and collisional simulations were performed. In collisionless simulations, electrons and ions are injected into the center of the trap at fixed current and initial temperature. For collisional simulations, electrons are injected into a fixed, uniform background of neutrals, with ions forming via ionization. Modeled collisions include Coulomb collisions between all charged species and electron-neutral collisions. Neutral particle excitation was not modeled. Both electron-neutral elastic and ionizing collisions were informed by experimental data.⁵⁶ 

Simulations were designed to be completed within a realistic time frame, enabling parametric investigations of instabilities within the Penning discharge. To achieve this goal, the simulation domain was reduced from 10 cm to 5 cm and the xenon gas mass was reduced to that of Helium-4. To relax the strict constraints on PIC simulations, the relative permittivity was increased to ε_r = 400, increasing the Debye length and reducing plasma frequency, thereby allowing for larger cell size and time step, respectively. The relative permittivity was scaled, since it is suspected to play a small role in the instabilities which may be responsible for spoke formation. This scaling, however, will influence higher frequency instabilities, modifying the micro-structure of the discharge.

These modifications led to a grid of 250 × 250 cells, with a cell size Δx = 200 μm and a time step Δt = 40 ps, suitable to resolve all relevant plasma length and time scales. Simulations of 100 μs of plasma time could therefore be completed within 2 days. Simulations were run with 28 cores on the Princeton University Perseus supercomputer and 24 cores on the Department of Energy's National Energy Research Scientific Computing Center (NERSC) Edison supercomputer. Post-processing was performed using the Interactive Data Language based LSP post-processing tool P4 and in-house Python codes.

Electrons and ions are injected into the center of an initially empty domain with fixed currents and temperature I_e, T_e,inj and I_i, T_i,inj, respectively (see Table I).

The negative discharge current I_d = I_i – I_e < 0 results in an initially non-neutral plasma until sufficient ions have been injected to provide a neutralizing background. Quasi-neutrality is achieved around 100 μs, after which a long-wavelength, single mode structure forms, rotating in the azimuthal direction. When viewed in the direction of the applied magnetic field, the structure rotates in the anti-clockwise direction.

Figure 2 shows contour plots of instantaneous electron density, plasma potential, and current streamlines, commencing at 250 μs and incremented by a π/4 phase shift of the rotating structure. Figure 2(a) shows that the structure does not appear to rotate as a rigid body, but rather as a density perturbation rich in micro-instabilities and turbulence. This is further evidenced by Fig. 2(b) which shows a noisy and highly fluctuating plasma potential, with dips in plasma potential only weakly correlated with the rotating density perturbation. Figures 2(b) and 2(c) show that there is neither a strong electric potential gradient nor a current channel associated with the front of the spoke; however, both indicate the presence of a turbulent plasma. Qualitatively, it appears that increased cross-field electron transport within the structure is the result of enhanced turbulence, rather than a strongly correlated azimuthal electric field.

The phase of the rotating structure is measured to obtain a rotation frequency of f_s = 66.0 kHz, significantly lower than the electron and ion plasma frequencies, as well as the electron cyclotron frequency and the lower-hybrid frequency. When computed at r = R₀/2, the structure rotation velocity is v_r = 5.18 km/s, around half of the ion-acoustic velocity v_s = 10.6 km/s.

Numerical convergence was verified by measuring the mode frequency for simulations with half the cell size, half the time step, and double the particle number, showing no more than a 3.8% discrepancy. Modeling convergence was checked by scaling the relative permittivity, which for ε_r = 100 agrees within 11% (see Table II). This discrepancy may be explained by a modification of the plasma micro-structure with relative permittivity, which may in turn effect global properties and therefore spoke frequency [see, for example, Fig. 3(a)].

The discharge remains in a quasi-steady state, with the rotating structure persisting until the end of the simulation at 500 μs. This comprises approximately 26 full rotations, from which time averaged statistics of the turbulent structure can be computed. Average radial profiles are computed by first azimuthally averaging each sample and then temporally averaging over all samples. To compare the effects of the relative permittivity ε_r on the plasma properties, averaging is performed for values of ε_r = {100, 400}.

Figures 3(a) and 3(b) show the averaged radial profiles of electron density n_e(r) and gradient length scale $Ln(r)=ne(r)/(dne(r)/dr)$⁠. The plasma density is peaked at the center of the trap due to the source of electrons and ions, decaying away radially as the electrons are transported across the applied axial magnetic field. Increasing the relative permittivity results in a heightened density peak at the center of the trap, possibly due to weaker electric fields and a reduced rate of cross-field transport. The gradient-length-scale is negative at all radial positions.

Figures 3(c) and 3(d) show the averaged radial potential ϕ(r) and radial electric field E_r(r) profiles. Since the unmagnetized ions are more mobile than electrons across the field lines, an ambipolar field forms in the negative radial direction. This results in a potential well at the center of the domain acting to confine ions and thereby maintain quasi-neutrality. Where quasi-neutrality is maintained, the potential and electric field profiles remain similar for different relative permittivities. This demonstrates that, on a global scale, the electrostatic forces remain similar despite this scaling of the physics. This similarity breaks down near the sheath at the edge of the simulation, for r > 0.8R₀. Since the plasma is non-neutral in this region, the sheath size and shape are strongly influenced by the relative permittivity. No physical explanation can be given for the noisy signal of average radial electric field within r < 0.1R₀; however, it is most likely either a physical or a numerical artifact of particle injection.

The average electron temperature was similarly calculated and found to vary only weakly with radius, remaining near the injection temperature with an average $⟨Te⟩=4.7 eV$⁠. Ions were found to be heated from their original temperature, with an average temperature of $⟨Ti⟩=0.58 eV$⁠. This heating is the mechanism by which ions are eventually able to escape the potential well.

To obtain estimates for the E × B and diamagnetic drift velocities, the radial electric field and gradient length scales are averaged away from the injection and sheath regions within r ∈ [0.2R₀, 0.8R₀], giving an average azimuthal E × B velocity of $v0=⟨|Er|⟩/B0=13.1 km/s$ and electron diamagnetic drift velocity of $v*=Te/(e⟨|Ln|⟩B0)=64.5 km/s$⁠. Within the same frame of reference as Fig. 2, both of these drifts occur in the anti-clockwise direction. The rotation velocity of the large scale structure is an order of magnitude smaller than the electron diamagnetic drift velocity and less than half of the E × B velocity. Since the large scale structure is shown to be low-frequency, long-wavelength, and rotating in the E × B direction, it exhibits all of characteristic behaviour of the plasma spoke, as observed within the Penning discharge and numerous other experiments.^3,5,7,28 Therefore, it is proposed that the rotating structure observed within these simulations is the plasma spoke.

Current and density probes were placed at various azimuthal locations near the discharge anode. Figure 4 plots the local electron density n_e at (x_p, y_p) = (0, –0.996R₀) and the magnitude of the local electron current I_e,p passing through a horizontal chord which intersects (x_p, y_p). The peaks in electron density have the same frequency to that of the spoke shown in Fig. 2(a) and correlate with the passage of the structure. Peaks in local electron current are strongly correlated to peaks in electron density and therefore to the passage of the spoke, indicating that electron transport is enhanced by the spoke, in agreement with previous experimental observations.^6,8

Since a majority of the current is conducted through the spoke structure, the maximum magnitude of the large scale self-generated magnetic fields can be estimated. Conservatively assuming that the spoke is a column of current with the same radius as the injection region R_i, the maximum self-generated magnetic field is δB ≈ 0.2 G. Therefore, the maximum Lorentz force field generated by this perturbation interacting with drifting electrons is approximately $|vd×δB|=(v0+v*)|δB|≈2 V/m$⁠, or around 2% of the ambipolar electric field strength, demonstrating the validity of the electrostatic assumption.

The average cross-field electron conductivity σ_⊥e(r) can be calculated via Ohm's law with the average electric field and gradient length scale profiles from Sec. III A

where j_re(r) is the radial cross-field electron current density and d is the simulation depth (in LSP d = 1 cm).

The conductivity is related to an effective turbulent collision frequency ν_t(r) via

where ω_c = q_eB/m_e is the electron cyclotron frequency and it is assumed that ω_c ≫ ν_t.

Figure 5 plots the average turbulent collision frequency as a function of radius. Values are similar for both the ε_r = 100 and ε_r = 400 cases. Within r ∈ [0.2R₀, 0.8R₀], the mean turbulent collision frequency for the ε_r = 400 case is $⟨νt⟩=139.5 MHz$⁠, giving an effective Hall parameter β_eff = ω_c/ν_t = 12.6. This value is comparable to those measured by Bohm within numerous experiments exhibiting anomalous transport.⁵⁷ The decay of the collision frequency for r > 0.8R₀ is most likely due to non-neutral behaviour within the sheath.

To better replicate conditions of the experiment, the collisionless simulation from Sec. III A is modified to allow for ionization via MCC processes. The simulation domain initially consists of a fixed uniform background of neutral particles with pressure P_n and temperature T_n (see Table I). A uniform axial beam of electrons is injected with energy V_b and ions are formed via collisional processes between electrons and neutrals; note that neutrals are not depleted during this process.

The plasma is initially non-neutral, resulting in highly excited electrons which rapidly ionize the background gas and render the system quasi-neutral. A large scale, low-frequency single mode develops after 20 μs, rotating in the E × B direction. Figure 6 shows four contour plots of electron density, commencing at 61.4 μs and incremented by a π/4 phase shift of the rotating structure. This structure exhibits all of the same physical properties as in the collisionless case and is similarly proposed to be the plasma spoke.

The spoke frequency is computed as f_s = 62.4 kHz with a rotation velocity at r = R₀/2 of v_r = 4.90 km/s, both within 6% of the values obtained for the collisionless case. The structure of the density contours reveals a similar behaviour to that of the collisionless case, although the spoke structure extends further towards the anode and appears to have a larger azimuthal extent. This most likely indicates that ionization is occurring not only within the injection region at the center of the trap but also within the spoke itself, enhancing plasma density and broadening the spoke.

The average radial density and average gradient length scale profiles [see Figs. 7(a) and 7(b)] and average electric potential and average electric field profiles [see Figs. 7(c) and 7(d)] are computed and compared with the collisionless case. The density profile and density gradients for the collisional case are shallower than the collisionless case. This also indicates that ionization is likely occurring outside of the injection region, since a more diffuse source of ions can sustain a flatter electron density profile. The electric potential and electric field profiles are similar between both cases, with the collisional case exhibiting a more linear potential and therefore flatter electric field profile.

The profiles in electric field and gradient length scale are used to compute the average effective electron collision frequency ν_eff = 107.2 MHz and the effective Hall parameter β_eff = 16.4, similar in order of magnitude to the turbulent collision frequency calculated for the collisionless case. The electron-neutral collision frequency for Xenon with T_e ≈ 5 eV is ν_en = 12.2 MHz, giving a classical Hall parameter of β_c = 144.3. This indicates that the transport is highly anomalous and dominated by turbulent fluctuations.

Despite these differences, the spoke frequency and structure as well as the average discharge profiles are remarkably similar to the collisionless case. This indicates that the same fundamental mechanism is likely responsible for the formation of the spoke. Since it was shown in Sec. III A that ionization is not necessary for spoke formation, it is unlikely that the spoke is caused by an ionization wave. Furthermore for the reduced mass xenon, the CIV³¹ v_civ = 17.1 km/s, significantly larger than the spoke rotation velocity.

This leaves the Collisionless Simon-Hoh Instability (CSHI) as a likely candidate for the formation of the spoke within these simulations. In both the collisional and collisionless cases $Er(r)/Ln(r)>0 ∀ r$ such that the instability criterion for the CSHI is satisfied. Keeping in mind the limitations of these simulations, it is therefore possible that the rotating spoke observed within the Penning discharge is the result of the CSHI, rather than an ionization wave.

Considering the relative magnitudes of the average diamagnetic drift velocity v_* = 64.5 km/s, E × B velocity v₀ = 13.1 km/s, and ion-sound speed v_s = 10.6 km/s, linear CSHI theory suggests the following scaling for the spoke angular frequency:⁵⁸ 

Assuming a single azimuthal mode propagating at r = R₀/2, we have k = k_θ = 2/R₀, and a theoretical estimate for the spoke frequency

The validity of this scaling is tested by modifying the discharge parameters of the collisionless simulation introduced in Sec. III A. For each simulation, the radial electric field multiplied by the gradient length scale $⟨|ErLn|⟩$ is averaged for r ∈ [0.2R₀, 0.8R₀]. The estimate for spoke frequency f_s,th is plotted with the measured spoke frequency f_s against different discharge parameters.

The discharge current I_d is varied by modifying the ion injection current I_i. Figure 8 shows how the average radial density profile changes with discharge current magnitude. Increasing the magnitude of the discharge current results in relatively fewer background ions available to neutralize injected electrons and therefore a shallower density profile. The density at the center of the trap n_e(0) correlates linearly with discharge current magnitude.

Figure 9 demonstrates a correlation between spoke frequency and discharge current magnitude. Additionally, the average electron cross-field conductivity is computed at r = R₀/2 via Ohm's law [Eq. (1)] and plotted with respect to the discharge current magnitude in Fig. 10.

Increasing the applied magnetic field strength reduces electron mobility and therefore results in an enhanced ambipolar electric field. Figure 11 shows a near linear correlation between the radial electric field and the applied magnetic field strength. Therefore, as per Eq. (4), the spoke frequency should scale as $fs∼B$⁠. This is demonstrated in Fig. 12 and consistent with experimental observations.²⁸ 

Ion-mass m_i is varied and Fig. 13 shows a near linear correlation between spoke frequency and $1/mi$ (normalized to the ion mass of Helium-4), consistent with experimental observation.²⁸ This also demonstrates that the spoke rotation velocity scales in an identical way to the ion-acoustic velocity with respect to ion-mass.

For each of these parameter scans, the approximate theoretical estimate for the spoke frequency shows good agreement in both magnitude and scaling to the measured numerical spoke frequency, providing strong evidence for this scaling relationship and that the CSHI is the driving instability for the rotating spoke. It should be taken into consideration that the observed structure is clearly highly non-linear and turbulent, making it surprising that a simple estimate based on linear theory provides such an accurate fit.

A highly non-linear turbulent structure rotating in the azimuthal E × B direction is observed to form within full-size two-dimensional kinetic simulations of a Penning discharge. This structure exhibits a characteristic behaviour very similar to that of the rotating spoke observed in experiments. Electron cross-field transport within the discharge is highly anomalous, with a majority of the current being channeled through the spoke structure. Similarity between simulations with and without ionization suggests that the spoke is not the result of an ionization wave. The generated ambipolar electric field and the density gradient are aligned so as to destabilize the collisionless Simon-Hoh instability. The correlation of the resulting spoke frequency with the average radial electric field, gradient length scale, and ion mass supports the claim that the collisionless Simon-Hoh instability is responsible for its formation.

The physics work was supported by the Air Force Office of Scientific Research.

The code development was supported by the Princeton University Program in Plasma Science and Technology.

This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

This research used resources of the Perseus cluster at the TIGRESS high performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering and the Princeton University Office of Information Technology's High Performance Research Computing Center.