Physics of E × B discharges relevant to plasma propulsion and similar technologies

Igor D. Kaganovich¹, Yevgeny Raitses¹, and Andrei Smolyakov²

¹Princeton Plasma Physics Laboratory, Princeton NJ 08543, USA

²University of Saskatchewan, 116 Science Place, Saskatoon, SK S7N 5E2, Canada

This perspective describes joint efforts by the low-temperature plasma community and the plasma propulsion community to develop a rigorous understanding of the rich phenomena observed in E  ×  B devices, as summarized in presentations given at the E  ×  B workshop 2019, Princeton.¹ 

In E × B configurations, the external magnetic field is applied perpendicular to the discharge current in order to generate a strong electric field. Such devices are widely used for electric propulsion; the most common are the Hall or Hall-effect thruster (HET) and DC magnetrons for plasma processing. Some of the discussed effects are also relevant to radio frequency (RF) plasma sources, e.g., helicon sources and thrusters which make use a magnetic nozzle (MN) to generate thrust. It is widely recognized that many of the complex physics effects observed in these devices are not well understood. Here, we share perspectives from 25 leading experts within the field, report on the state-of-the-art in eight major subtopics, and propose recommendations on directions for future research.

Like many spacecraft components, it would be highly desirable to design plasma thrusters via robust engineering models, which through step-by-step procedures can deliver predictable performance and lifetime. The status quo is far from this ideal situation. The design and development of plasma thrusters are not fully based on first principles physics models but rather rely on a semi-empirical approach, combined with long and expensive lifetime tests. This is because plasma propulsion is one of few remaining applications where engineering developments are seriously constrained due to the lack of a complete description of plasma properties. For example, in Hall thrusters (HTs), one of the most critical challenges is to self-consistently simulate the electron cross field current, which affects the efficiency of propellant ionization and ion acceleration. Electron scattering in turbulent processes also affects the plasma–wall interactions and erosion because it changes the sheath potential at the walls. Sixty years of experimental and theoretical research on Hall thrusters has provided much insight into the instability driven mechanisms which contribute to the enhanced electron cross field transport in these plasma devices. Nevertheless, we still do not have a sufficiently clear understanding of these mechanisms and thus cannot accurately describe their contributions to the cross field transport. Consequently, reliable and robust methods to predictively design these thrusters do not yet exist. Such methods are especially important for new applications, e.g., very low power thrusters for CubeSats or extremely high-power thrusters for deep space exploration. For these applications, purely engineering solutions, such as size-scaling or clustering of multiple thrusters, are not necessarily optimal for achieving the required thruster characteristics. Any variation of established thruster design, for example, achieving thruster compactness or higher thrust density, developing operational envelopes through throttling and variable thrust, or altering the propellant from xenon, requires significant experimental work which could be greatly aided by physics-based predictive modeling.

To facilitate a discussion on these complex processes between experimentalists, numerical simulation experts, and theorists, a series of specialized workshops started in 2012.² A focused workshop on “E × B plasmas for space and industrial applications” was organized in Toulouse by Boeuf and Smolyakov in June 2017.³ A follow-up to the workshop included a collection of special topic papers titled “Modern issues and applications of E × B plasmas” published in Physics of Plasmas in 2018.³ The next workshop was organized in Princeton by Raitses et al. in October 2018.¹ The workshop participants discussed the following topics: “Validation & Verification for discharges and Sheaths,” “Mechanisms of Electron Cyclotron Drift Instability (ECDI) saturation and turbulence,” “Kinetic vs Fluid, Hybrid models,” “Low-frequency phenomena in E × B discharges, modeling and experiments,” “Experiments in turbulence,” “Towards full 3D modeling and GPUs,” “E × B Mass-filtering,” and “Unusual effects in magnetic nozzles.” Workshop participants—about 40 leading experts in the field of E × B plasma physics and devices—agreed that there is a need to prepare perspectives on each of these research directions.

Correspondingly, this perspective article discusses nine topics, which represent major directions for the electric propulsion community:

1. Plasma–wall interactions in E × B discharges relevant to propulsion plasma devices

2. Low-frequency oscillations in E × B discharges

3. Experiments in turbulence in low temperature, E × B devices

4. Electron drift instabilities in E × B plasmas: mechanisms, nonlinear Saturation, and turbulent transport

5. Fluid and hybrid (fluid–kinetic) modeling of E × B discharges

6. Toward full three-dimensional modeling of Hall thruster E × B discharges

7. E × B configurations for plasma mass separation applications

8. Validation and verification procedures for discharge modeling

9. Magnetic nozzles for electric propulsion.

Section II describes plasma–wall interactions. It is a well-established experimental fact for Hall thrusters that the wall material can affect the discharge current, the electric field, and the plasma plume. Thruster performance can, therefore, be affected by the wall material. This wall material effect is commonly attributed to the electron-induced secondary electron emission (SEE), which is different between materials. The secondary electron emission causes additional transport across the magnetic field, the so-called near-wall conductivity (NWC), that can increase the electron current across the magnetic field and lead to a reduction of the thrust. Sheaths near the wall determine the ion energy of ions impinging the walls and therefore, the sheath structure and sheath potential affect the wall erosion and, correspondingly, the thruster lifetime. The channel of the Hall thruster is narrow, and plasma is collisionless with the electron mean free path with respect to collisions with ions and neutrals being much larger than the channel width. Correspondingly, the sheath structure and the voltage drop on the sheath are determined by many kinetic processes that control electron fluxes to the walls. For example, the electron emission from the wall becomes especially strong if the secondary electron emission yield, which is the ratio of the flux of emitted electrons to the flux of primary electrons, reaches or becomes larger than unity. Under such conditions, intense electron flux can cause a nonmonotonic potential profile in the sheath, creating either so-called space charge limited or inverse sheath structures. The sheath can also become unstable in the presence of intense electron emission and instabilities can lead to oscillations of the wall potential, the so-called relaxation sheath oscillations. The value of the sheath potential is strongly affected by the non-Maxwellian character of the electron energy distribution function (EEDF) due to the depletion of the EEDF tail. Many processes form the EEDF tail. They include scattering in electron–neutral collisions as well as in turbulent processes, such as the oblique two-stream instability. Correspondingly, the wall fluxes and the sheath potential can be affected by anomalous transport. Examples of the complex interplay between all these processes are given in Sec. II.

Section III describes low-frequency (typically <100 kHz) oscillations, namely the breathing and spoke-type oscillations observed in thrusters. These constitute the most intense oscillations in E × B devices and may reach up to 100% modulation of plasma parameters. The breathing oscillations manifest as oscillations of the plasma and neutral densities that are coupled through ionization. These oscillations develop mostly along the direction of the external electric field and resemble predator–prey type oscillations. However, recent studies have shown that an adequate description of these modes cannot be reduced to a simple zero-dimensional description of predator–prey type oscillations. Theoretical models and simulations of the breathing oscillations suggest a large sensitivity to values of the effective electron mobility along the applied electric field (and hence sensitivity to a poorly known anomalous transport), which is still not well understood, and the details of the ionization zone near the anode. In contrast, the spoke develops mostly in the azimuthal or E × B direction, perpendicular to the crossed electric (E) and magnetic (B) fields. The spoke may exist without ionization and could be caused by the Simon–Hoh (SH)-type instability, driven by the combination of the applied electric field and gradient in the plasma density. Ionization can also strongly affect the spoke, especially for devices with higher pressures such as DC magnetrons. With a few exceptions, modeling of the breathing and azimuthal spoke oscillations has been performed separately and without considering possible coupling effects. Some experimental data suggest that there is a coupling between them which makes analysis complicated and will require 3D modeling for a deeper analysis.

Section IV describes experimental observations of high-frequency (>1 MHz) and short-wavelength (<1 mm) waves. These oscillations are believed to be the main contributing mechanism responsible for the enhancement of the electron cross field current in these devices. The electron anomalous cross field current is directly responsible for power losses in E × B devices that heat electrons and do not contribute to propulsion thrust. The wave characteristics are measured with laser-based diagnostics, coherent Thomson scattering (CTS), fast Langmuir, and emissive probes. It is very difficult to measure inside the narrow thruster channel and therefore most measurements are performed outside, where access is possible. Measurements indicate that high-frequency and short-wavelength waves are driven by the electron drift instability (EDI), which under some conditions exhibit cyclotron resonances (for which the frequency of modes is proportional to the electron cyclotron frequency). Under these circumstances, this instability is considered to be an Electron Cyclotron Drift Instability (ECDI). In the far plume area where the magnetic field is small, a wide spectrum of excited turbulent oscillations develops and has properties typical for ion–acoustic turbulence (IAT). The relationship between ECDI and IAT is still actively debated because large scale 3D simulations are required to understand the spectrum of resulting turbulence in steady-state.

Section V continues the description of high-frequency and short-wavelength waves, but unlike Sec. IV it provides a focus on the theoretical review of the current understanding of high-frequency and short-wavelength waves. In the case of a purely 2D system in the axial–azimuthal plane perpendicular to the magnetic field and with finite electron temperature, there exists a set of multiple narrow bands of unstable ECDI modes near the cyclotron resonances. The linear theory of ECDI instability is well developed. However, nonlinear saturation is still the subject of active research. Some earlier numerical work suggests that at a certain amplitude, the cyclotron resonances are washed out by nonlinear effects and the instability proceeds similar to the ion sound wave instability in unmagnetized plasma. Other studies indicate that the instability retains the cyclotron resonance features even in the nonlinear stage. When the direction along the magnetic field is considered, an additional instability appears due to the finite value of the wavevector along the magnetic field, either the so-called modified two-stream instability or the modified Buneman two-stream instability (MBTSI). This instability leads to the heating and scattering of electrons along the magnetic field. Given the complexity of mode saturation in realistic 3D geometry, only large-scale 3D simulations can provide a full understanding of the mechanism of anomalous transport and heating. This is made more challenging by the fact that several simulations indicate that secondary nonlinear processes take place resulting in the appearance of long-wavelength waves. The large-scale waves are typically expected to provide large contributions to the anomalous transport and need to be resolved for realistic parameters.

Section VI summarizes the current state-of-the-art in fluid (fluid electrons–fluid ions) and hybrid (fluid electron–kinetic ions) simulations of E × B plasma discharges. The fluid simulations are being used as a relatively nonexpensive alternative for full kinetic simulations with realistic device parameters, which are outside the current capabilities of existing kinetic codes. The first principles 2D nonlinear fluid models are based on rigorous closures for the electron dynamics perpendicular to the magnetic field (in the azimuthal-axial plane), which are based on the low-frequency approximation, i.e., the mode frequency is assumed to be low compared to the electron cyclotron frequency. Such models show that fluid-type instabilities, such as the gradient-drift modes of the Simon–Hoh type and their generalizations including the lower-hybrid instabilities, result in large anomalous transport generally consistent with experimental values. A notable effect demonstrated in such simulations is a nonlinear energy transfer from the most unstable small-scale modes (of the lower-hybrid type) to the large-scale structures of the simulation-box size. These models, however, are highly simplified; typically, they are in the azimuthal-axial plane and neglect the electron motion along the magnetic field; they also neglect sheath and plasma–wall interactions, and most importantly, consider ions in the fluid approximation. Such models are not suitable for the design of practical systems but could prove useful for understanding complex nonlinear processes in anomalous transport. Motivated by engineering design needs, for many years the alternative hybrid (fluid–kinetic) approach has been employed. This method describes the ions and neutrals kinetically, while electrons are treated with the fluid model. Such 2D, radial–axial models do not include the azimuthal direction; therefore, they require ad hoc or experimentally based models for the anomalous transport, typically parameterized by the anomalous collision frequency. The formulation and verification of the anomalous electron transport model remains a great challenge here. It is also unclear whether the concept of the anomalous collision frequency remains a good parameterization for the anomalous transport in the presence of large-scale dynamical structures such as spokes and breathing oscillations.

Section VII explains why a predictive model of E × B discharges requires a kinetic three-dimensional treatment. Typical descriptions of instabilities and transport in two dimensions are performed in the axial–azimuthal or axial–radial plane. In the first case, the important effects of plasma–wall interactions, most notably wall losses, are not captured and neither is the convective transport of energy from unstable regions to stable regions. In the second case, the EDI is not captured and only the breathing oscillations and other axial, e.g., gradient, modes are resolved. Recent simulations of the spoke in the anode region of the thruster show their intrinsic 3D structure. In general, reduced two-dimensional models always show stronger instability characterized by large amplitude oscillations. This leads to a significantly overestimated cross field mobility as compared with that simulated in three dimensions and observed in experimental measurements. Another finding from the recent comparison between 2D and 3D simulations is that the spectrum of excited waves in 3D does not exhibit strong peaks at a few dominant frequencies as observed in 2D. This observation is in accordance with both experimental and 3D linear analytic kinetic theory predictions that the instability is excited at a wide range of wavenumbers. For 3D simulations, future development of 3D codes should embrace modern computer algorithms; legacy codes generally do not scale well on modern computing architectures and therefore cannot perform full-size 3D simulations. Future codes should implement efficient data structures for memory access as well as hybrid parallelism via vectorization, OpenMP, or Message Passing Interface (MPI) in order to improve scalability up to millions of processors (exascale computing). Notwithstanding the difficulty of full 3D simulations, the complete understanding of electron transport will lead to a new era in the technological development of E × B plasma devices: designs based on an empirical approach will give way to code-based refined optimization. As it has been done in many other engineering disciplines, predictive design and optimization via computer-based techniques will assist and eventually replace more expensive empirical methods.

Section VIII describes recent progress in developing E × B configurations for plasma mass separation applications. Notwithstanding a long history of crossed-field (or E × B) configurations to separate charged particles based on mass, the limitations of small throughput of current devices drive innovation in this area and have led to the development of new plasma isotope separators, where crossed-field configurations were used to produce plasma rotation in plasma centrifuges. Applications of plasma-based elemental separation based on mass include nuclear waste clean-up, spent nuclear fuel (SNF) reprocessing, and rare earth element (REE) recycling. The separation in plasmas is conditioned upon the ability to externally apply a high electric field in the direction perpendicular to the magnetic field. This is limited by anomalous conductivity, similar to electric propulsion devices. Demonstrating the practicality of crossed-field mass filter concepts, therefore, hinges on a comprehensive understanding of anomalous perpendicular conductivity, which calls for combined modeling and experimental research efforts. Another outstanding issue in the presence of neutrals is the possible upper limit set on the rotation speed by the critical ionization velocity phenomenon. However, the promise plasma separation holds for many outstanding societal challenges is a compelling motivation to tackle these questions.

Section IX is devoted to verification and validation procedures. The ultimate objective for developing computer simulations of complex physical systems is to use these simulations as a predictive tool for science and engineering design. Critical to these goals is the need to verify and validate codes used for the predictive modeling of low-temperature plasmas. Though rigorous verification and validation (V&V) procedures have been relied upon for decades in other fields, it is only recently that these procedures were applied to simulations of E × B devices. The section describes V&V efforts in low-temperature plasmas and specifically for 2D axial–radial simulations of Hall thrusters. For validation, a comprehensive set of measurements is needed. This is challenging for compact and energetic Hall thrusters where probes can strongly perturb the plasma and it becomes difficult to measure plasma parameters in the ionization and acceleration zones. A possible approach is to validate codes on specially designed plasma systems which allow for improved access for diagnostics, similar to the Gaseous Electronic Conference (GEC) cell for RF discharges. An example of such a study for E × B discharges is the penning discharge. For Hall thrusters, perhaps a scaled-up device with improved diagnostic access is desired for accurate validation of the simulation results. Another approach is to use a wall-less Hall thruster where the acceleration zone is outside the thruster channel. Comprehensive measurements by several laser diagnostics and fast probes are needed for the complete characterization of anomalous transport and oscillation spectra. Currently, there is a need to develop verification practices for low-temperature magnetized plasmas. It is important to emphasize that the definition of the test cases is a significant part of the problem. Indeed, it is important to develop a comprehensive set of test cases to provide benchmarking of codes for the regimes of interest that sufficiently characterize the relevant physics important in E × B plasmas. For low-temperature magnetized plasmas, these are anomalous transport, low-frequency oscillation such as breathing oscillations, and plasma spokes. A detailed description of validation and benchmark test cases will be the subject of future dedicated E × B workshops and the “Frontiers in Low-temperature Plasma Simulations” workshop. These workshops will facilitate further discussions and ideas.

Section X discusses magnetic nozzles. The magnetic nozzles are used to control the radial expansion of the plasma jet, in a similar way as a de Laval nozzle operates on hot ideal gas. The section discusses which part of the electron velocity distribution function (EVDF) is responsible for transforming energy from electrons into ion kinetic energy. This process can be affected by collisions in the magnetic nozzle. This phenomenon is more acute inside a vacuum facility due to the additional effects of the background pressure and of the electrical connection between the plasma beam and the metallic chamber walls, which can affect the electric potential profile, the EVDF of confined electrons, and the amount of electron cooling. Kinetic models are necessary to rigorously describe the magnetic nozzles. However, efforts were directed to develop fluid models with closures that can be sufficient to predict thrust. Complex sets of magnetic coils can create three-dimensional magnetic structures with adaptable shapes capable of steering the plasma jet and allow for a nonmechanical way of thrust vector control. The effectiveness of these approaches for weakly magnetized ions remains unexplored. A variety of plasma instabilities can develop and potentially affect magnetic nozzle operation and the resulting thrust. This research field has been little explored so far, both theoretically and experimentally.

Last but not least, we assembled a very comprehensive list of references with the goal to provide readers with a detailed list of references that should be credited for the original work. We ask readers to cite the original papers where appropriate. In addition, individual chapters of the perspective could be mentioned in a similar way to individual chapters in a book.

Eduardo Ahedo¹, Michael Keidar², Irina Schweigert², Pascal Chabert³, Yevgeny Raitses, and Igor D. Kaganovich⁴

¹'Equipo de Propulsión Espacial y Plasmas, Universidad Carlos III de Madrid, Leganés 28911, Spain

²George Washington University, Washington D.C. 20052, USA

³Laboratoire de Physique des Plasmas, CNRS, Ecole Polytechnique, Sorbonne Universités, 91120 Palaiseau, France

⁴Princeton Plasma Physics Laboratory, Princeton NJ 08543, USA

Plasma–wall interactions and magnetic field effects have a strong influence on the plasma discharge operation for many plasma applications, from electric propulsion to magnetic confinement fusion, material processing, and plasma diagnostics. There is a large body of evidence that wall materials affect Hall-effect thruster (HET) operation. The early observations⁴ showed that “when a ceramic stick is introduced into the channel radially, the discharge current increases significantly. The discharge current, propulsion, and efficiency coefficient of the thruster are very sensitive to the insulator material and contamination of the surface.” Following the fundamental theoretical predictions of Morozov and Savelyev,⁴ systematic experimental studies of these effects began in 1978 in Bugrova's laboratory at MIREA.⁵ Direct demonstrations of wall material effects on thruster performances^6–8 have shown that it can affect strongly the discharge current, both in terms of its magnitude, Fig. 1, and its fluctuation frequency (known as the breathing mode), Fig. 2(a). To a lesser extent, it affects thrust and thus thrust efficiency (the ratio of the thrust squared to the input electric power times twice the mass flow rate), Fig. 2(b).

A detailed comparison of plasma properties and discharge operations for Stationary Plasma Thruster (SPT)-type thrusters with boron nitride and carbon segmented walls was performed at Princeton Plasma Physics laboratory (PPPL), see Ref. 9 and references within. Figure 3 shows the V–I characteristics, the maximum electron temperature, and the maximum electric field measured for the boron nitride ceramic channel and the channel with the nonemitting carbon–velvet walls; a big effect of wall material on these properties is evident for voltages above 350–400 V. The effect was attributed to the electron-induced secondary electron emission (SEE) for the boron nitride surface as compared to the carbon walls and a corresponding increase in the near-wall conductivity leading to the increased electron conductivity in the channel for boron nitride ceramic walls. Detailed analysis of the discharge (and shown electron temperature) is complicated by the fact that the electric field is affected by the emission and strongly changes with the wall material as described in Ref. 9, shown in Fig. 3(c). However, the relationship between the emission properties and the discharge current is not always straightforward. For example, in Ref. 10, a floating graphite electrode and passive electrode configurations were implemented with two ceramic spacers made from a quartz and MACOR (machinable glass–ceramic). The measurements for the discharge currents show the following trends: for the MACOR spacer 1.63  A, for the graphite electrode 1.66  A, for the boron–nitride channel 1.7  A, and for a quartz spacer 1.94  A. Note that the secondary electron emission from quartz is lower than that for MACOR and boron nitride.¹⁰ The PPPL study⁹ showed that metal or ceramic spacers can significantly affect the plasma potential distribution and shape of equipotential lines in Hall thrusters, even though these spaces have sizes comparable but smaller than the acceleration region. A theoretical study of these effects was performed in Refs. 12 and 13. Mechanisms behind these effects possibly involve changes in the electron temperature and transport, which are affected by the secondary electron emission and conductivity of wall materials. The effects of spacers on operating conditions of the thruster depend on the precise placement of the electrodes relative to the magnetic field distribution.¹⁰ Another potentially important effect is the contamination of the surface; for example, quartz is not a porous material whereas boron nitride is, therefore that latter is more prone to contamination. Reference 11 showed that contamination can greatly affect SEE.

Similar effects of secondary electron emission were observed in magnetrons, another E × B discharge. The effects are due to ion-induced rather than electron-induced SEE another EB discharge. Reference 14 discusses the influence of the secondary electron yield on the property of the spokes in High Power Impulse Magnetron Sputtering (HiPIMS) plasma, by adding nitrogen and effectively changing the target surface from a metallic to compound target, affecting the secondary electron yield of the target surface.

All these experimental observations stimulated comprehensive theoretical and modeling studies aimed at the description of complex plasma processes in HETs as well as accurate measurements of secondary electron emission from dielectrics.^15–18 It became clear that the theories have to be kinetic, as fluid theories cannot explain such high electron temperatures observed in HETs due to fast losses to the channel walls, see Fig. 3(b). Such high electron temperatures were also recently confirmed by measurements with incoherent Thomson scattering (ITS) diagnostics¹⁹ in addition to the previous measurements by probes.⁹ 

In typical electric propulsion devices and conventional annular Hall-effect thrusters, the magnetic lines intersect the walls and therefore plasma confinement is not magnetic but mainly electrostatic (ES); walls are typically biased negatively with respect to plasma bulk in order to collect the correct electric current (for instance, zero in the case of a dielectric wall). The absence of magnetic confinement in a conventional HET leads to high particle and energy losses at the walls and surface erosion. Furthermore, the large secondary electron emission (SEE) induced by energetic electron impact on ceramic walls reduces the wall potential, therefore worsening the electrostatic confinement of energetic electrons, and increasing the energy loss from the plasma.

The mean free path for electron collisions with neutrals or Coulomb collisions is large compared to the HET channel width. Consequently, losses of energetic electrons are expected to be higher than losses in inelastic collisions with neutrals for the conditions of conventional HETs.^8,9,20,23 Moreover, the high-energy electrons that can escape to the walls with energy above the potential energy corresponding to the wall potential are not being replenished sufficiently fast. As a result, the electron velocity distribution function (EVDF) is non-Maxwellian, as shown in Fig. 4. The electron flux to the walls is proportional to the EVDF for these energetic electrons escaping to the wall; whereas the EVDF in this region is small compared to the Maxwellian EVDF. The uncertainty in the EVDF makes it difficult to model HETs using fluid or hybrid approaches which are based on the assumption of a Maxwellian EVDF.²¹ 

Indeed, Particle-In-Cell (PIC) simulation studies of the EVDF formation in HETs are showing a rich plasma material interaction phenomena that is non-amenable to simple scaling laws of wide use. These studies include all necessary effects in self-consistent treatment: collisional replenishment of the EVDF for energetic electrons, the influence of the SEE, and the Debye sheath formation. By contrast, scaling laws were derived for sheath potential and near-wall conductivity for a case of a magnetic field perpendicular to the wall.²³ 

It was also shown using PIC simulations that the average energy (or effective temperature) of the electrons decays in the sheath due to the non-Maxwellian distribution functions mentioned above. Therefore, the usual isothermal sheath theories used to evaluate the sheath potential drop and the secondary emission yield are not correct.^24,25 A polytropic sheath model has been proposed to overcome the limitation of isothermal models. The model works well but is not self-consistent and uses the PIC data to evaluate the polytropic index.

Two-stream instabilities of the SEE beams injected into the plasma can reduce the energy of the SEE due to exchange with colder bulk electrons and increase the number of energetic SEE electrons trapped by a potential well and therefore make the total EVDF closer to a Maxwellian.²⁷ The evolution of SEE within the plasma needs to take into account possible re-collection by the walls.^26,28 Asymmetries in the EVDF caused by cylindrical effects can be significant too.^29,64 The sheath structure and plasma–wall interaction processes are quite sensitive to the detailed description of SEE features (e.g., to the amount of true-secondary vs elastically or inelastically backscattered electrons).^27,30

At SEE yields of about 100%, the sheath becomes space-charge limited, and in this regime, the sheath may become unstable.^31–34 Furthermore, in the center of the acceleration region of a HET, the axial electric field is maximum and the resulting E × B drifts can be on the order of the electron thermal velocity, thus enhancing the impact energy of electrons and the transition to a space-charge limited unstable sheath.³⁰ 

The sheath instability (in the radial direction) may work as a trigger for the azimuthal fluctuations. An example of the temporal evolution of the total current collected on the outer wall and the corresponding floating potential is shown in Fig. 5. It should be pointed out that in this case the instability is detected only on the outer wall, where the secondary electron emission coefficient reaches a value larger than 1.

In the case of high SEE yield, the transition was observed from a space-charge limited sheath to an inverse sheath.³⁶ The important difference between the two regimes is that ions are not accelerated toward the wall in the inverse sheath case and thus improve ion confinement (but deteriorate energy losses). This phenomenon has been studied in the context of large thermionic emission, mainly related to emissive probes and arcs.^37–39 Excitation of ion–acoustic waves in the presence of a very intense SEE when an inverse sheath forms was observed in the HET channel for a high electric field.⁴³ 

Sheath oscillations were experimentally observed when a beam of electrons strikes a biased plate.⁴¹ The oscillating sheath near the floating emissive plate bombarded by a beam of energetic electrons from the cathode was also observed in 2D particle-in-cell (PIC) Monte Carlo Collisions (MCC) simulations in Ref. 42 [see Fig. 6(a)]. The virtual cathode appears, as expected, when the total electron flux from the plasma produces a larger flux of secondary electrons from the emissive surface. The potential profiles for different times of sheath oscillation cycles shown in Fig. 6(a) are related to the periodical accumulation of secondary electrons near the emissive surface. A fragment of oscillating floating potential is shown in Fig. 6(b). The numbers in Fig. 6(b) point out the time of snapshots of the potential and electron density profiles shown in Figs. 6(a) and 6(c). The sheath oscillation frequency of about 25 kHz is set by the rate of production of secondary electrons and the ion velocity.

One of the important aspects of HET application is the lifetime which is determined primarily by ion-induced erosion of the thruster channel walls.^4,44,45 It is important to note that wall erosion is quite visible (few mm in depth) after a few hundred hours of thruster operation in SPT-100 type thrusters,^46–48 leading to the exposure of magnetic circuits and eventually ill-functioning. The rate of erosion is governed by the ion flux to the wall and the sputtering yield of the material, which depends on the ion energy and the incidence angle.^4,44 An example is shown in Fig. 7(a), where wall erosion causes significant changes to the channel shape. Due to a big variation in the channel wall shape the erosion rate decreases, see Fig. 7(b). In addition to the above volumetric erosion, there exists abnormal erosion forming striations [shown in Fig. 8(b)] on the material surface which could be related to plasma or sheath oscillations but is largely unknown.⁴ Due to the high economical value of both increasing and qualifying lifetimes, simulation tools for lifetime characterization and extension are an important research effort.

Most models of plasma in HETs consider a simple configuration of a magnetic field that is normal to the channel walls. If the magnetic field is oblique to the wall, the characterization of the partially depleted tail of the EVDF becomes more complex and the E × B drift can be altered significantly.^4,40 Consequently, wall erosion might be mitigated by applying such an oblique magnetic field. This is particularly relevant for innovative designs of HETs bearing magnetically shielded topologies, (i.e., with magnetic lines near-parallel to the walls around the thruster exit). The magnetic-shielded topologies appeared after the discovery that the erosion of the BN walls in the BPT-4000 Hall thruster designed and developed by Aerojet had essentially stopped after 5600 h of operation during life testing;⁵¹ a similar observation was made in the Soviet Union,^52,53 see Fig. 7. A follow-up extensive investigation into the mechanism responsible for this effect was initiated at the NASA Jet Propulsion Laboratory (JPL). The magnetic shielding reduces, of course, the plasma ion fluxes to the walls, but, indirectly, the mean impacting ion energy as well, because the plasma acceleration region moves downstream.^49,50 Then, there is the THALES High-Efficiency Multi-Stage Plasma Thruster (HEMPT), a propulsion technology rather similar in operational principles to the HET where magnetic shielding is applied through magnet-based cusped magnetic fields.⁵⁴ Different conceptions of the cylindrical HET lie between HET and HEMPT designs.⁵⁵ Much remains to be known on characterizing the local E × B drifts, the EVDF, and plasma fluxes to walls in these complex topologies. As discussed above, there were many efforts on optimization of magnetic topology to better focus the plasma and minimize plasma–wall interactions, resulting in the development of high-performance Hall thrusters of SPT and Thruster with Anode Layer (TAL) types. However, the magnetic shielded thruster developed by JPL and Aerojet is currently a culmination of these efforts in terms of the extended lifetime.⁵⁶ 

Interesting effects are associated with the formation of filaments in a magnetic field under conditions similar to those in HETs.^57,58 The plasma–wall interactions in the external oblique magnetic field were studied in cases of different magnetic field strengths and incidence angles for the conditions similar to the Hall thruster but for a simpler plasma source. To this end, PIC MCC simulations were employed.^59,60 The results shown in Fig. 8 illustrate the electric potential distribution and alignment of electron density peaks along the magnetic field for different angles. The magnetic striation due to the instability was theoretically predicted in Refs. 61 and 62 and experimentally observed in Ref. 63. These striations can possibly explain striations in the erosion pattern observed at the exit of the thruster, shown in Fig. 8(b) from Ref. 66.

Secondary electron emission at the wall is predicted to be a key phenomenon in E × B discharges, particularly in Hall thrusters. However, there is still no direct evidence that it is the SEE alone that is responsible for the wall material effect. It will be important to validate modeling predictions by a comparison with measurements of electron velocity distribution function or energy distribution function (EDF) in the thruster channel using, for example, Laser Thomson Scattering (LTS) and/or electrostatic Langmuir probes, respectively. Implementation of both diagnostics inside the thruster channel where the predicted effects are important is challenging. For example, for the LTS, key challenges are (1) to get the laser beam and scattered light from the channel inside and (2) to resolve the tail of the EVDF. For probes, probe-induced perturbations of the plasma are hard to avoid unless the probe is installed on the channel wall that may change the wall effect.

For modeling, in order to appropriately model the SEE effects, the challenge is to correctly evaluate the EVDF. This is naturally done in PIC simulations, but unfortunately, these simulations are computationally expensive and therefore not yet adapted for thruster design optimization. It is, therefore, necessary to propose simplified sheath models, which account for non-Maxwellian EVDFs and would allow better evaluation of the secondary emission yield. This is a very challenging task and a large research effort is needed in this area.

As it was shown above, the oblique magnetic field adds considerable complexity in the form of magnetic striations. How the magnetic striations form in the HET geometry is not yet explored.

The biggest complexity is due to the mutual interaction of turbulence and plasma–wall effects because turbulence leads to electron heating and scattering and can populate depleted parts of the EVDF due to loss cone. These phenomena, though observed in PIC simulations, are not thoroughly studied to the degree that they can be accounted for in simplified sheath models as discussed above.

Finally, it would be important to study changes in the material properties due to exposure to the thruster plasma and how that change affects the thruster operation.

From the technological point of view, two possible directions could be followed in parallel. The first is to continue the research efforts in materials science to improve the performances of dielectric walls, both in terms of intrinsic properties (for example, by reducing secondary electron emission) while providing sufficiently long lifetime and aging properties. The progress in this direction has been already done and research in this area probably needs to be maintained. The second avenue is to propose new designs that would minimize the role of plasma–surface interactions. This is already happening with the development of magnetic shielded or wall-less thrusters. Other ideas that would allow controlling the sheath properties could be proposed.

From the scientific point of view, combined efforts are needed in theory, numerical simulations and experiments to validate predictions of theory and simulations. In theory, new sheaths models that account for non-Maxwellian EVDF are needed to be able to predict the total secondary electron emission yield at the surface.

On the simulation side, a large community effort is needed to develop massively parallelized and optimized codes able to simulate turbulence and plasma–wall interaction in 3D geometries. These codes will also have to be benchmarked and verified.

Finally, optical (noninvasive) diagnostics, able to achieve both short time and space resolution, are required to investigate the plasma–wall interaction effects in these high-density and high-energy plasmas where probes hardly survive. This is difficult and is a major challenge for the community. Moreover, it is also important to continue the development of electrostatic diagnostics suitable for measurements of EEDF in harsh plasma environments of the thrusters. As in high-temperature fusion research which pioneered and uses many sophisticated RF, optical and laser diagnostics of plasmas, probes and energy analyzers remain important diagnostics for edge physics and especially so at the divertors of fusion reactors.

There was a large progress in understanding plasma–wall interaction in the recent decade. A future challenge is to generalize the developed understanding for realistic devices in 3D accounting for self-consistent interaction of turbulence and plasma–wall interactions and provide experimental validation of these effects using advanced diagnostics.

Yevgeny Raitses¹, Andrei Smolyakov², Mark Cappelli³, Kentaro Hara³, Jean-Pierre Boeuf⁴, and Igor Kaganovich¹

¹Princeton Plasma Physics Laboratory, Princeton NJ 08543, USA

²University of Saskatchewan, 116 Science Place, Saskatoon, SK S7N 5E2, Canada

³Stanford University, Stanford, California 94305-3032, USA

⁴LAPLACE, University of Toulouse, CNRS, INPT, UPS, 118 Route de Narbonne, 31062 Toulouse, France

Low-frequency (typically <100 kHz) oscillations, including spokes and breathing modes, are some of the most prominent examples of self-organization in partially magnetized quasineutral plasmas of crossed-field devices at a wide range of pressures ∼0.1–100 mTorr, such as Hall thrusters, Penning discharges, and sputtering magnetrons.^7,65,67–70 The spoke mode manifests itself as strong perturbations in plasma density that propagate in the E  ×  B direction perpendicular to the crossed electric (E) and magnetic (B) fields, generating substantial components in electric field in this E × B direction.^65,67–69 The breathing mode propagates in the direction of the external electric field, and is one of the powerful modes observed in Hall thrusters. It reveals itself in oscillations of the discharge current [see Figs. 9 and 2(a)], often reaching amplitudes comparable to the mean discharge current itself.^7,65,71 While these modes have been known and studied for some time, only recently have coordinated efforts been made toward self-consistent modeling,^72–77 detailed experimental validation,^78–82 and their control.^83,86–88

For cylindrical devices, the spokes rotate azimuthally in the direction of the E  ×  B drift, but with a speed that is an order of magnitude smaller than the E  ×  B drift velocity.^{67,68,78,80,89–92} The mode number of these spokes is usually low, m = 1 − 8.^{67,68,81,93,94} It has been reported that anomalous (turbulent) electron current may be enhanced in the spoke, carrying 20%–90% of the total discharge current.^91,92 Although the mechanism for spoke formation is still debated, one candidate is the Simon–Hoh (SH)-type instability,^95,96 driven by the combination of the applied electric field and the gradient in plasma density. A modified theory of this instability for partially magnetized collisionless plasmas was developed^{74,79,97,98,107} and experimentally verified for some conditions.^{67,79,92,93,97}

Recent results of large-scale Particle-In-Cell (PIC) simulations of a Penning discharge^75,108 were found to be in good agreement with experimental data,⁹² showing the formation of a m = 1 spoke rotating with a frequency of a few kHz (Fig. 10), generating anomalous current due to fluctuations. The scaling of the spoke frequency deduced from both the simulations and experiments was consistent with theoretical predictions for SH instability in collisionless plasmas,⁹⁷ i.e., $f∝ E r L n / m i$⁠, where E_r is the radial electric field, L_n is the gradient density scale, and m_i is the ion mass.⁷⁵ According to these PIC simulations, the ionization of the working gas has a minor effect on the spoke formation.

Unfavorable magnetic field gradients can lead to violent perturbations, but typical magnetic profiles in modern Hall thrusters, are designed to partially stabilize the most disruptive spoke modes.⁶⁵ The theory of collisionless SH instability has been modified to include the gradients of the magnetic field and electron temperature. With these modifications it is generally referred to as the gradient-drift instability.⁷⁴ Linear analysis^109,110 as well as recent particle-fluid hybrid simulations show that gradient-drift effects are critical for the formation of azimuthally rotating spoke-like structures in the typical conditions of Hall thrusters.⁸¹ 

Although simulations have shown that spokes can form in collisionless plasmas without ionization,⁷⁵ there is strong evidence from numerous experiments in high power pulsed magnetron discharges (High Power Impulse Magnetron Sputtering—HiPIMS discharges) and in lower power, DC magnetron discharges^99,100,103 that ionization can play an important role in the formation and dynamics of spokes. An important similarity between magnetron and Hall thruster discharges is the existence of a region where the radial component of the magnetic field decreases axially toward the anode. Recent experiments in HiPIMS and DC magnetron discharges have shown that the observed spokes are associated with an ionization instability rotating in the azimuthal (E × B) direction. The measurements reveal the existence of a double-layer structure with a large electric field at the leading edge of the ionization zone and this double layer plays a crucial role in the energization of electrons.¹⁰³ References 101 and 102 discuss the experimentally determined potential structures of spokes in HiPIMS. The presence of a double layer and enhanced ionization at the spoke front is reminiscent of the critical ionization velocity model of Ref. 105 and of the PIC simulations of Ref. 63 and 106 although the velocity of the spoke measured in HiPIMS and magnetron discharges does not necessarily match the theoretical critical ionization velocity. Moreover experiments,^{79,99,100,103} and recent modeling¹⁰⁴ have evidenced that spokes can rotate in the −E × B direction (at low power) as well as in the +E × B direction (at higher power).

The very detailed and recently published measurements of the spoke properties in HiPIMS and DC magnetron discharges should open the way for a better understanding of these structures and of their dynamics based on modeling and simulations. Although the physics of HiPIMS is more complex (e.g., ionization of the sputtered atoms), it seems very likely that spokes in Hall thruster and magnetron discharges share several common properties.

Ionization plays a key role in the axial (breathing) mode oscillations.^{70,81,82,111–113} Experiments demonstrated the dependence of these oscillations on the thruster operating conditions such as discharge voltage, magnetic field strength and topology, gas flow, cathode operation, and background pressure.^7,81,82,86 A zero-dimensional predator–prey model for the coupled evolution of the ion and neutral density¹¹² predicts the oscillations, however, the equivalent one-dimensional model shows no oscillations. It was also realized that self-consistent electron dynamics is important.^69,113 The resistive instability^114,115 due to the phase shift in the response of the electron current (resistive) and ion current (inertial) to the perturbation of the electric field was proposed to be a triggering element of the breathing mode.^81,116,117 Besides, the electron temperature evolution was also shown to be important and needs to be included into fluid or hybrid models.^69,81,114 Recent studies that include fluctuations of the electron temperature have shown that temperature may render zero-dimensional predator–prey models unstable.⁸¹ Without the variation of the electron temperature and ionization rate taken into account, the one-dimensional modeling shows no excitation of the breathing oscillations.^117,118

There is a growing realization that in conditions relevant to the practical operation of E × B plasma devices, the low-frequency azimuthal and axial modes result from a complex interplay of various phenomena. An understanding of these complex interactions remains far from complete and they present a critical challenge for the development of predictive modeling tools needed for existing and future applications. For stronger magnetic fields and faster ion and electron rotations, it is important to investigate the onset of a Simon–Hoh instability taking into account centrifugal forces^119,120 that are especially relevant to novel mass separation devices discussed in Sec. VII.

As for experiments, it is challenging to independently control the discharge parameters such as local electric field and density gradients—the properties important to the instability. This coupling makes it difficult to identify the cause for the ubiquitous presence of large-scale azimuthal disturbances while the smaller scale azimuthal modes typically have larger growth rates: for the low-m modes the mode growth rate increases almost linearly with the wavenumber.^67,97 Fluid simulations have shown that the low-m modes can be formed as a result of nonlinear inverse energy cascade from small-scale modes.^74,107 Furthermore, magnetic field, temperature, and plasma gradients, which are present in many engineering devices (e.g., Hall thrusters and magnetrons),^65,77,121 and the presence of physical boundaries and associated sheaths⁹² all need to be taken into account.¹²² Combined effects of gradients, nonlocal electron kinetics, and the coupling between large and small-scale plasma structures (energy cascade) should continue to be explored theoretically, numerically and experimentally. It is also important to further explore the effects of ionization and neutral depletion in such complex, highly nonuniform plasma systems.^123–125 

Theoretical models and simulations of the breathing mode suggest a large sensitivity to the values of the effective electron mobility along the applied electric field or even wall material as was shown in Fig. 2(a). Empirical values of the anomalous mobility are typically used for modeling of the breathing oscillations with values adjusted to achieve reasonable agreement with experiments.¹²¹ In addition to anomalous mobility, the models should include other anomalous transport coefficients, including anomalous heating and energy losses. For the most part, physical mechanisms governing these anomalous effects are not well-understood. For example, it is unknown whether the anomalous heating can be described by the same effective collision frequency as the anomalous electron mobility. Understanding these anomalous effects is required for further progress in this research field.

With a few exceptions, modeling of the breathing and azimuthal spoke modes has been performed separately and without considering possible coupling effects. Some experimental data suggest that there is a coupling between them.¹²⁶ A coupling between the two modes is expected because the azimuthal mode is driven by axial gradients in plasma parameters such as the density and electric field, which experience large spatial and temporal variations during breathing mode oscillations. Control of the axial oscillations modifies the driving forces for the azimuthal modes.¹²⁷ Recent control experiments have used varying cathode electron emission⁸⁶ as well as external voltage modulation,¹²⁷ demonstrating a suppression of the oscillations associated with both of these modes.

Modern computational tools to study E × B plasma devices include fluid, hybrid (PIC for ions and atoms and fluid for electrons), and full PIC (PIC for all plasma species) simulation codes. Benchmarking should be performed for simulations with identical conditions, specifically relevant to low-frequency phenomena. Codes should be benchmarked against each other to better understand the limitations of numerical algorithms and physics approximations. For instance, a fluid model with electron pressure closure¹¹¹ has been proposed to eliminate the numerical uncertainties in conventional quasineutral models. Alternative approaches such as continuum grid-based kinetic models, or direct kinetic simulations,^128,129 can be used to understand the issue of numerical noise in PIC codes. PIC codes can be used to verify the numerical diffusion issue in direct kinetic simulations.

Crucial for validating simulation results are experiments capable of measuring plasma properties with spatial and temporal resolution, including energy distribution functions (EDFs) of the electron, ion, and neutral species in different directions along with the electric and magnetic fields in the E × B direction, and time-resolved (<10⁻⁴ s) electric field measurements. Experiments involving active control of low-frequency oscillations^86,87,111 are also important as they may reveal underlying physical mechanisms of instabilities needed to be captured by truly predictive models.

The use of fast-sweeping electrostatic probes and energy analyzers for measurements of electron and ion EDFs, and plasma potential is appropriate when these invasive diagnostics induce minor plasma perturbations (e.g., Penning systems). However, for Hall thrusters and magnetron discharges, with nonuniform magnetic fields and strong potential gradients, nonintrusive diagnostics, such as Laser-Induced Fluorescence (LIF) of electronically excited ions and atoms^130–135 and Laser Thomson Scattering (LTS),¹³⁶ may be more quantitative in measuring species EDFs in these devices. A critical challenge for time-resolving LIF of the ion EDF is accounting for electronically excited states produced during spoke and breathing oscillations by direct ionization of neutral atoms as well as ions in other electronically excited states.^137,138 For measurements of the electron EDF using LTS, a key challenge is a relatively low electron density detection limit of 10¹⁰ cm⁻³, making it difficult to characterize the EDF high energy tail that may develop and be affected in low frequency spoke and breathing mode cycles.

Low-frequency oscillations occurring in E × B discharges are the most powerful oscillations and therefore, may significantly affect the performance of these devices in several different ways including but not limited to power losses on electron transport and heating, defocusing of ions, and mismatch between the device and the power supply. Therefore, it is important to understand these oscillations and their control. The progress made in understanding the low-frequency phenomena, including spoke and breathing oscillations, has been driven by combined efforts of modern experimentation, theory, and simulation studies. However, there is an emerging need for dedicated efforts to compare and benchmark various numerical codes. Experiments capable of measuring spatially and temporally resolved plasma properties during low-frequency oscillations will be crucial to advance our understanding of these phenomena. Continued advances in computational capabilities will eventually allow simulations to be carried out at scale and over the times needed to resolve these low-frequency structures. The ultimate challenge and the long-term goal would be the development of experimentally validated predictive computational tools that self-consistently model anomalous electron mobility and heat conduction, because these transport phenomena play a critical role in low-frequency oscillations as observed during nominal operation of E × B plasma devices. Then, the next step would require the implementation of modeling including electric circuits with passive and active control of these oscillations and other means of their control such as segmented electrodes,⁸³ cathode gas flow,^84,85 and electron emission.⁹⁰ 

Benjamin Jorns¹ and Sedina Tsikata²

¹Department of Aerospace Engineering, University of Michigan, Ann Arbor, MI, USA

²ICARE, Electric Propulsion Team, Centre National de la Recherche Scientifique, Orléans, France

While there are many types of plasma oscillations in low-temperature, E  $×$ B devices, including the longer-wavelength modes linked to plasma inhomogeneities and ionization (see Secs. III and VI),^70,140 there is a growing interest in experimentally characterizing the role of short wavelength (<1 mm) plasma turbulence in these devices. This interest largely has been motivated by numerical studies (starting with the work of Ref. 142 and more recently, in work such as Ref. 141 and references therein) that have shown these oscillations may be a dominant driver of anomalous electron transport in low-temperature E  $×$ B systems. To this end, experiments have investigated which waves are present in these devices, how much energy is contained in these waves, and how these measured wave properties can be related to anomalous transport. These experimental studies have been facilitated by the development of new diagnostic techniques in the past fifteen years—such as coherent Thomson scattering (CTS), developed in response to numerical work of Ref. 142—that have allowed for unprecedented levels of noninvasive access and resolution.

Most experimental efforts on turbulent measurements to date have focused on fluctuations formed in the E  $×$ B direction (although other short wavelength fluctuations—the axially directed two-stream instability— have been detected and characterized as well¹³⁹). For example, Fig. 11 shows a set of representative experimental results for a Hall thruster discharge. In these experiments, which are described in Refs. 143 and 147, CTS was employed to interrogate the short-wavelength oscillations in the electron density in the direction of the Hall drift. Figure 11(a) illustrates the geometry of the system and the measurement location, downstream of the thruster exit, with investigations performed at an observation wave vector k with variable orientation in the (E  $×$ B, E) plane. Figure 11(b) shows the resulting measured wave dispersion from CTS in the Hall (azimuthal) direction. It was proposed in Refs. 143 and 147 that the features of the measured dispersion—specifically, of MHz, mm-scale density fluctuations in the azimuthal direction—were indicative of the presence of the electron cyclotron drift instability (ECDI). The characteristic physics of the ECDI (also referred to more generally as the electron drift instability, or EDI) are discussed in more detail in Sec. V.

To further support the link between experiment and ECDI, Refs. 143 and 147 drew parallels between their measurements and previous 2D numerical and linear kinetic theory studies^142,154,182 that had predicted that the ECDI would propagate in the Hall thruster plasma. Subsequent measurements on other Hall thruster configurations using different probing techniques¹⁴⁴ as well as studies on a similar E  $×$ B device, the planar magnetron operating in the pulsed, high-current, high-density regime known as HiPIMS¹⁵⁵ also have shown dispersions which the authors have attributed to the presence of ECDI. The growing evidence that the ECDI exists in these devices is physically impactful as this mode was the first one linked (in simulations¹⁴²) to turbulence-driven anomalous cross field transport.

With that said, as discussed in more detail in Sec. V, there is some debate in the community about the appropriateness of labeling the measured waves as so-called ECDI/EDI. This stems from the fact that the linear dispersion in Fig. 11(b) exhibits features that could also be described simply as ion sound (group velocities in the Hall direction commensurate with the local ion sound speed). This is an intriguing result as the classical ion sound is unmagnetized and therefore would not be expected to appear in the Hall direction. To reconcile this apparent incongruity, some theories (Sec. V) posit that when wave amplitudes become sufficiently large for the ECDI, the waves become effectively unmagnetized, leading to ion sound-like behavior. Alternatively, experimental measurements have motivated a different explanation. In particular, measurements have shown the existence of a finite wave component along the magnetic field.¹⁵⁷ It subsequently has been demonstrated (Refs. 156 and 158 and in Sec. V) that the presence of 3D effects like this parallel propagation also may account for the sound-like dispersion.

Another question that experiments have tried to address is how much energy is contained in the excited waves. This is a critical consideration as the anomalous transport depends both on the wavelengths of the excited modes and the energy in each mode (Sec. V). Early simulation work based on 1D and 2D kinetic theory indicated that energy could be concentrated at the electron cyclotron resonance corresponding to maximum growth (e.g., see Fig. 14 in Sec. V) with large amplitudes (exceeding >10% of the thermal background). This type of energy concentration is manifested in numerical results as a well-defined, characteristic length scale for the instability (e.g., Fig. 16 in Sec. V). The first experimental measurements, in contrast, indicated that the shape and amplitude of the power spectrum measured in actual thrusters is significantly different.

Figure 11(c) illustrates this contrast with a plot of the distribution of energy over wavelengths from the CTS measurements performed on the same device where the dispersion was measured in Fig. 11(b). This representation shows the energy spectrum scaling [using the calibrated density fluctuation amplitude, known as the dynamic form factor S(k, ω)] over the range of wavenumbers k and frequencies ω found in experiments in Ref. 143. Such information has been used to establish scaling laws (see Ref. 143) relating the integrated fluctuation intensity S(k) to k. An exponential decay with wavenumber has been observed, with an e-decrement on the order of the electron Larmor radius (<1 mm), and amplitudes < 1% of the thermal background. The shape of this spectrum is consistent with a classically turbulent distribution in which the growth of the energy of the waves is governed by nonlinear processes, and significantly, wave energy is not limited to a single wavenumber. This departure from the predictions of early simulation results in part motivated more recent numerical studies. These new simulations have led to an updated hypothesis for the energy content of the ECDI where it has been proposed that the energy spectrum may be subject to an inverse energy cascade.¹⁶³ 

In an effort to better elucidate the processes that could lead to a nonlinear state, Ref. 144 explored the spatial evolution of the energy spectrum of oscillations in the Hall direction. Since the Hall thruster plasma convects at high speed and the ECDI moves with the plasma, it was argued in this work that measurements of the energy spectra at different distances from the thruster might show the transition of the ECDI from its initial growth upstream in the plasma to a nonlinear state. To this end, a set of translating ion saturation probes were used to measure the energy spectra of ion density fluctuations in a 9-kW class Hall thruster as a function of normalized axial distance (z/L) measured from the thruster anode.

The results are shown in Fig. 12(a), where for reference, the energy spectrum of the discharge current oscillation, I_D, is also plotted. At low frequencies (∼10 kHz), there is a peak in the energy spectra that is consistent with the so-called breathing mode (Sec. III). At higher values of frequency, there are number of marked features that are qualitatively consistent with the nonlinear growth processes outlined Ref. 163. Closer to the thruster, there are a series of discrete peaks (denoted with dotted vertical dashed lines). These were shown in Ref. 144 to correspond to the resonant frequencies of maximum growth predicted from linear theory (Fig. 14 in Sec. V). With increasing distance from the thruster exit plane, the amplitude of these peaks (1 MHz–20 MHz) decreases. Concurrently, there is an increase in energy from 100 kHz to 500 kHz where the spectrum is broadband, showing an inverse decay in frequency. The exchange in energy from discrete structures at high frequency to these broadband lower frequencies is captured by Fig. 12(b), which shows the percentage of total wave energy in each frequency range as a function of position.

The transition of energy from the high frequency peaks near the thruster to a more broadband, lower frequency spectrum is consistent with the interpretation that the ECDI waves are born in the upstream region and grow as they are convected downstream. The growth of the spectrum transitions from linear to nonlinear, giving rise to an inverse energy cascade to lower frequency (longer length scales). With that said, we note that although the data shown in Fig. 12 is consistent with the idea of an inverse energy cascade, experimental measurements of this transition to date have been confined to frequency and not wavelength (a limitation of the measurement technique). To fully explore the mechanisms leading to the nonlinear energy spectrum, future measurements will need to measure the full range of wavelengths.

In addition to studying the dispersion and energy content of the drift-driven waves, experimental efforts to date also have focused on trying to link the measured wave properties to anomalous transport known to exist in these devices. To this end, most experimental methods have relied on a quasilinear approximation (cf. Ref. 150), which can be leveraged to relate transport coefficients such as an anomalous collision frequency to wave amplitude and growth rate. Ideally, this method requires simultaneous measurements of the electric field and density fluctuations associated with the waves. In practice, this type of measurement has not been feasible in these devices, and instead, linear approximations are employed to relate what can be measured (typically density) to the potential fluctuations in the waves.^151,152 For example, in low temperature electrostatic modes, it is common to use the Boltzmann relationship.^148,149,158

Subject to these simplifying assumptions about transport and the wave properties, experimental measurements have been used to calculate effective collision frequencies for electrons in the plumes of E  $×$ B devices. It was shown in Ref. 149, for example, that depending on the assumptions about the properties of the ECDI, the electron transport from waves can account for 25%–100% of the measured electron transport in a Hall thruster. Similarly, an attempt was made in Ref. 158 to determine the electric field directly from CTS by using calculations on the electron density fluctuation to determine a corresponding fluctuating potential. This gave an ECDI field amplitude that was 25% of the background electric field in the far-field thruster plume. 2D numerical simulations performed in Refs. 144 and 142 demonstrated that fields this large could account for all of the anomalous electron transport.

With that said, there are a number of reservations about the validity of this preliminary, quasilinear analysis. For example, the evident turbulent spectrum of the fluctuations already suggests that nonlinear effects cannot be ignored. Similarly, the simplifying assumption that the electron density fluctuations can be related to the potential fluctuations through a Boltzmann relationship is not consistently supported by findings from available PIC simulations. The use of the Boltzmann relationship also relies on the assumption of thermal electrons, which does not apply to much of the region where the instability is generated. Recent PIC work (Sec. V) also shows that the applicability of quasilinear analysis requires further study. Linking experimental measurements of turbulence to anomalous transport, in light of the complex nature of the particle properties and behavior in the plasma region of interest, therefore poses a considerable challenge and remains an active area of research.

While the exact relationship between turbulence measurements and transport remains an open question, experimental measurements of plasma turbulence are informing a new, correlational understanding of how turbulence-induced transport may impact self-organized behavior in thrusters and magnetrons (Sec. III). For example, experimental measurements have indicated that the ECDI can coexist with lower-frequency (kHz) breathing-mode oscillations,^145,146 and experimental studies on magnetrons have revealed the presence of the ECDI and its coexistence with kHz-frequency rotating oscillations or spokes discussed in Sec. III.

An example of such a result obtained in recent years from CTS experiments on pulsed planar magnetrons is shown in Fig. 13, from Ref. 147. In Fig. 13(a), this image shows the averaged symmetric MHz frequency spectrum characteristic of the azimuthally propagating ECDI measured at Larmor-radius observation scales. A time-resolved analysis of the ECDI over a single 150 μs operating pulse, shown in Fig. 13(b), illustrates how these MHz fluctuations exhibit an additional kHz range modulation (measurable due to the fluctuations in density as spokes traverse the measurement region). Recent experimental work on hollow cathodes operating in conjunction with a Hall thruster also showed that not only does plasma turbulence dominate the electron transport in the cathode, but it also may be linked to lower-frequency oscillations.^149,150 Taken together, these experimental results from cathodes, Hall thrusters, and magnetrons now point to the conclusion that turbulence-driven transport has a critical role in governing both the steady-state plasma properties as well as the transient, large scale oscillations discussed in Sec. III.

The most pressing need for future experimental efforts is to clarify the relationship between different instabilities and anomalous transport. Although the ECDI is now considered a key factor in this process, supported by simulations showing its appearance leading directly to an axial electron drift (cf. Ref. 154), it is unlikely to be the only instability concerned. The balance of contributions from different unstable modes may depend on the nature of the operating regime. However, evaluating the contribution of each mode to transport poses a number of experimental and theoretical challenges.

One of the most significant challenges is estimating the electron transport directly from plasma wave properties measured experimentally. As discussed in the preceding, the nonlinear state of the measured instabilities calls into question¹⁴⁶ the current practice of using simplified expressions for transport coefficients derived from quasilinear theory. The least ambiguous approach to overcome these challenges is to measure the amplitude and relative phase of electric field and density fluctuations simultaneously across the energy spectrum. Ideally, future diagnostics would offer this capability.

This need also could be addressed by finding methods or validated theories to map more precisely the fluctuating electric field amplitude to density. While the very nature of the phenomena under study (weak fluctuations, high-frequency modes) renders this a significant challenge from a theoretical perspective, detailed simulations may offer a potential solution. For example, one possible future technical path could involve developing higher fidelity simulations, validating the predictions of these simulations against measurements that can be made with current techniques (e.g., CTS or probe-based measurements of density fluctuations), and then using the simulation results to determine the appropriate map from density to electric field for the experiments. To be effective, this type of approach would need to address the role of significant nonlinear effects (such as the coupling of wall emission to the ECDI,¹⁴⁶ which affects the overall electron current) and 3D effects present in experimental devices that are often absent from the linear kinetic theory and certain model descriptions.

Another key future challenge for experimental measurements will be to document the growth and energy flow of the excited instabilities. As discussed in the preceding and Sec. V, the earliest notions regarding ECDI-driven transport as an electron Larmor-scale phenomenon are now being revised in light of observations in simulations and experiments of the development of large (cm-scale) modes, which have been linked to an energy transfer from smaller scales. Future efforts will need to assess this type of energy flow experimentally. This will require improvements in both the spatial and wavelength resolution of current diagnostics.

Building on the work in Ref. 144, characterizing the spatial evolution of the turbulence as it convects out of the Hall thruster can provide critical clues about the growth and saturation of the turbulence. However, this type of spatial interrogation poses a problem for current, noninvasive optical techniques that cannot reach the internal parts of the discharge chamber in thrusters. This latter issue regarding optical diagnostic access may potentially be alleviated through the study of modified architectures, such as so-called magnetically shielded thrusters, which shift much of the acceleration and ionization region beyond the confines of the thruster channel. However, such architectures induce modifications to the plasma which further complicate physical interpretation.

Future efforts to study energy flow also will require an expanded capability to interrogate a wider range of length scales (1 mm < $λ$ < 1 cm) than currently can be accessed with state-of-the-art methods. These measurements similarly should have the ability to characterize all three dimensions of propagation of the waves in order to assess the spatial direction of energy flow. Given that there is evidence of coupling between the turbulence and low-frequency, self-organized structures, these future diagnostics also should have the capability to measure turbulence properties on the timescale of the lower frequency modes discussed in Sec. III (i.e., ∼100 kHz).

As the sophistication of diagnostics continues to grow, the insight that emerges from experiments must continue be leveraged to inform improved, predictive kinetic modeling. The comparison between experiments (measuring 3D phenomena) and numerical simulation results can only truly be considered valid if the simulations can eventually reach a sufficient level of sophistication, accounting for (i) the time scales needed for the study of fast and slow waves, (ii) the full range of spatial scales (from electron Larmor radius-scales to centimetric scales of large-scale fluid turbulence), and (iii) the 3D plasma environment itself. An example of the importance of this convergence is illustrated in the following way. Numerical work from a large number of groups (involving 1D, 2D PIC codes) uniformly show two persistent features for the ECDI-like fluctuations developing in the region of the E  $×$ B drift: (i) the establishment of a single dominant length scale for the mode (<1 mm), and (ii) the alignment of these fluctuations primarily along the azimuthal direction (but with a positive axial wave vector component). In contrast, more recent 3D simulations (see Sec. VII) show clearly (i) the weakening of the coherent nature of the ECDI fluctuations (which become more compatible with a broadband structure, already established in CTS experiments) and (ii) the disappearance of a preferred propagation direction outwards. The differences already observable between 1D, 2D, and 3D simulation results are strong evidence of the need for 3D simulation codes which can provide a more solid basis for future numerical/experimental comparisons.

Beyond kinetic simulations, experimental measurements also will be critical for guiding fluid-based attempts to approximate the effects of the turbulence on background plasma properties (Sec. VI). Indeed, while fluid simulations cannot capture the kinetic effects of the waves directly, these processes can be represented with approximate closure models. The fidelity of these closure models in turn depends directly on an understanding of the evolution and growth of macroscopic properties of the turbulence such as the total wave energy. Recent work on the study of anomalous electron transport in the hollow cathode, the electron source for Hall thrusters, could serve as a potential roadmap for future efforts.

In these previous works, experimental measurements of the evolution of a turbulence and anomalous transport in both space and time were leveraged to guide the development of fluid codes that included equations to approximate the effects of plasma turbulence on electron resistivity.^{84,149,150,153,164,165} These models in turn have been successfully used to simulate the plasma in both the interior and exterior of the cathode. A similar approach could be applied for modeling the effects of turbulence-induced anomalous transport in E  $×$ B devices, though, unlike in the cathode studies where the turbulence lent itself to a simple 1D description, closure models for the effects of E  $×$ B modes likely will be more nuanced. Before it is possible to develop approximations for these effects in a fluid framework, many unresolved questions related to the turbulence discussed above (3D propagation, inverse energy cascades, etc.) will need to be resolved.

Many of these technical challenges outlined in Sec. IV B may be addressed by building on existing diagnostics techniques. For example, the current implementation of CTS was designed for the purpose of accessing shorter wavelengths (<2 mm) as this is where PIC simulations suggested the ECDI should exist. In principle, the CTS technique can be modified to access a wider range of wavelengths and therefore capture the relevant content at longer wavelengths. This increase in wavelength range comes at the expense of spatial resolution (as there is an approximate inverse relationship between measurable wavelength and spatial resolution). Future efforts should focus on expanding the capabilities of CTS while maintaining high spatial resolution. With that said, assuming CTS still will have a practical lower bound in wavelength, probes can be used to fill in information about the energy spectra of oscillations at longer length scales.

The recent development of a sensitive incoherent Thomson scattering (ITS) platform¹³⁶ allowing investigations of cathode,¹⁵⁹ thruster,¹⁶⁰ and magnetron plasmas¹⁶¹ gives access to information that has long been lacking regarding background electron properties (temperature and drift) in these sources. This provides crucial information that can impact the predicted growth and density fluctuations in the measured waves. It is hoped that the coupling of improved CTS and ITS can be an important tool for validating basic physical understanding regarding the conditions in which certain instabilities arise. An understanding of the velocity distribution of electrons also may help lead to refined expressions for relating measured wave properties to transport.¹⁷⁶ Similarly, there have been recent advances in techniques based on laser induced fluorescence that have allowed, for the first time, the noninvasive measurement of anomalous electron transport inside Hall thrusters.¹⁶⁶ These data will be crucial for providing “ground truth” for the local transport that can be used as a point of comparison to determine if the contributions from instabilities, if any, are sufficient to explain the electron dynamics.

As discussed in the preceding, in order to directly estimate transport from the measured instabilities, it is necessary to make measurements of the electric field fluctuations and electron density fluctuations simultaneously. Both CTS and probe-based measurements that have been employed to date, however, have focused on characterizing oscillations in the plasma density. In principle, probes could be used to measure both potential and density (with limited spatial resolution), though probes are known to perturb the local plasma state. Alternatively, there are diagnostic techniques from the study of higher energy plasmas such as Stark broadening that may offer a potential technical path for characterizing the electric field fluctuations noninvasively and at small wavelengths. There are several technical hurdles stemming from accessibility and signal to noise ratio that must first be overcome, however, before these methods can be applied to low temperature systems.

Numerical simulations (and, as described above, 3D simulations in particular,) may offer a key capability for bypassing the need to make direct measurements of electric field. Before these can be used, these simulations must first be validated against measurements of both the background and wave properties. With this in mind, establishing the error bars of invasive and noninvasive techniques will be important for future progress in the use of models to inform experiment. It will be necessary to explain and quantify differences observed between probes and optical techniques for the measurement of electron properties, and to establish on what data future simulation efforts should be validated, given the limitations and challenges of different diagnostics. As an example, measurements of the electron energy distribution function using ITS rely on obtaining a sufficiently high signal-to-noise ratio, which is difficult for plasma densities in the range of 10¹⁶ m⁻³, even now with recent diagnostics improvements. This requirement would ultimately affect the detection of features pertaining to the spectral wings such as high-energy electron populations.

Combining information from complementary diagnostics—on electron density fluctuations, EEDFs, absolute electron density, electron temperature, and ion properties—is required to refine models and theory, and gradual progress is being made on this front. The combination of multiple diagnostics running simultaneously on the same test device would provide an ensemble of useful data for simulations. As noted above, a concerted attempt to bridge the gap between measurements at large scales using probes, and small-scale measurements using CTS, could be an important contribution toward establishing the presence (or absence) of dominant wavenumbers in anomalous transport. As these experimental capabilities become available, it will be necessary to analyze experimental findings coupled with sufficiently advanced numerical codes, still under development, which are capable of capturing 3D plasma features.

There has been substantial progress in the past decade in experimental methods and diagnostics that have yielded new insights into the role of plasma turbulence in E  $×$ B devices. The detection and in-depth study of different types of instabilities present in such devices, such as the ECDI, ion acoustic turbulence, spokes and high-frequency, long-wavelength modes^170,171 have contributed to our understanding of basic physics. Several challenges remain, including reconciling experimental data with simulation predictions, understanding the limitations of invasive and noninvasive diagnostics, and studying the growth and nonlinear evolution of the observed turbulence and anomalous transport in E  $×$ B devices. Efforts in diagnostic development and analysis are under way to help address these challenges. Looking to the future, the development of standardized test devices, particularly ones that can be shared internationally, could enable coordinated and more impactful efforts that leverage the diagnostic capabilities of multiple institutions. Similarly, closer collaborations between experimentalists and modelers will be critical for establishing the validity of codes currently under development, understanding the role played by different instabilities in thruster operation, and guiding the design of future E  $×$ B devices and modeling efforts to study plasma turbulence. It has become evident that understanding features of plasma turbulence measurable experimentally, and the contribution of different instabilities to transport, will also ultimately require 3D numerical code results for comparison. Parallel progress in both diagnostic development and numerical simulation capability is a clear objective for the future.

Andrei Smolyakov¹, Jean-Pierre Boeuf², and Trevor Lafleur³

¹University of Saskatchewan, 116 Science Place, Saskatoon, SK S7N 5E2, Canada

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

³PlasmaPotential, Entry 29, 5/1 Moore Street, Canberra ACT 2601, Australia

Electron transport across the magnetic field in devices employing magnetic filter configurations is typically anomalous and exceeds the transport due to classical collisions of electrons with other particle species in the plasma. Typically, the inclusion of wall-collisions (near-wall conductivity) is not sufficient to explain experimental values inside the channel. Moreover, the electron current significantly exceeds the classical value outside the channel, where the near-wall conductivity is absent. Thus, it is widely believed that the anomalous transport enhancement is due to convective transport and scattering of electrons due to turbulent plasma fluctuations. The nature of such fluctuations is neither well understood, nor does there exist any validated theoretical model that can be used to predict the turbulent electron transport for specific thruster conditions.

The electron drift instability (EDI) due to the electron drift has attracted a lot of attention recently as a possible instability and anomalous transport mechanism in Hall thrusters.^{141,142,172–176} Such interest was also stimulated by experimental observations of small scale fluctuations where the wave-frequency is found to be linear with the wave-vector, which is consistent with the ion sound wave dispersion relationship, see Fig. 11.

Studies of this instability were started much earlier in relation to the problem of anomalous resistivity to explain the width of collisionless shock waves in space and turbulent heating experiments.^177–180 The instability occurs as a result of the differential drift of electrons with respect to unmagnetized ions, when the Doppler frequency shift results in an overlap of the kinetic resonances of electron cyclotron (Bernstein) type modes with the ion sound branch. As a kinetic instability, it does not require plasma or magnetic field gradients, and is purely based on the electron E × B drift in crossed magnetic and electric fields. When applied to Hall thrusters, the EDI is of great interest in the region of large electric field, where one can expect that the kinetic instability due to this large electric field will dominate over other instability mechanisms, such as plasma gradients and collisions.⁷⁴ 

The genesis of EDI mechanisms can be tracked down to the magnetized Buneman instability in cold plasmas.¹⁸¹ This is the reactive instability which occurs as a result of the interaction of two stable modes: the upper hybrid resonance, and the low-frequency ion oscillations, $ω 2 = ω p i 2$ (which is the short wavelength limit of the ion sound mode). The upper hybrid mode is Doppler shifted by the electron $v E = E × B / B 2$ drift so it moves into the ion plasma $( ω p i )$ frequency range. In the plasma with a finite electron temperature, the overlapping modes are the Bernstein and ion sound modes, $ω c e − k y v E ≃ ω$⁠, so that the approximate resonance condition is $k y v E ≃ ω c e$⁠, where $k y$ is the wavevector in the direction of the $E × B$ drift, which corresponds to the azimuthal direction in Hall thruster and magnetron geometries, and $ω c e$ is the electron cyclotron frequency. When one allows electron motion along the magnetic field, an additional instability appears due to the finite value of the wavevector along the magnetic field, $k z ≠ 0$⁠; the so-called Modified Two-Stream Instability. This instability was also considered in the original paper,¹⁸¹ so below it will be called the modified Buneman two-stream instability (MBTSI), which has the largest growth rate for relatively small values of $k z ≪ k y ≃ ω c e / v E$⁠.

For purely perpendicular propagation and finite electron temperature, there a exists a set of multiple narrow bands of reactively unstable modes near the resonances $ω− k y v E−m ω c e=0$⁠, $m=1, 2,….$⁠. When electron parallel motion is included, for finite $k z$⁠, kinetic (Landau) resonances become possible leading to dissipative instabilities at $ω− k y v E 0− k z v | |−m ω c e≃0$⁠. For $k z v T e ≥ ω c e ,$ the kinetic resonances broaden and eventually overlap, resulting in an instability equivalent to the ion sound instability in an unmagnetized plasma, but driven instead by the electron E × B beam flow,^156,174,182 see Fig. 14. The various limits of this linear instability (usually referred to as the electron cyclotron drift instability or ECDI) are well studied; however, nonlinear behavior and saturation are much less understood, especially for small or zero $k z$ when linear Landau interactions due to particle motion along the magnetic field are not possible. Some earlier studies¹⁷⁹ have concluded that at a certain amplitude, the cyclotron resonances are washed out by nonlinear effects and the instability proceeds as an ordinary ion sound wave instability as found in unmagnetized plasmas until it is saturated by linear or nonlinear Landau damping by ions. Other studies have argued^183,184 that electron trapping in the magnetic field remains important, making the nonlinear stage different from that of an unmagnetized plasma ion sound instability. Alternative mechanisms (different from cyclotron and kinetic resonances) have also been proposed¹⁸⁵ to explain the nonlinear instability; hence, below, we refer to the general nonlinear regime of this instability as the electron drift instability (EDI).

The ECDI/EDI is a very robust instability and can easily be seen even in the simplest versions of 1D PIC simulations of magnetized plasmas, e.g., see the cold plasma example in Ref. 186. In the context of the anomalous transport in Hall thrusters, it has been studied in 1D simulations,^141,163,174 2D axial–azimuthal simulations,^{142,173,176,187} and 2D radial–azimuthal simulations.^146,188,190 Typically, in such simulations, significant anomalous transport in the axial direction as well as very effective electron heating are observed. In 1D versions, the instability occurs for perturbations propagating strictly in the direction of the electron drift and typically shows up in the simulations as a fairly coherent nonquasineutral mode with strongly peaked ion density perturbations which are noticeably larger in amplitude than those of the electron density which remains relatively smooth,¹⁴¹ see Fig. 15.

Quasilinear models of electron transport due to the EDI have been proposed^{175,193–195} based on the assumption that the underlying turbulence is ion-acoustic as is in the case without the magnetic field. In the quasilinear approximation, the magnitude of the turbulent flux was estimated as $Γ =$ $⟨ n ̃ E ̃ y ⟩ / B 0$⁠, using relationships between the density and electric field fluctuations from the linear theory of ion–sound instability in unmagnetized plasma. The wave amplitude at saturation in these models,^{175,193–195} was estimated from the linear wave kinetic equation for the evolution of the wave energy $W k = | E | k 2$ convected axially by the ion flow (with the wave group velocity in the axial direction) and a combination of the linear Landau damping from warm ions and nonlinear ion trapping. However, since the wave kinetic equation is linear in the wave energy $W k ,$ the result is sensitive to the choice of the initial “reference” value at some given location. A further refinement was done in Ref. 193, based on the comparison with PIC simulations which show substantial non-Maxwellian distribution functions for the electrons. Using non-Maxwellian distribution functions¹⁹³ reduces the ion-sound growth rate, bringing the estimate of the fluctuation level down and in closer agreement with magnitudes expected experimentally.

Experimentally, it was found that fluctuations near the acceleration region exhibit discrete cyclotron resonances characteristic of the ECDI while downstream the turbulence evolves into the linear dispersion relation typical of the long-wavelength ion sound modes.^144,196

There are significant challenges in numerical simulations of the electron drift instability as well as in the theoretical understanding and interpretations of the results of such simulations. Because of computer resource limitations, many simulations have to be performed for parameters far realistic: reduced dimensionality, low resolution, and limited spatial coverage. Several model approximations for ionization sources and energy particle losses are also often used. The EDI has proved to be a very effective mechanism of electron heating, leading to a fast rise in the electron temperature. In the 1D periodic domain (and in some 2D) simulations, a “virtual” axial acceleration length and/or inelastic collision losses have to be introduced to achieve a stationary temperature state. The level of anomalous transport may become dependent on the mechanism of how the electrons are “cooled down” in these simulations. Concerns have been raised¹⁶⁸ that the strong electron heating may be in part due to numerical instabilities. Many simulations are performed over a limited azimuthal region (e.g., in a domain around 1 cm wide vs 30 cm in reality). Such narrow widths may prevent the development of long-wavelength azimuthal modes that typically dominate the anomalous transport. Investigations of the effects of larger azimuthal simulation boxes and extending them to realistic values are difficult because of computer resource limitations, but will be required to become relevant to experiments. The inclusion of elastic and inelastic collisions, ionization, as well as full 3D effects make such simulations very resource demanding.

Among basic questions are the effects of the magnetic field on the electrons, in particular, whether the description in terms of fully unmagnetized ion sound is appropriate for modeling the EDI in the nonlinear regime and on the saturation mechanisms. The differences between the ion sound type turbulence in unmagnetized plasmas and in a plasma with a magnetic field have been much debated previously^183,198,199 without definitive conclusions that can be directly applied to Hall thruster parameters.

In the linear case, the transition to the ion sound occurs for finite values of the wave vector along the magnetic field with $k z v T e ≃ ω c e$⁠. Several PIC simulations were performed however, in 1D (azimuthal) and 2D (azimuthal–axial) geometry, where the direction along the magnetic field is ignored, $k z = 0$⁠, and thus the linear mechanism of the transition to the ion sound is absent. The cyclotron resonances can be smoothed out by collisions for $( ν / Ω c e ) k 2 ρ e 2 > π / 2$⁠, forcing the instability into the ion–sound regime. However, for most typical regimes of interest, simulations were collisionless or almost collisionless. For the most part, numerical noise, which can play the role of collisions, was estimated to be in the range of $ν n = 5 × ( 10 − 5 − 10 − 6 )$ $ω p e$⁠;^197,198 therefore, according to the criteria above, it is not expected that numerical noise can result in electron demagnetization. The statistical fluctuations $∼ 1 / N$ however can be essential in simulations with low numbers of particles per cell, $N$⁠, and can affect the level of anomalous transport.

Another reason for smoothing out the cyclotron resonances is nonlinear resonance broadening²⁰⁰ which can effectively demagnetize electrons. For short-wavelength regimes, with $k ρ e > 1$⁠, a simple estimate for the electron demagnetization has the form $Ξ > ( k ρ e ) − 1$⁠,¹⁶³ where $Ξ ≡ ( ω p e 2 / Ω c e 2 ) W / ( n 0 T e )$⁠, $W = E 2 / 8 π$⁠. It has been reported that in some 1D and 2D simulations^163,190 this criterion is not satisfied (though not with a large margin). Consistent with this, the nonlinear evolution of the electron cyclotron instability in these simulations shows distinct electron cyclotron resonances at $k y v E=m ω c e,$^{163,190–192} as well as the ion sound modes features, e.g., mode propagation with a phase velocity of the order of the ion sound speed and oscillations (and nonlinear harmonics) at $ω = ω p i$ in the short wavelength part of the spectrum.

Azimuthal-radial simulations, in which the $k z$ (in the radial direction) is finite, and the linear mechanism to the transition to ion sound should remain operative, also showing inconclusive results. Quasicoherent nonlinear waves at the cyclotron resonances moving with roughly the ion sound velocity have been observed in highly resolved simulations,¹⁹⁰ similar to the 1D case, with strongly nonlinear cnoidal waves peaked at short wavelengths. The important role of the modified two-stream instability, due to a finite $k z$⁠, leading to strongly anisotropic heating and large scale radial structures in the anomalous current, was emphasized.¹⁹⁰ Various values of the effective $k z$ were reported in different simulations.^24,193–195 Effects of Secondary Electron Emission (SEE) and different propellants were studied^24,189 and it was shown that the EDI activity is affected by the reduced electron temperature due to sheath cooling. Besides, SEE induces large electron transport due to near-wall effects.

The 2D azimuthal–axial simulations, while remaining relatively simple, are the closest to real Hall thruster configurations (in some ways) and allow one to test the effects of the axial profile of the electric field (the main EDI driver) and the effects of axial mode propagation. In these simulations, the direction along the magnetic field is ignored and the linear mechanism of the transition to the ion sound does not apply. Nevertheless, these simulations reveal strongly coherent nonlinear waves propagating with a velocity of the order of the ion sound velocity, a wavelength that scales with the Debye length, and a wave amplitude that seems to be consistent with estimates from the ion trapping mechanism.¹⁸⁷ The azimuthal–axial simulations in Refs. 193 and 176 report transport levels that are generally consistent with the model in Refs. 174 and 201 and some experimental results.

It appears that in the nonlinear stage, the EDI typically reveals itself as a coherent strongly nonlinear wave, akin to the periodic cnoidal wave with characteristic features of wave breaking manifested by sharp peaks in the ion density and with smaller and much smoother electron density perturbations. This is consistent with the regime of unmagnetized ions, which therefore exhibit the tendency of wave breaking (similar to the neutral gas sound wave modes), while the electrons remain largely magnetized and show a smoother density. The coherent, almost single mode, emerging in many numerical simulations, e.g., see Fig. 16, presents a challenge to the application of quasilinear theory, which assumes an ensemble of weakly nonlinear, wide spectrum overlapping modes. The anomalous transport estimate in the form $Γ =$ $⟨ n ̃ E ̃ y ⟩ / B 0$ assumes that the dominant electron current is due to the E × B drift of magnetized electrons, which appears to be inconsistent with the regime of unmagnetized electrons for the ion sound instability in the absence of a magnetic field, although simulations^187,193 seem to show that this estimate is roughly valid.

For typical parameters in the acceleration region, the electron drift velocity $v E = E 0 / B 0$ becomes comparable to the electron thermal velocity (and much larger than the ion sound velocity, $v E≃ v T e≫ c s$⁠). In such regimes, nonlinear electron and ion dynamics become strongly nonlinear²⁰² and numerical modeling is critical for theoretical advancement and validation. This regime is far from the marginal stability criteria and normally a broad spectrum of excited modes is expected. Nevertheless, many PIC simulations demonstrate the excitation of a highly coherent mode and very effective electron heating. The mechanism of anomalous transport and heating in the electron interaction with quasicoherent waves and the role of numerical noise have to be understood. Several simulations indicate that secondary nonlinear processes take place resulting in the appearance of long-wavelength modes (similar to modulational instabilities) on top of the quasicoherent mode. The large scale modes are typically expected to provide large contributions to the anomalous transport and need to be resolved for realistic parameters. Therefore, simulations with larger spatial boxes representative of realistic geometrical dimensions are important and will require substantial computer resources.

High-performance large scale simulations need to explore how the current results (with low resolution and narrow simulation boxes) can be extrapolated to real-sized devices. Many reported simulations from different groups have been performed with varying approximations and assumptions, which makes the comparison more difficult and inconclusive. Recently, a broad collaborative effort between several groups has been initiated to investigate the accuracy and convergence of the results of different numerical codes^203,204 under the same conditions. A 2D axial–azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas instabilities has been recently completed.²⁰⁴ The results obtained for this benchmark show good agreement between the different codes from several groups.

One has to note also that the ECDI is only one of a number of instabilities that could be relevant to Hall thrusters, and more generally, to E × B devices; e.g., resistive axial and azimuthal instabilities can be important as well as azimuthal Simon–Hoh type modes^74,75,107 that are described in Secs. III and VI.

Andrei Smolyakov¹, and Ioannis G. Mikellides²

¹University of Saskatchewan, 116 Science Place, Saskatoon, SK S7N 5E2, Canada

²Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, 91109, USA

A major challenge in the modeling of E × B discharges such as those in the HET has been that the governing processes are inherently three-dimensional (3D) and span multiple scales. Specifically, the spatiotemporal resolution must span the device length down to the electron Larmor radius (ρ_e) as well as the long timescales associated with the motion of atoms down to those for the electron motion. Moreover, it has been argued that the species velocity distribution function can be far from Maxwellian, implying that the nature of plasma instabilities in these discharges may be strongly affected by kinetic effects. If true, fully kinetic simulations would be the most comprehensive approach to model such conditions. However, in many cases, such simulations are as complex as experiments and can be difficult to interpret. Though petascale computing is now possible, kinetic simulations with realistic parameters that span the entire HET domain remain computationally very expensive. Although particle methods for all three species in z–θ^142,173,201 and r–θ^30,146,206 domains, and even in 3D,^168,169 have indeed been used successfully and have yielded critical insight in some cases, they are constrained to restricted spatial and/or temporal domains, which limits our ability to capture the entire range of physical processes. In HETs for example, though on one hand simulations in the z–θ plane can provide a detailed insight into instabilities with wave vectors in the E × B and axial directions, they provide no information about plasma–wall interactions which are known to affect the operation of the device. On the other hand, in the r–θ domain the simulations cannot take into account the axial variation of the plasma properties and therefore cannot provide insight into how the transport physics evolves throughout the different regions of the channel.

Fluid theory predicts several instabilities that may be responsible for turbulence and transport in E × B discharges. In general, fluid simulations that account for plasma turbulence and thus produce self-consistent transport are computationally faster and cheaper compared to kinetic simulations. Such simulations are also easier to interpret, provide much greater flexibility in separating various physics elements, and are vital in developing intuition of the complex processes that occur in plasma discharges. Over the years, fluid models have been providing indispensable contributions to the development of nonlinear physics of hot and dense magnetized plasmas such as in space and laboratory devices for fusion applications.^222,226,231

Generally, in rarefied plasmas where the particle collision mean free paths are large, the fluid models suffer from the fundamental problem of closure. When the magnetic field is strong enough such that the wavelength (λ) of the modes of interest and other characteristic length scales (L) are large compared to ρ_e, (ρ_e/L, ρ_e/λ) ≪ 1, the dynamics in the plane perpendicular to the magnetic field is well represented by the fluid model. In fact, reduced fluid models taking into account the higher-order gyroviscosity tensor can capture well the behavior of the plasma even in regimes of short wavelengths (ρ_e > λ).²³³ Many phenomena in E × B discharges occur with a characteristic frequency well below ω_ce. Therefore, another restriction of fluid models, namely, the low-frequency approximation, ω ≪ ω_ce, is typically well satisfied for many processes of interest.

Specifically, for E × B discharges, fluid models and insights from the theory of fluid instabilities have been used successfully to characterize a wide range of fluctuations, especially in HETs.²³⁴ For example, using the proper magnetic field profiling the stabilization of the most violent large-scale azimuthal modes driven by plasma density and magnetic field gradients^68,97 have been predicted.^{65,224,228,229} The generalization of the fluid theory to shorter length scales (of the order of ρ_e and smaller) leads to the lower-hybrid modes which are destabilized in E × B plasmas by density gradients and collisions.¹¹⁸ The gradient-drift instabilities (i.e., modes driven by density gradients and electron E × B drift) have been invoked to explain turbulence in magnetron configurations²³² and various regions along the acceleration channel of HETs.^{74,107,227,230,236} A theoretical model of such waves has been developed⁷⁴ that accounts for plasma density, temperature, and magnetic field gradients as well as electron-neutral collisions and sheath boundary effects.^122,225 This model also predicts another type of unstable modes: resistive axial instabilities that propagate in the direction of the applied electric field.^115–117 These modes saturate due to ion nonlinearities (similar to the ion sound waves) and therefore can grow to large amplitudes,¹¹⁷ though they have generally smaller growth rates compared to the azimuthal modes.¹²¹ It has been suggested that the resistive axial modes play an important role as a trigger for the low-frequency ionization modes involving plasma and neutral density dynamics also known as the breathing modes.¹¹⁶ Nonlinear 2D z–θ simulations of gradient drift modes demonstrate complex self-organization of turbulence coexisting at small and large scales with a large level of anomalous transport with the effective Hall parameter, $H = ω c e / ν eff ≃ 10 − 30$⁠,⁷⁴ thus providing a self-consistent first principles anomalous transport without any additional closures. Fluid turbulence of partially magnetized plasmas (at the scales larger than ρ_e) demonstrates another important property: the inverse energy cascades toward longer wavelengths when the fastest instabilities occur on small scales which subsequently merge into large scale nonlinear structures such as shear flows, vortices, and streamers.^74,107 Figure 17 shows one typical example of plasma turbulence at an intermediate state in which the most unstable modes at the length scale of the order of λ ≤ 10ρ_e gradually develop into the large-scale structures (eventually saturating at the box length scales). The large-scale structures typically provide dominant contributions to the anomalous transport and may lead to the intermittent (bursting) avalanche-like transport events.²³⁵ These simulations suggest that large anomalous transport is possible, yielding an effective (time-averaged) Ω_e ≤ 15. However, because they have been performed in simplified 2D geometry and have neglected many important factors such as 3D effects, sheath boundaries, temperature, and neutral component evolution, they cannot be used yet to model real devices at the full scale.

In some E × B discharges, the disparate length and time scales between the plasma species allow for modeling methodologies that circumvent some of the aforementioned limitations of fully kinetic and fully fluid approaches. In HETs for example, the Knudsen number for the heavy species can, in many cases, exceed unity and the ions are not magnetized. Thus, the ion and neutral dynamics may be tracked easily and relatively inexpensively with particle-based methods, within a computational framework that still treats electrons as a distinct fluid. These so-called “hybrid” (fluid-particle) methods gained substantial popularity, especially during the early attempts to model the HET discharge over three decades ago.^207–210 Typically, the electron fluid equations employ formulations that include an anomalous contribution to account for the effects of turbulence on the transport of mass, momentum, and heat. One of the first hybrid codes was developed by Fife in the late 1990s. The code Dubbed HPHall,²⁰⁷ assumed Bohm's 1/B scaling for the cross field mobility.²¹¹ As more plasma measurements became available, however, it became clear that this scaling could not explain the behavior of the plasma in most regions' interior to the thruster channel. The argument against Bohm-driven transport was later strengthened by Mikellides²¹² and Lopez Ortega.²¹³ They employed a multivariable spatial model of the anomalous collision frequency in a 2D axisymmetric (r–z) code called Hall2De²¹⁴ and combined it with detailed plasma measurements in a HET to obtain a near-continuous, empirically derived, piecewise spatial variation of the frequency everywhere in that thruster. The latest mathematical formulation of the model is described in detail in Ref. 220. Typical Hall2De solutions for the electron collision frequency along the channel centerline (CL) of a magnetically shielded Hall thruster and comparisons with ion velocity measurements for two different magnetic field strengths are illustrated in Fig. 18.²³⁸ The comparisons underscore the dominance of the anomalous contribution to the total collision frequency over that from (classical) electron–ion (e–i) and electron–neutral (e–n) collisions. They also demonstrate the importance of the transport model in capturing spatial shifts of the discharge that have been well observed in the laboratory as operating conditions and/or the background in ground test facilities are varied (e.g., see Ref. 239 and references therein). Similar attempts with empirically derived anomalous resistivity models were made in other codes in the past as reported for example by Hagelaar²¹⁵ and by Scharfe.²¹⁶ In the Hall2De work, however, the combination of simulations and extensive measurements yielded some of the closest agreements we have achieved with time-averaged plasma measurements so far, at a relatively low computational cost and without statistical noise for neither ions nor neutrals. That is because, in this code, ions are treated using multifluid conservation laws and neutrals are tracked using line-of-sight formulations, not particle methods. The concern with statistical noise has been that it may interfere with the true plasma dynamics in these devices, which, as we alluded to earlier, span a wide range of spatiotemporal scale lengths. On the other hand, hydrodynamic approaches are unable to capture the details of the ion velocity distribution function (IVDF). The multifluid approach in Hall2De was an attempt to provide a marginally better approximation of the effects of the IVDF. More recently, Jorns¹⁶⁷ increased the efficiency and speed of this approach by employing machine-learning methods to obtain a solution for the anomalous resistivity. The approach takes advantage of the many advances made recently in machine-learning algorithms to replace with computer iteration the human iteration required in the determination of the anomalous collision frequency from plasma measurements. At this time, such 2D r–z (axisymmetric) solvers with empirical or machine-learning models for the electron transport offer the only feasible approach in providing practical support to the design and development of HETs.

Limited description of the kinetic phenomena is still possible within some fluid models by adding linear kinetic closures for the higher-order terms in moment equations. For example, it has been argued that the effects of linear Landau damping can be included in the fluid equations via the Hammett–Perking type closures for the viscosity terms,²²⁶ and effective schemes have been developed to solve them numerically.²²³ Such closures were shown to be effective in fluid modeling of turbulence and anomalous transport due to temperature-gradient modes in magnetized plasmas as well as in the treatment of ion–acoustic turbulence (IAT) and transport of heat in laser plasmas.²²¹ It is worth mentioning that in these closure schemes the turbulence and anomalous transport are modeled self-consistently with modified fluid equations. Other kinetic processes however have proven to be more challenging. In recent years, for example, the electron cyclotron drift instability (ECDI)^{142,177–179,182} has been proposed by Adam¹⁴² and later Coche¹⁷³ as a potential source of fluctuations and anomalous transport in HETs due to the strong E × B drift in these devices; see also Sec. IV. This is a microinstability which, in principle, requires a kinetic description though, to the best of our knowledge, no attempts have been made to apply Hammett–Perkins type closures within a fully fluid framework. An alternative closure approach was attempted recently in 2D (r–z) simulations by Lopez Ortega²⁰⁵ who modeled the ECDI in Hall2De²¹⁴ with an additional fluid-like equation for the evolution of the wave action associated with the propagation of wave packets. The equation allowed for the convection of the wave energy by the mean ion drift and a wave production term that was proportional to the growth of ion–acoustic waves and accounted for Landau damping. Closure to the hydrodynamics equations in Hall2De was achieved by relating the wave action to an anomalous collision frequency. The simulations predicted well the location of the acceleration region in unshielded and magnetically shielded versions of a HET. However, finer details such as changes in the plasma potential gradient within the acceleration region were not captured.²¹⁷ The work was an extension of previous attempts by Mikellides²¹⁸ to incorporate the ECDI in Hall2De simulations based on the hypothesis that the instability excites IAT which, in turn, enhances the effective collision frequency. Other larger scale hybrid approaches are being pursued, such as that by Joncquieres et al.²¹⁹ who are developing a 3D unstructured massively parallel particle-in-cell/fluid solver. In this work, the fluid part for electrons and ions is based on a 10-moment model, while the kinetic simulations are used as a reference solver for the sheath boundary conditions.

Fully kinetic simulations in three dimensions probably hold the greatest promise of resolving explicitly the wide-ranging spatial and temporal scales that persist in these devices. The major challenge here is computational resources which continue to limit our ability to perform global multiscale simulations. Fully fluid models as well as hybrid models require the information and input (possibly from supplementary kinetic models in reduced dimensions) on sheath boundary conditions and related particle and heat fluxes to the walls, as well as possible kinetic distortions of the electron distribution function at high energies. However, the challenge here is how to properly incorporate such kinetic effects within the fluid formulations.

A physical basis for the underlying mechanisms of the electron transport in E × B devices such as the HET has not been established yet. Hybrid simulations that have had the greatest impact on the design of these devices are performed largely using empirically derived models of the anomalous resistivity. In recent attempts to employ first principles closure models for the ECDI modes, the influence of the magnetic field on the transport was not fully captured and the effects of linear and nonlinear ion Landau damping remained unclear. Moreover, the entire concept of anomalous transport as a diffusive process that can be characterized by an effective collision frequency may be challenged by some experimental and computational data that indicate anomalous transport is strongly intermittent and avalanche-like. The characterization of nondiffusive but “blobby” transport events would then require different approaches such as self-organized criticality.²³⁷ It is also worth noting that the electron transport in regions of low electric field has been found to remain highly anomalous. This is evident, for example, as shown in Fig. 18. In such regions, instabilities that are excited by a strong E × B drift are not expected and, hence, it has been argued that mechanisms other than the ECDI must be acting.^236,240

Finally, as is the case with any transport theory, a major challenge is experimental validation. In the HET specifically, it has been shown in 2D (r–z) simulations that large changes in the cross field anomalous transport produce only small changes in the plasma in some regions of the acceleration channel and near-plume (Fig. 19). Because such changes are typically too small or impossible to detect by the current state of the art in plasma diagnostics, the validation of a transport model by experiment then becomes, very challenging, at best, ambiguous at worst. This underscores the need for more advanced plasma diagnostics in laboratory investigations of these devices.

Francesco Taccogna¹, Johan A. Carlsson², Andrew Tasman Powis², Sedina Tsikata³, and Konstantin Matyash⁴

¹CNR-Institute for Plasma Science and Technology, via Amendola 122/D 70126-Bari, Italy

²Princeton Plasma Physics Laboratory, Princeton NJ 08543, USA

³ICARE, CNRS, 1C avenue de la Recherche Scientifique, 45071 Orléans, France

⁴University of Greifswald, Greifswald, D-17487, Germany

A predictive model of E × B discharges requires a kinetic three-dimensional representation; see Sec. IV and Ref. 241. The electron transport across the magnetic field lines involves all three coordinates in a complex way. In particular, Hall thrusters (HTs),^4,44,141,242 characterized by a large electron current flowing along the azimuthal direction, are subject to the electron drift instability (EDI); see Sec. IV and Refs. 175 and 201. Recent work has shown how this azimuthal instability is influenced by the non-Maxwellian character of the electrons¹⁴⁷ and connected to the behavior of the plasma along the radial (parallel to the magnetic field line) and axial (accelerating field) directions. A collective Thomson scattering experiment¹⁴⁷ has demonstrated that the azimuthal mode has a non-negligible wave vector component along the magnetic field, while theory predicts that the presence of nonzero components along the magnetic field direction will modify the resonant comb-like nature of the dispersion relationship toward the ion–acoustic type; see Sec. IV and Refs. 175 and 201. In addition to the broadening of cyclotron resonances in the mode dispersion relationship, the inherently 3D nature of the instability is associated with lower growth rates in comparison to the 1D and 2D cases in linear kinetic theory. This has profound consequences: in order to correctly determine the contribution of this particular instability to electron transport (particularly in comparison to other instabilities), a 3D treatment of the instability, both in theory and in simulations, appears to be critical. Although it would seem that the Near-Wall Conductivity (NWC) results in a minor contribution to the anomalous electron transport across the magnetic field lines, secondary electron emission, when modeled in Particle-in-Cell (PIC) simulations,^{24,146,169,243} may represent a saturation mechanism for EDI. Resolution of the axial direction represents a more natural way for the instability to be convected into the near plume region.^142,176,187 Finally, a fully kinetic representation is required due to the fact that ion kinetics contribute to the saturation of the wave by ion-wave trapping and ion axial acceleration. The situation is further complicated by the emergence of large scale low frequency azimuthal spoke modes^{73,75,76,78,93,108,168} (see Fig. 20). Due to their highly nonlinear, turbulent, and global nature, these structures have continued to evade a clear theoretical understanding. However, their observed correlation with enhanced electron transport makes them an important feature that must be captured by any numerical model.

Along with the need to capture the physics in three dimensions, there is also evidence that low-dimensional (1D and 2D) models often show their limitations by leading to inconsistent and controversial results. 1D azimuthal^106,163,191 and 2D azimuthal–radial^{24,146,190,243} models must mimic the particle transport along the axial coordinate: to avoid unphysical energy growth, particles are reinjected as new colder particles when they have traveled a prescribed axial distance. This numerical artifact strongly distorts the final results and is a possible candidate to induce discordant results about the transition from electron cyclotron resonances toward the ion–acoustic type. Moreover, the accelerating axial electric field E_z is externally imposed and fixed, while one expects that an increase in the local conductivity due to the anomalous transport would lead to a decrease in the local axial electric field, which in turn would reduce the local electron azimuthal drift velocity and thus significantly change the instability behavior. 2D (axial–azimuthal) models^142,176,187 are also subject to important shortcomings, such as neglecting components of the wave parallel to the magnetic field and effects of the nonlinear coupling of the EDI with secondary electron emission from the walls. Generally, reduced dimensional models always show stronger instability characterized by large amplitude oscillations, as shown in Fig. 21 comparing the electron and ion densities from 3D and 2D PIC models. This finally leads to an artificially larger cross field mobility (more than five times) with respect to 3D model estimations and experimental measurements. Another observation of note is that the single dominant length scale for the instability, shown in Fig. 21 for the 2D case, is no longer sharply defined in the 3D case. This observation may be compatible with both experimental and 3D linear kinetic theory analyses showing the instability to be excited nondiscretely, across a range of wavenumbers.

The quantitative description of the intrinsically nonequilibrium, nonlocal nature of anomalous transport in HTs requires the development of an efficient kinetic numerical tool capable of capturing both electron and ion time scales under the electrostatic (ES) approximation in a three-dimensional spatial domain which includes both the discharge channel and the near-field plume regions.

In order to resolve the Debye length under typical HT physical parameters (λ_D ≈ 50 μm), a three-dimensional PIC model requires ∼10³ mesh nodes per coordinate and time steps on the order of picoseconds to satisfy the particle-cell-transit Courant–Friedrichs–Lewy condition. The total simulation time necessary to capture the nonlinear and saturation phases of the E × B EDI is about 50 μs; an even longer simulation may be required to resolve the temporal dynamics of large scale azimuthally propagating structures^{73,75,76,78,93,108,168} (up to hundreds of μs: 10⁸ time steps). In addition, an acceptable statistical noise level corresponds to a number of particles per cell N_ppc >10³;^163,187,191 numerical heating^244–246 is often responsible for the destruction of electron cyclotron resonances and facilitates an artificial transition to the ion acoustic-type instability. A 3D ES-PIC algorithm, using linear interpolation, requires ∼200 floating-point operations (FLOPs) per particle (Lagrangian phase) and ∼40 FLOPs per grid cell (Eulerian phase), which leads to a total number of required calculations on the order of ∼10²¹ FLOPs. Finally, simulation of a time-dependent six-dimensional problem produces a large and complex set of output data whose analysis (spectral), visualization and animation require highly efficient and dedicated tools.

The natural progression of increasing computational power (exascale, i.e., 10¹⁸ FLOPs per second) and improved algorithm efficiency will enable the development of full-resolution 3D ES-PIC models of HTs within the coming decade, as is routinely performed within the laser-plasma accelerator PIC modeling community with OSIRIS,²⁴⁷ PICADOR,²⁴⁸ SMILEI,²⁴⁹ and PICSAR.²⁵⁰ 

A major hurdle that has been already overcome in the past decade is the vastly improved scalability of Poisson equation solvers and their efficient implementation.^251–255 However, in PIC simulation, the Lagrangian phase (loops over particles, i.e., projector, pushing, interpolator, and Monte Carlo collisions) is usually much more expensive than the Eulerian phase (loops over grid nodes). In particular, the projection operator (deposition of the charge density from the particle position to the mesh nodes) can take up to 60% (depending on the order of the particle interpolation used) of the whole particle advancing time. In order to speed up the Lagrangian phase, two different approaches can be adopted.

The first approach can already be applied to serial versions of PIC codes by implementing an efficient data structure and vectorization ideal for single instruction multiple data (SIMD) registers, in order to optimize memory access. This goal can be achieved by particle tiling and sorting methods.^256,257 The first consists of the fact that particles are placed in tiles that fit in cache, while the second allows particles belonging to the same cell to be contiguous in memory in order to maximize cache efficiency (avoiding random memory access and facilitating data reuse) and to be easily vectorized.

The second approach is an efficient use of hybrid parallelism: the particle decomposition with Open Multiprocessing (OpenMP) within a node and the domain decomposition with a Message Passing Interface (MPI) between the nodes, in order to improve the scalability up to 10⁶ processors. This will be crucial on future architectures that will have less available memory per core. Strictly related to efficient parallelism is the dynamic load balancing technique in order to balance the computational load among the different compute units. As seen above, since most of the computational load is proportional to the number of particles, the effort mainly consists of balancing the number of particles per compute unit.

Generally, most ES-PIC algorithms need to be completely rewritten and data structures need to be completely revolutionized in order to use the advantages of new architectures.

Immediately, two simpler approaches are readily available. One way to reduce the highly demanding computational requirements of a 3D PIC model of HTs is to use ion-subcycling routines (pushing ions with a larger time step) and to decrease the azimuthal length by imposing a smaller periodicity (Fig. 22 shows results corresponding to an azimuthal sector reduced by 40). However, this can go against the evidence that E × B EDI evolves toward low wavenumbers (e.g., wavelength comparable to the entire azimuthal size) through an inverse cascade process.^142,190 Another option is to perform simulations on a miniaturized HT²⁵⁸ that also meets the need of the community to scale down HTs for CubeSat micropropulsion.²⁵⁹ 

The physics characterizing E × B plasmas, and in particular, Hall thrusters has a nonequilibrium, nonlocal character and involves self-organized three-dimensional structures. A predictive numerical model has to describe electron and ion time scales over a region that extends from the discharge channel up to the near field plume region.

One- and two-dimensional models are affected by their limitations and often show contradictory results. Progress in high-performance computing (HPC) and the availability of more powerful supercomputers and algorithms will enable the development of 3D fully kinetic PIC models, representing real numerical experiments devoid of any artificial hypothesis. The complete understanding of electron transport will lead to a new era in the technological development of E × B plasma devices and designs based on an empirical approach will give way to code-based refined optimization. As has been done in many other engineering disciplines, predictive design and optimization via computer-based techniques will assist and eventually replace empirical methods.

Renaud Gueroult¹, and Nathaniel J. Fisch²

¹LAPLACE, Université de Toulouse, CNRS, INPT, UPS, Toulouse, France

²Princeton Plasma Physics Laboratory, Princeton University, Princeton, NJ 08543, USA

There is a long history of crossed-field (or E × B) configurations to separate charged particles based on mass. Although crossed-field mass spectrometers²⁶⁰ first demonstrated these capabilities, it was quickly realized that the throughput achievable with non-neutralized ion beams (such as those employed in the calutron²⁶¹) was severely limited by space charge effects and instabilities. The smaller throughput necessary for isotope separation, however, led to the development of plasma isotope separators,²⁶² where crossed-field configurations again played a role. Specifically, crossed fields and their associated drift were used to produce plasma rotation²⁶³ in plasma centrifuges,^264,265 as illustrated in Fig. 23(a).

More recently, there has been a growing recognition that elemental separation based on mass offers highly promising conceptual solutions to several societal high-impact applications, including nuclear waste cleanup,²⁶⁶ spent nuclear fuel (SNF) reprocessing,^267–271 and rare earth element (REE) recycling.²⁷² For these applications, conventional chemical separation techniques are inefficient since one needs to separate feeds of varying chemical composition (REE recycling and nuclear waste clean-up) or elements with similar chemical properties (SNF reprocessing). Further to these inefficiencies, chemical separation often leads to the production of secondary, and possibly hazardous, waste. On the other hand, as illustrated in Fig. 24, it happens that the feed conveniently splits into two groups based on atomic mass, supporting the interest of elemental mass separation. Compared to chemical techniques, another anticipated advantage of plasma separation lies in its much smaller environmental footprint. Thanks to these characteristics, plasma separation may be economically competitive for certain applications,^266,272 particularly considering that energy costs from solar may become significantly cheaper. However, the separation needs of these new applications differ significantly from those of isotope separation, thereby offering new challenges.²⁷³ 

To be useful for these new applications, plasma separators not only have to separate elements with large mass differences (≥tens of amu), but also, and importantly, to operate at high throughput (≥10⁴ kg/yr). It is also advantageous if the separation methods do not perform well on isotopes to lessen proliferation risks. Since these capabilities are well beyond those accessible to isotope separators, new plasma mass filter concepts are called for.

Here, also, the rich physics of crossed-field configurations offers many opportunities and various conceptual solutions have been suggested in the last decade.^267,275 In regimes where both ions and electrons are magnetized, crossed-field configurations can be designed to drive rotation. By exploiting rotation in novel ways, conceptual solutions for mass separation can then be found beyond plasma centrifuges^275–277 as shown in Figs. 23(b) and 23(c). Alternatively, in regimes where ions are unmagnetized, but electrons are magnetized, crossed-field configurations can in principle be designed to separate the ions of a rotating annular ion beam,^278,279 or to further exploit differences in gyroradius.²⁸⁰ 

In all of these concepts, large mass differences are advantageous, and, in view of proliferation risks, usefully necessary. On the other hand, high-throughput operation has been shown to bring a new set of physics and engineering challenges,²⁸¹ which are yet to be addressed.

Upstream of the elemental separation stage discussed above, one needs to consider how to produce the high-density plasma required for separation from solids or liquid input streams. Many challenges are found at this stage. First, how to feed into a vacuum chamber and then ionize, grams of material per second while minimizing energy cost? Laser evaporation of solid targets and injection of micrometer size powder have been suggested so far, but this question ought to be addressed in detail. Another constraint on the plasma source is the need to maximize charge uniformity. Indeed, since separation often results from differences in the ion charge to mass ratio Ze/m_i and not in the ion mass alone, effective separation requires producing a plasma with Z as close as possible to homogeneous and uniform across ion species. Furthermore, low Z materials are to be preferred to minimize radiation losses, which translates into electron temperature T_e of at most a few eV for most atoms. Finally, related to this challenge is the open question of how molecules will impact separation. While significant molecular ion components are expected due to the low T_e, the detailed composition of plasmas formed from complex mixtures is far from understood and so is the dynamics of these compounds. In addition, cross section data for molecular processes is likely not readily available for some of the uncommon elements (e.g., the sesquioxides) found in these applications.

The separation stage also offers challenges. Assuming a singly ionized atomic plasma could be produced, separation in crossed-field concepts is conditioned upon externally applying a specific potential profile in the direction perpendicular to the magnetic field. Fundamentally, the ability to support this perpendicular potential profile depends on the ratio of perpendicular to parallel conductivity σ_⊥/σ_∥.²⁸² However, many different driving mechanisms (e.g., collisions with neutral, instabilities and turbulence,²⁸³ magnetic fluctuations,²⁸⁴ ion viscosity²⁸⁵) are known to contribute to σ_⊥ depending on the operating plasma conditions, and new driving mechanisms continue to be uncovered.^286,287 Demonstrating the practicality of crossed-field mass filter concepts hence hinges on a comprehensive understanding of perpendicular conductivity, which we believe will be best obtained by a combination of modeling and experiments.^274,288 Another outstanding issue in the presence of neutrals is the possible upper limit set on the rotation speed by the critical ionization velocity phenomenon.²⁸⁹ 

The challenges briefly discussed above are only a small subset of the many open questions that remain to be addressed to demonstrate high-throughput plasma separation as a practical process (see Refs. 267 and 281 for a detailed discussion). Yet, we believe that the promise plasma separation holds for many outstanding societal challenges is a compelling motivation to tackle these questions. In addition, progress toward many of these scientific and engineering goals will benefit applications beyond plasma separation. For instance, the basic question of perpendicular conductivity in a magnetized plasma is also central to recently proposed promising schemes for magnetic confinement fusion.^290,291 Understanding how perpendicular electric fields can be externally applied will more generally benefit the large number of applications making use of cross field configurations.

Anne Bourdon¹, Pascal Chabert¹, Alexander V. Khrabrov², Andrew Tasman Powis², Ioannis G. Mikellides³, and Igor D. Kaganovich²

¹Laboratoire de Physique des Plasmas, CNRS, Ecole Polytechnique, Sorbonne Université, 91120 Palaiseau, France

²Princeton Plasma Physics Laboratory, Princeton NJ 08543, USA

³Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, 91109, USA

The ultimate objective for developing computer simulations of complex physical systems is to use these simulations as a predictive tool for science and engineering design. This would reduce the number of costly experiments required for the development of applications involving low-temperature plasmas. For example, currently, the design and development of electric thrusters are semi-empirical with long (on the order of 10 000 h) and expensive lifetime tests. Applications that employ low-temperature plasmas generally involve complex multiphysics and multiscale processes. Developing accurate and predictive simulation tools remains an active area of research. Critical to these goals is the need to verify and validate codes used for the predictive modeling of low-temperature plasmas.

In other fields, rigorous verification and validation (V&V) procedures have been implemented for decades. Computational fluid dynamics (CFD) codes, which are used for the design of commercial airplanes, cooling systems, and to simulate wind loads around buildings and bridges apply tests of code reliability by regulation to adhere to professional engineering standards.^292,293 Prominent engineering journals require V&V testing as part of their editorial policy before acceptance of papers.²⁹⁴ 

In order to develop predictive capabilities for low-temperature plasmas, several necessary steps have to be performed.

• Develop a minimum complexity self-consistent mathematical model that describes the effects to be predicted.

• Develop numerical codes and check these codes for bugs by benchmarking with other codes, and with analytical or manufacturing solutions, and verify that convergence is achieved.

• Assemble a reliable set of atomistic physics and plasma–surface interaction data required for modeling.

• Perform validation tests of the assumptions and simplifications of the model by comparing simulation results with available or dedicated experiments capable of exploring a wide range of parameters. The comparisons must be sufficiently comprehensive^296,297 in order to avoid the bias which may occur when “measurements data and modeling results agree if they are performed in the same Lab” and because a narrow set of measurements data can be always matched by variation of adjustable parameters available in the models.

Code verification²⁹⁸ is used to assess the numerical accuracy of a mathematical model. The first test should include convergence tests with number of grid points, values of time steps for fluid, hybrid and particle-in-cell codes, and number of particles for particle-in-cell and hybrid codes. Note that some of the studies, for example the number of particles in particle-in-cell (PIC) codes, are challenging due to slow convergence, see, e.g., Ref. 295. Proving the convergence of a solver can be accomplished through comparison with simple limiting test cases or analytical solutions. Another way to verify the code is to use manufacturing solutions, see, e.g., Ref. 295.

Code validation proves the predictive capabilities of the developed model by comparing simulation results with experimental data.²⁹⁸ This is usually more challenging and demanding because it requires a well-organized systematic campaign of combined experimental and numerical studies in order to convincingly demonstrate validation.^296,297

An early example of such an effort for low-temperature plasmas is the work of Surendra²⁹⁹ for a capacitively coupled radio frequency discharge, in which results from twelve different codes (including four particle-in-cell codes) were compared with each other (benchmarking) and with experimental data (validation). At the end of the 1990s, a standardized RF discharge experiment—a so-called GEC reference cell³⁰¹—was also specifically designed for validation and is still used as shown in the recent study performed in Ref. 302. Recently, there have been different works on the benchmarking of codes. In particular, Turner et al.³⁰⁰ performed a parametric study of several 1D particle-in-cell codes for verification of the Poisson solver, the particle pusher, the wall boundary conditions and the ionization, and Monte Carlo collision modules. The comparison of the codes was performed until an error of less than 1% was obtained. Therefore, the results of this benchmark can be considered as a reference solution for other 1D PIC codes. An international effort is under way to preform 2D benchmarks for (nonmagnetized) low-pressure low-temperature plasmas.³⁰³ Recently, there has also been a renewed interest to carry out studies on plasma sheaths in low-pressure low-temperature plasmas and to use PIC simulations to guide the development of theoretical models.³⁰⁴ 

The benchmarking of fluid codes could be more involved than that of PIC codes. Indeed, a large variety of fluid models are used, ranging from drift-diffusion models to those which maintain a larger number of moments for fluid closures (13 and above). Furthermore, there is a large variety of numerical schemes used to solve the set of fluid equations. For example, six fluid simulation codes were recently benchmarked for simulation of positive streamers at atmospheric pressure.³⁰⁵ Three test cases of increasing complexity have been studied. A reasonable agreement between the results of the different codes has been obtained, also showing the difficulty to of obtaining a reference solution under given conditions.

For low-pressure low-temperature magnetized plasmas, of particular interest for electric propulsion or magnetrons, a significant international effort is currently underway to benchmark PIC and fluid codes on test cases³⁰⁶ corresponding to typical conditions of Hall thrusters. The problem has, however, been sufficiently simplified for ease of benchmarking. The first results obtained by different international groups on the three test cases have been presented at the E × B workshop organized at Princeton in November 2018.³⁰⁷ In particular, the most significant effort has been placed on the 2D test-case dedicated to the simulation of microinstabilities induced by the large E × B electron drift in the acceleration zone of a Hall thruster. The test-case is two-dimensional in the axial–azimuthal direction, close to that proposed in Ref. 187. Seven international groups have compared their independently developed PIC code results for this 2D test-case and obtained good agreement within 5% difference between all the codes; the results of this study are published recently in Ref. 204. Therefore, these results can be considered as a reference solution for other PIC codes.