1. Introduction
Plasmas subject to mutually perpendicular electric and magnetic fields – commonly known as E × B or cross-field plasmas – constitute a class of partially magnetised discharges that underpin a wide variety of applications and technologies. They are encountered in plasma-assisted material processing systems such as magnetrons, in spacecraft propulsion devices including Hall thrusters and in laboratory plasma sources such as Penning discharges. From these low-temperature discharges to magnetic confinement fusion devices, a common challenge shared across magnetised plasma systems – despite the wide disparity in operating conditions – is cross-field transport of charged species, particularly electrons, driven by microinstabilities. In these plasma systems, coupled phenomena exist that span multiple spatial and temporal scales (Kaganovich et al. Reference Kaganovich2020; Boeuf & Smolyakov Reference Boeuf and Smolyakov2023). Their global behaviour – and, consequently, the confinement and transport dynamics in general – arise from the intricate interplay between instabilities, nonlinear wave interactions and turbulence. The coupling among these processes governs how plasma species are transported, often in ways that are difficult to predict or control. Advancing the understanding of these coupled dynamics and developing means to regulate them remain central goals across both fundamental and applied plasma physics, with direct implications for the efficiency, stability and longevity of magnetised plasma systems.
In partially magnetised E × B plasma devices, a key factor determining overall operational performance is the confinement of electrons and the extent to which their motion across magnetic field lines can be limited. Effective confinement is essential for maintaining efficiency in systems that rely on electron-impact ionisation, such as Hall thrusters, where excessive cross-field transport leads to power losses and reduced thrust efficiency. It is now well established that this transport is not purely collisional, but is largely instability-driven, arising from fluctuations that develop along the azimuthal direction coinciding with the electron E × B drift (Adam, Héron & Laval Reference Adam, Héron and Laval2004; Lafleur, Baalrud & Chabert Reference Lafleur, Baalrud and Chabert2017; Boeuf & Garrigues Reference Boeuf and Garrigues2018; Brown & Jorns Reference Brown and Jorns2023). Correlations between oscillations of electron density and azimuthal electric field give rise to an anomalous momentum exchange between electrons and ions (Cavalier et al. Reference Cavalier, Lemoine, Bonhomme, Tsikata, Honore and Gresillon2013; Lafleur et al., Reference Lafleur, Baalrud and Chabert2016a,b ), often interpreted as an instability-enhanced friction force. These azimuthal instabilities, including the electron cyclotron drift instability (ECDI) and the modified two-stream instability (MTSI), therefore play a key role in governing cross-field electron mobility. Understanding how such modes evolve, interact and can potentially be mitigated remains central to advancing both the predictive modelling of plasma transport and the practical optimisation of the operation of E × B plasma devices.
Among various approaches proposed to influence plasma instabilities and associated transport, modulating the anode voltage has emerged as an effective means of controlling low-frequency macroscopic oscillations, as observed, for example, in Hall thrusters. Simmonds et al. (Reference Simmonds, Raitses, Smolyakov and Chapurin2021) investigated the use of a modulated anode voltage to manipulate the breathing mode of a Hall thruster, demonstrating modest performance improvements. Romadanov et al. (Reference Romadanov, Raitses, Diallo, Hara, Kaganovich and Smolyakov2018b ) achieved active control of the so-called spoke mode using a similar approach. In both cases, a modulated discharge voltage was applied to the plasma at relatively low driving frequencies – of the order of kilohertz – corresponding to the characteristic frequencies of the breathing and spoke modes, respectively.
Externally driven breathing-mode oscillations under discharge-voltage modulation have been characterised in more detail in both experiments and simulations. Romadanov et al. (Reference Romadanov, Raitses and Smolyakov2018a ) reported that voltage modulation can drive a highly coherent breathing-mode response, with distinct linear and nonlinear response regimes and a resonant dependence on the driving frequency. Perales-Díaz et al. (Reference Perales-Díaz, Domínguez-Vázquez, Fajardo and Ahedo2023) performed axisymmetric hybrid fluid/particle-in-cell simulations of a magnetically shielded Hall thruster under discharge-voltage modulation, showing that modulation can control the oscillation frequency when driven near the natural breathing mode and clarifying the role of the electron-temperature/ionisation-cycle response in determining the driven dynamics and net performance trends.
Tejeda et al. (Reference Tejeda, Reza, Faraji and Knoll2022) investigated high-frequency RF forcing in a Hall-thruster configuration by imposing a sinusoidal oscillation of the cathode reference potential in a pseudo-two-dimensional (pseudo-2-D) axial-azimuthal PIC model. Their analysis is primarily focused on electron-cyclotron resonance and RF power deposition (heating), including cyclotron-related spectral features and frequency-locking behaviour. While they reported improved global performance under near-resonant forcing, it is not possible to unambiguously separate whether these gains arise from reduced cross-field electron transport or from enhanced ionisation, and this attribution is further complicated by the non-physical ionisation treatment and simplified pseudo-2-D PIC formulation adopted in that work. A comparable strategy has also been adopted to suppress drift-wave instabilities in Q-machines, as demonstrated in the experimental and theoretical works of Nishida, Tanibayashi & Ishii (Reference Nishida, Tanibayashi and Ishii1970) and Weiss & Morrone (Reference Weiss and Morrone1976).
Complementary studies have examined the effects of azimuthal neutral-density variations and azimuthal magnetic-field strength on electron transport. Bak et al. (Reference Bak, Kawashima, Komurasaki and Koizumi2019) showed experimentally that azimuthal non-uniformities in the propellant density can induce axial-azimuthal electric fields that enhance electron transport. Bak et al. (Reference Bak, Kawashima, Romanelli and Komurasaki2022) further showed that imposing an azimuthally non-uniform
$B_{r}(\theta )$
generates a correlated
$\mathrm{z}$
–
$\theta$
structure and finite
$E_{\theta }$
, such that drift-driven transport (notably
$E_{\theta }\times B_{r}$
and diamagnetic contributions) becomes dominant, and the effective axial electron mobility increases relative to an azimuthally uniform magnetic field.
An alternative theoretical approach was proposed by Kapulkin & Prisnyakov (Reference Kapulkin and Prisnyakov1995), who investigated suppressing the electron-drift instability by machining grooves into the channel walls to modify the electron surface-collision frequency. Most recently, Rose & Knoll (Reference Rose and Knoll2023) demonstrated, in one-dimensional simulations, that an axial electric field modulation at the fundamental frequency of the ECDI can influence the resulting electron transport.
In this study, we investigate how externally applied electrostatic modulation can change plasma response in one-and two-dimensional E × B discharge configurations, representative of Hall thruster geometry. The objective is to assess how modulations applied at various frequencies and amplitudes affect the instability spectra and the resulting wave-induced electron transport, and to identify modulation conditions under which electron mobility is reduced or enhanced. This is supported by a bicoherence analysis to characterise nonlinear mode interactions and their role in shaping the observed spectra, and by extension, electron mobility variations.
The study hence provides: (i) a systematic mapping of plasma response to external modulation across a wide range of frequencies and amplitudes; (ii) quantitative evaluation of the associated changes in instability characteristics and electron mobility; and (iii) assessment of the underlying coupling processes through phase-correlation analysis of azimuthal electric field fluctuations. Collectively, these investigations provide new perspectives on the mechanisms of instability control and transport regulation in E × B plasma systems.
2. Theoretical background
2.1. Overview of the characteristics of dominant azimuthal instabilities and resulting transport in the studied set-ups
In E × B plasma configurations, a range of azimuthal electrostatic instabilities exists. A key azimuthal instability in these discharges like in a Hall thruster is the ECDI – a high-frequency (1–10 MHz) kinetic instability with wavelengths on the scale of the electron Larmor radius. This instability propagates mainly along the azimuthal (E × B) direction. In typical conditions of Hall-thruster-like E × B discharges, the phase velocity of the waves is lower than the azimuthal drift of the electrons, allowing the waves to grow by extracting energy from the drifting electrons. The ECDI excites as a result of coupling between Doppler-shifted electrostatic Bernstein modes and the ion acoustic waves (Lafleur et al. Reference Lafleur, Baalrud and Chabert2017). Extensive research has been conducted on ECDI over the years through theoretical models (Ducrocq et al. Reference Ducrocq, Adam, Héron and Laval2006; Cavalier et al. Reference Cavalier, Lemoine, Bonhomme, Tsikata, Honore and Gresillon2013), numerical simulations (Taccogna et al. Reference Taccogna, Minelli, Capitelli and Longo2012; Boeuf Reference Boeuf2014; Coche & Garrigues Reference Coche and Garrigues2014; Reza et al., Reference Reza, Faraji and Knoll2023b,c , Reference Reza, Faraji and Knoll2024a,b ; Villafana, Cuenot & Vermorel Reference Villafana, Cuenot and Vermorel2023) and experimental studies (Tsikata et al. Reference Tsikata, Honore, Lemoine and Gresillon2010; Tsikata, Honore & Gresillon Reference Tsikata, Honore and Gresillon2013). The significance of the ECDI roots in the multitude of effects it has on plasma behaviour, in particular, increased electron cross-field transport (Lafleur et al. Reference Lafleur, Baalrud and Chabert2016a ; Croes et al. Reference Croes, Lafleur, Bonaventura, Bourdon and Chabert2017; Brown & Jorns Reference Brown and Jorns2023; Reza et al. Reference Reza, Faraji and Knoll2023c ) and significant plasma heating (Lafleur et al. Reference Lafleur, Baalrud and Chabert2016b ; Janhunen et al. Reference Janhunen, Smolyakov, Chapurin, Sydorenko, Kaganovich and Raitses2018b ; Reza et al. Reference Reza, Faraji and Knoll2024a,b ).
The 2-D approximation of the ECDI’s dispersion relation, perpendicular to the applied magnetic field (in the axial-azimuthal plane), reveals that the ECDI unstable modes develop close to the resonance condition specified by (Ducrocq et al. Reference Ducrocq, Adam, Héron and Laval2006; Lafleur et al. Reference Lafleur, Baalrud and Chabert2016a )
where,
$k_{z}$
represents the wave’s azimuthal wavenumber (in
$\mathrm{m}^{-1}$
),
$\omega _{ce}$
is electron cyclotron frequency (in
$\mathrm{rad}\,\mathrm{s}^{-1}$
),
$V_{de}$
is electrons’ azimuthal drift velocity (in
$\rm m \, s^{-1}$
), and
$E$
and
$B$
are the magnitudes of the axial electric field (
$\mathrm{V}\,\mathrm{m}^{-1}$
) and radial magnetic field (in
${\rm T}$
), respectively.
Another influential instability in E × B plasmas is the MTSI. McBride (Reference McBride1972)’s theoretical analysis demonstrated that this mode behaves as a fluid-like instability, unaffected by the electron-to-ion temperature ratio (
$T_{e}/T_{i}$
). Numerical studies by Janhunen et al. (Reference Janhunen, Smolyakov, Sydorenko, Jimenez, Kaganovich and Raitses2018a
) and Taccogna et al. (Reference Taccogna, Minelli, Asadi and Bogopolsky2019) identified this instability as a relatively long-wavelength mode with a non-zero wavenumber in the radial direction. The presence of this mode in the systems studied has been shown to contribute to an enhanced inverse energy cascade phenomenon (Koshkarov et al. Reference Koshkarov, Smolyakov, Raitses and Kaganovich2019).
Boeuf & Smolyakov (Reference Boeuf and Smolyakov2023) presented the MTSI dispersion relation by adapting the kinetic ECDI dispersion relation, specifically for the case of cold electrons (
$T_{e}\rightarrow 0$
) and including a finite wavenumber component along the magnetic field direction. Petronio et al. (Reference Petronio, Tavant, Charoy, Alvarez-Laguna, Bourdon and Chabert2021, Reference Petronio2023) obtained a simplified 2-D fluid dispersion relation for the MTSI in the radial-azimuthal plane, assuming isothermal electrons with isotropic temperature and treating ions as cold and unmagnetised. They then derived a simple relationship between the radial and the azimuthal wavevector components of the MTSI’s fastest growing mode as (Petronio et al. Reference Petronio, Tavant, Charoy, Alvarez-Laguna, Bourdon and Chabert2021)
\begin{align}{k_z} = {\frac{1}{2\pi }}\sqrt {{\frac{e}{{m_e}}}{\frac{{B^2}}{E}}{k_x},} \end{align}
where
$k_{z}$
and
$k_{x}$
denote the wave’s azimuthal and radial wavenumber components, respectively.
The azimuthal instabilities generate fluctuating azimuthal electric fields
$\tilde{E}_{z}$
and density perturbations
$\tilde{n}_{e}$
, whose average correlation induces a non-zero time-averaged azimuthal force
$q \langle \tilde{n}_{e}\tilde{E}_{z}\rangle$
on the electrons. This fluctuation-induced term acts as an effective friction between electrons and ions, and enhances cross-field (axial) electron transport much beyond the classical (collisional) value. This instability-driven cross-field electron mobility can be described by (Lafleur et al. Reference Lafleur, Baalrud and Chabert2016b
, Reference Lafleur, Baalrud and Chabert2017)
where
$n_{e}$
is the average background electron number density.
2.2. Wave excitation in plasmas and energy transfer mechanisms
Let us consider an externally applied, temporally oscillatory electric field, often termed electrostatic excitation or modulation, superposed on the stationary axial electric field, as will be described in §§ 3.1 and 3.2. Such perturbations couple to the plasma linearly – by directly engaging natural plasma wave spectrum – and nonlinearly – by mediating multi-wave interactions and wave–particle resonances that redistribute energy across scales. These couplings determine how injected energy is distributed among the natural eigenmodes of the system and, consequently, impact the instabilities spectrum and transport pathways.
Electrostatic modulation can also modify the phase relation between the number density and the azimuthal electric field fluctuations, changing transport behaviour. Modulation-induced phase detuning can weaken the cross-correlation
$\langle \tilde{n}_{e}\tilde{E}_{z}\rangle$
that underpins instability-enhanced electron transport as per (2.3). We will further examine this point in § 5.2.
2.2.1. Linear response and near-resonant wave coupling
In the linear limit, the plasma response to an external oscillatory electric field can be described by small perturbations superposed on the steady state. Considering a harmonic excitation with angular frequency
$\omega _{M}$
and wavenumber vector
$\boldsymbol{k}_{M}$
, the first-order perturbation of the electrostatic potential satisfies the linearised Vlasov–Poisson system, yielding the response function
where
$D(\omega , \boldsymbol{k})=0$
represents the plasma’s natural dispersion relation and
$\tilde{S}$
denotes the source associated with the imposed modulation. The transfer function
$T(\omega , \boldsymbol{k})=D^{-1}(\omega , \boldsymbol{k})$
peaks near the eigenfrequencies of the system; as a result, near-resonant modulation at
$(\omega _{M}, \boldsymbol{k}_{M})$
selectively amplifies natural modes with
$(\omega , \boldsymbol{k})\approx (\omega _{M}, \boldsymbol{k}_{M})$
. When
$\omega _{M}$
lies close to an instability’s peak (e.g. the ECDI band), the modulation can add coherent energy to that band.
2.2.2. Nonlinear wave–plasma coupling and spectral energy pathways
At finite excitation amplitudes, the plasma response becomes nonlinear and the imposed oscillation can no longer be treated as a small perturbation described by the linear model.
The nonlinear response of the system can lead to ‘frequency locking’. In this state, the oscillation frequency of a natural mode becomes synchronised to that of the applied excitation, even when the external modulation frequency
$\omega _{M}$
differs slightly from the natural eigenfrequency
$\Omega_{j}$
.
From a general dynamical-systems perspective, frequency locking between two coupled oscillators can be described by the general Adler equation (Adler Reference Adler1946)
where
$\Delta \theta$
is the instantaneous phase difference,
$\Delta \omega =\Omega_{j}-\omega _{M}$
is the detuning between the natural and modulation frequencies, and
$K$
is a nonlinear coupling strength.
Frequency locking occurs when the rate of change of the phase difference becomes zero, i.e.
${\text{d}(\Delta \theta )}/{\text{d}t}=0$
. In this condition, the phase difference settles to a constant value, which means the oscillator and the external modulation signal are oscillating at the same frequency. From the Adler equation (2.5), this stable condition is satisfied when
Because the sine function’s value is constrained between −1 and 1, a steady-state solution exists only if
$|\Delta \omega |\leq K$
, corresponding to a so-called ‘locked regime’. In frequency-locked states (Adler Reference Adler1946), the spectral response exhibits a narrowband structure centred at the modulation frequency (Kurokawa Reference Kurokawa1973), often accompanied by a reduction in broadband turbulence within the locked frequency band. This behaviour reflects the conversion of stochastic phase fluctuations into coherent oscillations synchronised with the external modulation (Razavi Reference Razavi2004).
In an E × B discharge,
$K$
quantifies the strength of phase pulling exerted by the imposed modulation on the dominant azimuthal mode. Since the modulations act primarily through perturbing the electron drift
$V_{E}=E_{y}/B_{x}$
, the coupling is expected to increase with the drift perturbation amplitude,
$\delta V_{E}\sim \delta E_{y}/B_{x}$
(i.e. stronger modulation and/or smaller
$B$
). For a given
$\delta E_{y}/B_{x}$
, the effective
$K$
is reduced when broadband turbulence produces strong phase diffusion and is enhanced when the response remains relatively coherent and single-mode-dominated. Strong phase coherence, hence ‘locking’ behaviour, is expected when the effective coupling induced by modulation is sufficient (meaning sufficiently large
$K$
) to overcome the frequency mismatch between the dominant instability and the modulation.
The modulating wave can interact nonlinearly with both waves and particles, opening additional energy-transfer channels:
-
(i) Three-wave (and higher-order) interactions: quadratic nonlinearities in the fluid and kinetic plasma descriptions permit three-wave (triad) couplings that satisfy matching conditions
(2.7)
\begin{align}\omega_M\approx\omega_1\pm \omega_2,\quad \boldsymbol{k}_M\approx \boldsymbol{k}_1\pm \boldsymbol{k}_2.\end{align}
Hence, quadratic nonlinearities multiply existing waves and generate sidebands at sum or difference frequency and wavenumber, populating harmonics and sub-harmonics, enabling spectral broadening.
-
(ii) Nonlinear wave–particle interaction: when
$\omega _{M}$
approaches a particle resonance (e.g. electron cyclotron,
$\omega _{ce}$
), the modulation can alter single-particle dynamics: cyclotron heating increases electron Larmor radius, weakens magnetisation and changes collisionality-free cross-field mobility. In the E × B discharge turbulence context, this effect serves to modify the background upon which ECDI or MTSI modes develop, shifting growth/damping and re-routing energy from coherent waves to thermal motion. At sufficiently large excitation amplitudes, this pathway can dominate, producing strong spectral broadening, elevated temperatures and transport increases by orders of magnitude.
3. Set-up of the simulations performed
We describe here the set-up of the one-and two-dimensional simulations performed to investigate the plasma response under externally applied modulation. The simulation parameters adopted for each configuration are based on well-established benchmark cases and/or reference studies widely used in the E × B plasma community, ensuring that the baseline behaviour of both systems is well understood and directly comparable with existing literature. The one-dimensional (1-D) and 2-D set-ups are therefore assessed primarily on their own merit, rather than for direct quantitative cross-comparison, with the aim of examining modulation effects within representative and known discharge conditions.
3.1. Simulations’ set-up and conditions in the 1-D azimuthal configuration
The configuration of this case corresponds to the azimuthal coordinate of a Hall-thruster-representative discharge with the settings that closely follow those reported by Lafleur et al. (Reference Lafleur, Baalrud and Chabert2016a
). The radial, axial and azimuthal directions of the coordinate system for this simulation set-up are denoted by x, y and z, respectively. The domain’s size is
$L_{z}$
= 0.5 cm, which is discretised using
$N_{i}=100$
computational nodes. A uniformly distributed and stationary axial electric field (
$E_{y,0}$
) with the magnitude of
$20\, \mathrm{kV\,m}^{-1}$
and a constant radial magnetic (
$B_{x}$
) with an intensity of 20 mT are applied. The simulation is initialised by loading electron and ion particles with a uniform density of
$1\times 10^{17}\,\mathrm{m}^{-3}$
inside the domain. The initially loaded electrons and ions are sampled from Maxwellian distributions at the temperatures of 2 and 0.1 eV, respectively. The number of macroparticles per cell for either electron or ion species is 200. The time step is
$5\times 10^{-12}\, \mathrm{s}$
, and the plasma properties are averaged and recorded every 20 time steps. Noting that the dominant physics of interest along the azimuthal coordinate, namely the excitation, growth and saturation of the involved azimuthal instabilities, is primarily collisionless in nature, the collisions are neglected in the simulations.
To represent the axial convection of the instability waves and, hence, to limit their growth, a fictitious axial extent with the length of 1 cm is considered, which is consistent with the adopted value of Lafleur et al. (Reference Lafleur, Baalrud and Chabert2016a ). The electrons and ions that leave the axial boundary on one side are resampled from their initial distribution and reinjected into the domain from the opposite axial boundary with a randomly sampled azimuthal position. Within the finite axial extent of the domain, ions are allowed to accelerate under the imposed axial electric field. Along the azimuthal direction, a periodic boundary condition is applied on the particles, while along the radial coordinate, particles are allowed to move freely. Note that even though a finite axial length is assumed, the simulation remains a 1-D problem, as the Poisson’s equation is solely solved along the azimuthal coordinate. To guarantee an azimuthally periodic solution of the plasma potential, zero-value Dirichlet conditions are imposed at both azimuthal boundaries of the domain.
The described set-up represents the baseline simulation. To study the response of plasma to external modulation, a temporally oscillatory electric field is applied on top of the constant axial electric field (
$E_{y,0}$
). Therefore, in these ‘modulated’ simulations, the axial electric field at each time step is
$E_{y}=E_{y,0}+A_{M}\!\sin (2\pi f_{M}t)$
, where
$f_{M}$
is the modulation frequency in Hz and
$A_{M}$
is the modulation amplitude. The modulated simulations are initialised from the quasi-steady state of the baseline case at the end of
$10\,\mu \mathrm{s}$
and run for an additional 10
$\mu \mathrm{s}$
. Hence, the total simulated time is 20
$\mu \mathrm{s}$
.
The modulation amplitude is chosen to be
$50\, \%$
of the stationary axial electric field (
$A_{M}=0.5E_{y,0}=10\, \mathrm{kVm}^{-1}$
). Modulated simulations are performed for various modulation frequencies across the range of
$1{-}800\, \mathrm{MHz}$
. For each modulation frequency, simulations are repeated eight times (each for
$20\, \mu \mathrm{s}$
duration) to provide a statistically representative result.
The electrostatic 1-D PIC code used for the simulations has been developed by the authors (M.R. and F.F.), and verified and benchmarked in Faraji, Reza & Knoll (Reference Faraji, Reza and Knoll2022) and Reza, Faraji & Knoll (Reference Reza, Faraji and Knoll2022).
3.2. Simulations’ set-up and conditions in the 2-D radial-azimuthal configuration
The simulations’ configuration resembles the radial-azimuthal cross-section of a typical Hall thruster with the numerical and physical set-up similar to the conditions of the community benchmark problem reoprted by Villafana et al. (Reference Villafana2021). The computational domain represents a 2-D Cartesian plane with equal lengths of 1.28 cm along both simulation directions (
$L_{x}=L_{z}=1.28$
cm). In terms of axis notation, x, y and z denote respectively the radial, axial and azimuthal directions. The stationary applied axial electric field (
$E_{y,0}$
) has a magnitude of
$10\text{ kV}\mathrm{m}^{-1}$
and the external radial magnetic field intensity (
$B_{x}$
) is
$20\text{ mT}$
. The schematic of the simulation domain, and the directions of the imposed electric and magnetic fields relative to the coordinate system are illustrated in figure 1.

Figure 1. Schematic of the computational domain of the 2-D radial-azimuthal simulations, the coordinate system, and the applied stationary electric and magnetic fields.
The computational cell size is 50
$\mu \mathrm{m}$
, which corresponds to 256 nodes along both simulation dimensions. The time step is
$1.5\times 10^{-11}\, \mathrm{s}$
. The azimuthal electric field signal, which is used to analyse the instability spectra, is averaged and reported every 20 time steps.
Initially, electrons and ions are sampled from Maxwellian distribution functions at the temperatures of 10 and 0.5 eV, respectively, and are loaded uniformly in the domain with a density of
$1.5\times 10^{16}\, \mathrm{m}^{-3}$
. The initial number of macroparticles per cell is 100. The simulations are collisionless and the establishment of a steady state in the simulation set-up is achieved through introducing a particle injection source. The injection source is azimuthally uniform and follows a cosine profile along the radial direction, spanning from
$x$
= 0.09 cm to
$x$
= 1.19 cm, with the peak value of
$8.9 \times 10^{22}\, \mathrm{m}^{-3}\mathrm{s}^{-1}$
. The electron–ion pairs are sampled from Maxwellian distribution functions at the initial temperatures of the respective species and injected into the domain following the distribution prescribed by the injection source.
Regarding particle boundary conditions, any particle reaching the wall boundaries is eliminated and no secondary electron emission is considered. To reflect periodicity condition along the azimuthal coordinate, particles crossing the azimuthal boundaries are reinjected from the opposing boundary while retaining their velocity and radial position. As the simulations do not resolve the axial direction, a finite artificial extent of 1 cm is assumed along the
$y$
direction on both sides of the radial-azimuthal simulation plane (Villafana et al. Reference Villafana2021). Particles reaching the axial boundaries are resampled from the initial Maxwellian distributions and reloaded onto the simulation plane at the same radial and azimuthal positions. Ions are advanced under the axial electric field over the finite axial length assumed.
As for the boundary condition for the electric potential, a zero-volt Dirichlet condition mimics the grounded walls along the radial coordinate. A periodic condition is applied along the azimuthal direction.
The simulations are performed using a reduced-order PIC code developed by M.R. and F.F. The code uses a domain decomposition (Reza et al. Reference Reza, Faraji and Knoll2023a ) of 50 regions along both the radial and azimuthal directions. At this approximation level, the simulations have been demonstrated to faithfully reproduce the traditional full-2-D PIC results (Faraji, Reza & Knoll Reference Faraji, Reza and Knoll2023).
The baseline simulation features only the stationary axial electric field (
$E_{y}=E_{y,0}$
). Also, the total simulated time in the baseline condition is 30
$ \mu\text{s}$
. In the modulated simulations, the imposed axial electric field involves an oscillatory component. The total electric field is thus represented as
$E_{y}=E_{y,0}+A_{M}\sin (2\pi f_{M}t)$
. The state of the plasma at the end of the baseline simulation at
$30\, \mu\text{s}$
is used as the initial condition in the modulated simulations. We assessed three different modulation frequencies including ion plasma frequency (
$f_{pi}=5.8\, \mathrm{MHz}$
), electron cyclotron frequency (
$f_{ce}=560\, \mathrm{MHz}$
) and a mid-range frequency in between the previous two (40 MHz). Each frequency is assessed at four different modulation amplitudes of 5 %, 10 %, 25 % and 50 % of the stationary axial electric field’s magnitude. We performed three 15
$\mu$
-long simulation repetitions (from 30 to 45
$\mu \mathrm{s}$
) for every modulation condition to ensure statistical significance of the observed behaviours.
3.3. Considerations and notes regarding the simulation configurations adopted for this study
The configurations considered in this work are best viewed as a local slab idealisation of an E × B discharge: a representative cross-section of the plasma in which mean properties vary weakly in time compared with those in a Hall thruster operating under strongly time-dependent conditions such as breathing cycles. This idealisation is intentionally adopted to isolate the coupling between externally imposed axial electric-field modulation and azimuthal instability dynamics and transport. Within this framework, modulation-induced changes in fluctuation spectra, nonlinear coupling signatures and the transport-driving correlation
$\langle \tilde{n}_{e} \tilde{E}_{z}\rangle$
can be examined without conflation with global axial discharge evolution and strong axial gradients inherent to a full E × B discharge configuration.
As a direct consequence of this modelling choice and the absence of axial resolution, the axial electric field
$E_{y}$
is prescribed rather than evolved self-consistently. Moreover, coupling between low-frequency axial oscillations (including the breathing mode and the ion transit-time modes) and the high-frequency azimuthal instabilities studied here are not captured. A fully self-consistent treatment would require resolving axial plasma dynamics alongside transverse instability evolution, ultimately leading to a three-dimensional model that lies outside the scope of the present study. The results reported here should therefore be interpreted as idealised, mechanism-focused responses to field modulation, rather than as predictive of global axial field evolution or discharge self-organisation in a real Hall thruster.
Despite these idealisations, the qualitative mechanisms identified in this study – namely, spectral redistribution of instability energy, modulation-induced modification of nonlinear coupling pathways and the resulting changes in instability-driven electron cross-field transport – are expected to remain relevant in more complete, self-consistent three-dimensional configurations representative of Hall thrusters. Exact quantitative assessment of the observed trends, however, will require self-consistent inclusion of axial gradients and field evolution, collisional processes, time-dependent discharge dynamics (e.g.) the breathing mode (Boeuf & Garrigues Reference Boeuf and Garrigues1998; Barral & Ahedo Reference Barral and Ahedo2009; Sekerak et al. Reference Sekerak, Longmier, Gallimore, Brown, Hofer and Polk2013), and additional physics such as secondary electron emission (SEE) (Hobbs & Wesson Reference Hobbs and Wesson1967; Sydorenko et al. Reference Sydorenko, Smolyakov, Kaganovich and Raitses2005; Raitses et al. Reference Raitses, Smirnov, Staack and Fisch2006) and the associated near-wall conductivity (Morozov Reference Morozov1968; Ivanov, Ivanov & Bacal Reference Ivanov, Ivanov and Bacal2002; Barral et al. Reference Barral, Makowski, Peradzyński, Gascon and Dudeck2003).
Accordingly, the present work should be regarded as a preliminary feasibility investigation of instability and transport control in simplified E × B geometries, providing a controlled framework for isolating fundamental modulation effects prior to incorporating the full complexity of realistic E × B discharges. More specific remarks on the configuration assumptions and/or physics idealisations for the purposes of this study follow.
-
(i) Finite axial extent and particle reinjection: both the 1-D and 2-D configurations employ a finite fictitious axial length combined with particle reinjection. This approach is known to influence instability saturation and absolute transport levels by affecting particle energy balance, wave convection and by introducing an effective artificial collisionality, as demonstrated in previous studies (Lafleur et al. Reference Lafleur, Baalrud and Chabert2016a ; Asadi, Taccogna & Sharifian Reference Asadi, Taccogna and Sharifian2019; Tavant Reference Tavant2019; Villafana et al. Reference Villafana2021). In the present work, this treatment is not intended to reproduce physical axial kinetics, but rather to provide a controlled and widely used means of limiting instability growth and maintaining a statistically stationary turbulent state in the absence of an explicitly resolved axial coordinate.
Crucially, the reinjection strategy and axial length are held strictly identical across all baseline and modulated simulations. As a result, the reported changes in instability spectra and electron transport reflect relative modulation-induced effects, rather than differences arising from axial closure or boundary treatment. Fully eliminating such artificial effects would ultimately require self-consistent resolution of the axial coordinate, leading to a fully three-dimensional treatment, outside the scope of the current study.
-
(ii) Finite azimuthal extent: a finite azimuthal domain inherently limits the smallest accessible azimuthal wavenumbers
$k_{z}$
and can therefore influence the development of long-wavelength modes. This limitation is most relevant in regimes where energy transfer induced by modulation is directed towards wavelengths comparable to the domain extent; in such cases, increasing the azimuthal length
$L_{z}$
may permit additional long-wavelength modes, potentially altering the saturated spectra and associated transport. Larger azimuthal domains can also enable the formation of global coherent structures (Koshkarov et al. Reference Koshkarov, Smolyakov, Raitses and Kaganovich2019), such as the rotating spoke (Raitses et al. Reference Raitses2012; Powis et al. Reference Powis, Carlsson, Kaganovich, Raitses and Smolyakov2018).In the present study, however, the emphasis is on high-frequency, short-wavelength instabilities, primarily the ECDI and MTSI, and on how external modulation modifies their spectra and transport impact, rather than on global large-scale azimuthal dynamics. To ensure adequate spectral representation, the azimuthal domain in both 1-D and 2-D configurations contains multiple wavelengths of the dominant modes. In all cases examined, the dominant spectral modes remain well below half of the domain size, indicating that the key instability dynamics analysed here are not constrained by the finite azimuthal extent.
-
(iii) Wall interactions, SEE and collisions: in the 2-D radial-azimuthal configuration, radial boundaries impose grounded-wall (Dirichlet) conditions for the electrostatic potential and wall interactions are simplified to particle absorption without SEE. In a Hall thruster, channel walls are typically dielectric and reach a floating potential. The resulting sheath dynamics and SEE can modify near-wall electric fields, electron energy balance and effective cross-field mobility through near-wall conductivity.
SEE and near-wall processes introduce additional energy-loss and redistribution channels. Likewise, physical collisions provide momentum and energy relaxation, and can modify both electron mobility and heating–cooling balance. Omission of SEE and collisions can therefore lead to different absolute electron temperature levels compared with configurations that include realistic wall interactions. Furthermore, collisions and wall losses can alter instability regimes and saturation pathways by modifying electron temperature, effective dissipation and confinement. Consequently, absolute fluctuation amplitudes and saturation levels may differ when collisions and SEE are included.
Prior unmodulated radial-azimuthal kinetic simulations by Reza et al. (Reference Reza, Faraji and Knoll2023c ) demonstrate that increasing effective SEE produces a pronounced cooling-and-confinement trend: stronger SEE cools bulk electrons while increasing plasma density, consistent with cold secondaries modifying sheath drops and wall-loss balance. That study further shows that SEE alters the relative amplitude and content of the fluctuation spectrum. At moderate SEE, the familiar ECDI/MTSI harmonic structure persists whereas, at higher SEE, a large-scale coherent long-wavelength mode can emerge as dominant, with the ECDI signature becoming more broadband and the MTSI contribution weakening or disappearing. These spectral changes are accompanied by corresponding changes in electron mobility, including strengthened near-wall contributions and non-monotonic trends as dominant modes shift. The electron energy distribution also narrows with increasing SEE, reflecting progressive depletion of the high-energy tail due to net cooling by low-energy secondaries.
These established SEE-driven modifications motivate a careful interpretation of the present modulation results. Because the modulation response observed here is mediated by spectral redistribution and changes in phase and correlation structure (e.g. via
$\langle \tilde{n}_{e} \tilde{E}_{z}\rangle$
), any physics that materially changes the baseline dominant instability content and its spectral amplitude can be expected to renormalise quantitative transport levels and potentially shift the parameter window in which modulation most effectively suppresses or enhances transport. For example, the strongest suppression identified in the present 2-D configuration occurs in a mid-frequency regime where modulation weakens the ECDI peak, whereas cyclotron-resonant forcing at sufficiently high amplitudes produces strong heating and transport enhancement. Since stronger SEE tends to cool electrons and modify near-wall conductivity pathways, it may alter the heating balance and the relative importance of near-wall versus bulk contributions to net mobility, thereby affecting both the magnitude and optimal tuning of modulation-based transport control.At the same time, the central transport-control mechanism identified in this work – modulation-induced reshaping of instability spectra and phase correlations – is expected to remain qualitatively relevant in the plasma bulk when SEE is included, even though quantitative transport levels and optimal modulation parameters may change.
-
(iv) Geometric curvature: the 2-D radial-azimuthal configuration is modelled on a planar Cartesian domain. As a result, geometric curvature terms associated with a cylindrical
$(r,\theta )$
formulation are neglected. This ‘straightened’ geometry is commonly adopted in benchmark and local PIC studies of E × B microinstabilities, such as the ECDI and MTSI, for which characteristic wavelengths are small compared with the device radius. For the high-frequency, short-wavelength instabilities that dominate the present study, the planar Cartesian approximation is therefore overall sufficient and curvature is not expected to significantly alter the instability mechanisms or modulation impacts reported here.Curvature effects become more important for larger-scale azimuthal structures and global low-mode-number dynamics, such as rotating spokes, as well as for full-annulus simulations. These effects are not captured in the present configuration and serve as a natural extension to the current work.
4. Results
In this section, we discuss the impact of electrostatic modulation on the frequency and wavenumber spectra of the azimuthal instabilities and, in turn, on the axial electron mobility induced by the excited instabilities. We first present the outcomes for the 1-D azimuthal case and then proceed to the results in the 2-D radial-azimuthal setting. The 2-D results enable us to observe the effects of electrostatic modulation on the plasma response in the presence of radial physics.
4.1. 1-D azimuthal configuration
Figure 2 provides a comparison between the frequency spectra of the instabilities in the baseline and the modulated simulations at various modulation frequencies. For the spectral analysis, temporal fast Fourier transform (FFT) is applied to the azimuthal electric field (
$E_{z}$
) signal over 2–10
$ \mu \mathrm{s}$
.

Figure 2. Variation of the frequency spectrum relative to the baseline case from 1-D azimuthal simulations with various modulation frequencies. In each case, the spectrum represents the average of the temporal FFT of azimuthal electric field (
$E_{z}$
) signal over all azimuthal positions and over eight simulation repetitions. The characteristic frequencies displayed in grey dashed lines from left to right correspond to theoretical ion acoustic frequency (
$f_{IA}$
), ion plasma frequency (
$f_{pi}$
), Hall circulation frequency (
$f_{E}=E_{y,0}/B_{x}L_{z}$
), and first and second harmonics of electron cyclotron frequency (
$f_{ce},2f_{ce}$
).
This duration refers to the time interval after arriving at quasi-steady state at 10
$\mu\text{s}$
, hence, the absolute time of 12–20
$\mu \mathrm{s}$
. The same is true whenever referring to the time interval in this section. The spectra for each simulation scenario shown in figure 2 is the average of the FFTs of
$E_{z}$
across all azimuthal positions and over eight simulation repetitions.
In addition, the azimuthal wavenumber (
$k_{z}$
) spectra of the instabilities in the modulated simulations are compared against the baseline
$k_{z}$
spectrum in figure 3. The
$k_{z}$
spectrum in each case is obtained by calculating the spatial FFT of the
$E_{z}$
signal across the entire domain at a specific time instance between 2 and 10
$\mu \text{s}$
. The plots then represent the average of the spatial FFTs over the full 8-
$\mu\text{s}$
duration and across eight simulation repetitions.

Figure 3. Variation of the azimuthal wavenumber (
$k_{z}$
) spectrum relative to the baseline case from 1-D azimuthal simulations with various modulation frequencies. In each case, the spectrum represents the temporally averaged (over 2–10
$\mu\text{s}$
) spatial FFT of azimuthal electric field (
$E_{z}$
) signal averaged over eight simulation repetitions. Panels (a)–(c) denote varying frequency ranges.
In the baseline 1-D set-up, the ECDI establishes within the simulation as the dominant instability with the frequency of
$4.5\, \mathrm{MHz}$
(figure 2) and azimuthal wavenumber of
$619\,\, \mathrm{m}^{-1}$
(figure 3) – wavelength of approximately
$1.6$
mm.
The plots in figures 2 and 3 show that electrostatic modulation significantly changes the spectral amplitude of the instabilities. The specific impact varies depending on the frequency of modulation (
$f_{M}$
). At very low frequency (
$f_{M}=1\, \mathrm{MHz}$
), the ECDI’s peak in the frequency spectrum is flattened and its amplitude is nearly halved compared with its original value in the baseline.
Slightly increasing the modulation frequency (cases with
$f_{M}=4.5$
and
$5.8\, \mathrm{MHz}$
) results in a shift of the peak of the frequency spectra towards the frequency of the ion acoustic instability (
$f_{IA}\approx 3.35\, \mathrm{MHz}$
, as calculated from
$f_{IA}\approx ({k_{z}\lambda _{D}\omega _{pi}}/{(2\pi \sqrt{1+k^{2}\lambda _{D}^{2})}})$
). From figure 2(a), the peak wavenumber also migrates towards
$k_{z}=412\,{\rm m}^{-1}$
corresponding to a longer wavelength of approximately 2.4 mm. In these cases, the spectral amplitude of an extended range of frequencies across
$15{-}250\, \mathrm{MHz}$
is increased (figure 3).
Figure 2 further shows that modulation at
$20{-}70\, \mathrm{MHz}$
excites instabilities with frequencies in the range of
$\geq f_{M}$
and
$\leq f_{ce}=560\, \mathrm{MHz}$
, with the distribution of energy being mostly concentrated on the discrete peaks. At modulation frequencies of
$f_{M}=150$
and
$500\, \mathrm{MHz}$
, a significant reduction in the amplitude of ECDI is evident. Meanwhile, the oscillations within the low-frequency range of the spectra (
$f\lt f_{IA}$
) and longer wavelengths are amplified. In these cases, the amplitudes of the high frequency range of the spectra experience a slight increase. Interestingly, these variations are not observed in the spectral amplitudes in the case with
$ f_{M}=300$
MHz.
Modulation at higher frequencies near the electron cyclotron frequency (
$f_{M}=f_{ce}=560\, \mathrm{MHz}$
) intensifies the high-frequency range of the spectrum while having similar effects to cases with
$f_{M}=150$
and
$500\, \mathrm{MHz}$
on the frequency range below
$30\, \mathrm{MHz}$
. These variations in the spectral amplitude of the instabilities appear to diminish at higher modulation frequencies such as
$f_{M}=600\, \mathrm{MHz}$
and nearly disappear when the modulation is applied at
$f_{M}=800\, \mathrm{MHz}$
.
As an interesting observation, in all cases with modulation frequency
$f_{M}\leq 70\, \mathrm{MHz}$
(and at
$f_{M}=300$
MHz), a peak between 1.8 and 2.2 MHz appears in their frequency spectrum. Whereas in cases with
$150\, \mathrm{MHz}\leq f_{M}\leq 600\, \mathrm{MHz}$
, the entire range of low-frequency spectrum below the ion acoustic frequency is enhanced.
Looking at figure 3, it is evident that modulation in most cases tend to shift the
$k_{z}$
spectra towards smaller wavenumbers, effectively increasing the wavelength of the dominant instabilities. The consistent shift of the dominant azimuthal mode – the ECDI – towards smaller
$k_{z}$
under modulation is consistent with how a time-dependent E × B drift modifies mode excitation in drift-cyclotron-type instabilities.
In the present configuration, the modulation is imposed on the cross-field electric field, so the electron drift becomes time-dependent,
and the linear resonance/phase-matching underlying the ECDI-like response is controlled by the Doppler-shifted frequency
$\omega -k_{z}V_{E}$
(cyclotron-harmonic resonance in kinetic descriptions) (Lafleur et al. Reference Lafleur, Baalrud and Chabert2017). When
$V_{E}$
oscillates, each mode experiences a periodic detuning
$\delta \omega (t)\sim k_{z}\delta V_{E}\!\sin (\omega _{M}t)$
. This detuning grows linearly with
$k_{z}$
, meaning that larger
$k_{z}$
modes are more strongly dephased/detuned over a modulation cycle. In practice, that reduces their net amplification (or broadens their spectral peak) because they spend a larger fraction of time away from the instantaneous resonance condition and/or accrue a larger phase mismatch. Conversely, lower-
$k_{z}$
modes have smaller
$\delta \omega$
and therefore remain closer to resonance for a greater fraction of the modulation period, giving them a higher effective growth relative to high-
$k_{z}$
modes. As a result, the peak of the
$k_{z}$
-spectrum shifts to smaller
$k_{z}$
.
The variations in the spectral amplitude of the azimuthal instabilities will inherently translate into change in the electron cross-field (axial) transport caused by these fluctuations. To quantify the change in electron transport, the electrons’ axial current density (
$J$
) and axial mobility (
$\mu$
) in the presence of electrostatic modulation at various frequencies are compared in figure 4 against their values in the baseline simulation. These quantities are obtained from the spatiotemporal signals of the electron number density (
$n_{e}$
) and the electrons’ axial drift velocity (
$V_{e,y}$
) from the simulations according to the following relations:
\begin{align} J &= {\frac{1}{N}}\sum\limits_{k = 1}^N {{\frac{e}{{L_z}({t_2} - {t_1})}}} \int_{{t_1}}^{{t_2}} {\int_0^{{L_z}} } n_e^k\left( {z,t} \right)V_{e,y}^k(z,t){\mkern 1mu} {\rm{d}}z{\mkern 1mu} {\rm{d}}t,\\[-10pt]\nonumber \end{align}
\begin{align} \mu &= {\frac{1}{N}}\sum\limits_{k = 1}^N {{\frac{1}{{L_z}({t_2} - {t_1}){E_{y,0}}}}} \int_{{t_1}}^{{t_2}} {\int_0^{{L_z}} } V_{e,y}^k(z,t){\mkern 1mu} {\rm{d}}z{\mkern 1mu} {\rm{d}}t, \end{align}

Figure 4. Electrons’ axial current density (
$J$
) and axial mobility (
$\mu$
) from the 1-D azimuthal simulations with various modulation frequencies. The colour bars represent the spatiotemporal mean value (over
$2{-}10\, \mu\text{s}$
and entire domain) over eight simulation repetitions. Panels (a) and (c) display the absolute values of
$J$
and
$\mu$
, and panels (b) and (d) show the normalised change in
$J$
and
$\mu$
in the modulated simulations with respect to the corresponding quantities in the baseline case. The error bars indicate one standard deviation among the eight simulation repetitions in each case. The solid and the dotted black lines respectively represent the average and one standard deviation of
$J$
and
$\mu$
corresponding to eight simulation repetitions for the baseline case.
where
$t_{1}=2 \mu\text{s}$
(12
$\mu\text{s}$
in absolute time) and
$t_{2}=10\, \mu\text{s}$
(20
$\mu\text{s}$
in absolute time), and
$N$
is the total number of repetition instances for each modulation scenario, which in this case is 8. Where (4.2) and (4.3) are used to derive
$J$
and
$\mu$
from the 2-D simulations’ results in § 4.2, the spatial integration has been carried out over the full radial-azimuthal cross-section.
To determine the statistical significance of the observed variations in the electron axial current and mobility, the first standard deviation among the eight-simulation set in each modulation case is indicated in figure 4.
The normalised variation percentages of the electrons’ axial current density and axial mobility with respect to their baseline values – also shown in figure 4 – are calculated according to
\begin{align}|\overline {\Delta J} | = {\frac{1}{N}}\sum\limits_{k = 1}^N {\left| {{\frac{J_M^k - J_B^k}{J_B^k}}} \right|} ,\quad |\overline {\Delta \mu } | = {\frac{1}{N}}\sum\limits_{k = 1}^N {\left| {{\frac{\mu _M^k - \mu _B^k}{\mu _B^k}}} \right|}. \end{align}
In the relations in (4.4), the subscripts ‘
$M$
’ and ‘
$B$
’ denote the quantities in the modulated and baseline simulations, respectively.
The plots in figure 4 show that, overall, modulation at frequencies below
$150\, \mathrm{MHz}$
leads to a reduction in the electrons’ axial current density and axial mobility by an average of approximately 15 %–30 % relative to the baseline case. The maximum reduction occurs at
$f_{M}=40\, \mathrm{MHz}$
, with approximately
$29\,\%$
reduction. The decrease in electron transport is attributed to two factors: (i) diminishing amplitude of the ECDI mode as a major contributor to the axial transport of electrons; (ii) reduced correlation between azimuthal density and electric field fluctuations (
$\langle \tilde{n}_{e}\tilde{E}_{z}\rangle$
) – as will be further discussed in § 5.1.
At frequencies
$\geq$
150 MHz, modulation either enhances axial electron transport or results in no significant change relative to the baseline. Modulation at 150 and 500 MHz increases the electron current by approximately 26 % and 21 %, respectively, while increasing the electrons’ axial mobility by approximately 13 % in both cases. When exciting at the electron cyclotron frequency (
$f_{M}=f_{ce}=560\, \mathrm{MHz}$
), the electrons’ axial current and mobility increase by approximately 12 % and 8 %, respectively.
In these cases, the enhanced transport is due to the combined effect of the following processes: (i) increased correlation between the electron density and azimuthal electric field fluctuation
$\langle \tilde{n}_{e}\tilde{E}_{z}\rangle$
(§ 5.1); (ii) amplification of the waves close to electron cyclotron frequency, which can resonate with the electrons’ cyclotron motion and increase the electron Larmor radius, thereby weakening magnetic confinement; (iii) amplification of the longer wavelength modes that also contribute to transport.
The variations in both the electrons’ current and mobility for the remaining modulation frequencies, namely
$f_{M}=300, 600$
and
$800\, \mathrm{MHz}$
, are not statistically significant because the changes in these cases fall within the range of variations observed among different simulation runs with the baseline condition.
4.2. 2-D radial-azimuthal configuration
In the simulated radial-azimuthal configuration, the plasma dynamics is predominantly influenced by the existence of two instabilities, namely the ECDI and the MTSI. The MTSI is reported to have a radial wave vector component (Janhunen et al. Reference Janhunen, Smolyakov, Chapurin, Sydorenko, Kaganovich and Raitses2018c ; Petronio et al. Reference Petronio, Tavant, Charoy, Alvarez-Laguna, Bourdon and Chabert2021), which prevents it from being captured in purely azimuthal simulations. Its presence notably modifies the fluctuation spectra, significantly influencing the coupling of modulation waves to the plasma and the energy cascade process.
Figures 5 and 6 present the variations in the frequency and wavenumber spectra compared with the baseline case for various modulation frequencies and amplitudes.

Figure 5. Variation of the frequency spectrum relative to the baseline case from the radial-azimuthal simulations with various modulation frequencies and amplitudes. In each case, the spectrum represents the average of the temporal FFT of azimuthal electric field (
$E_{z}$
) signal over all azimuthal and radial positions and over three simulation repetitions. The characteristic frequencies displayed in grey dashed lines from left to right correspond to theoretical ion acoustic frequency (
$f_{IA}$
), ion plasma frequency (
$f_{pi}$
), Hall circulation frequency (
$f_{E}=E_{y,0}/B_{x}L_{z}$
), and first and second harmonics of electron cyclotron frequency (
$f_{ce},2f_{ce}$
).

Figure 6. Variation of the azimuthal wavenumber (
$k_{z}$
) spectrum relative to the baseline case from the radial-azimuthal simulations with various modulation frequencies and amplitudes. In each case, the spectrum represents the spatiotemporally averaged [over
$5{-}15\, \mu\text{s}$
(35–45
$\mu\text{s}$
absolute time) and all radial positions] spatial FFT of azimuthal electric field (
$E_{z}$
) signal averaged over three simulation repetitions.
In baseline conditions (without modulation, solid black lines), the wavenumbers of the ECDI and the MTSI modes that develop in the simulation are
$1186$
and
$237\,\, \mathrm{m}^{-1}$
, respectively. These values will be referenced in the follow-on discussions surrounding figures 5 and 6. They also align well with the theoretical predictions from (2.1) and (2.2), which yield wavenumbers of
$1114$
and
$207\, \mathrm{m}^{-1}$
for the ECDI and the MTSI, respectively.
The different modulation frequencies examined in figures 5 and 6 (5.8, 40 and 560 MHz) exhibit distinct effects on the plasma instability spectra, with varying impacts based on the amplitude of modulation. Modulation at the plasma frequency (5.8 MHz) primarily excites high-frequency oscillations within the 10–200 MHz range, particularly at higher modulation amplitudes (
$A_{M}\geq 0.25E_{y}$
). At low modulation amplitudes (e.g.
$A_{M}=0.05E_{y}$
), the frequency spectrum closely resembles the baseline case, with minor deviations. However, as the modulation amplitude increases, spectral energy spreads to higher frequencies beyond the modulation frequency, activating modes that were less prominent in the baseline.
Modulation at 40 MHz introduces waves at the modulation frequency itself and excites higher harmonics as the modulation amplitude increases. At an intermediate modulation amplitude, particularly
$A_{M}=0.25E_{y}$
, energy is notably redistributed, with reduced spectral power in the mid-frequency range (2–20 MHz). This frequency range overlaps with the typical ECDI modes, leading to a significant suppression of this instability.
Modulation at the electron cyclotron frequency (560 MHz) and sufficiently high amplitudes (
$A_{M}\geq 0.25E_{y}$
) induces substantial changes in the instability spectrum, significantly altering plasma dynamics. At these amplitudes, 560-MHz modulation effectively disrupts the ECDI’s wave structure, as reflected in the reduced spectral power in the ECDI frequency range. Concurrently, high-frequency oscillations are strongly enhanced, with pronounced spectral energy around the modulation frequency and its harmonics. Moreover, remarkable growth of lower frequencies and larger spatial structures is observed in both the frequency (figure 5) and wavenumber spectra (figure 6). This behaviour is indicative of an enhanced coupling across scales, where high-frequency modulation affects the large-scale plasma organisation.
From figure 6, a noteworthy observation is that at modulation frequencies of 5.8 and 40 MHz, the wavenumber of the ECDI shifts towards lower values, with the extent of this shift increasing as the modulation amplitude rises. In contrast, the wavenumber of the MTSI remains unaffected across these modulation conditions. These observations can be justified following the same logic presented with regards to figure 3 for the 1-D configuration – that the shift of the ECDI mode towards smaller wavenumbers under modulation is consistent with the modification induced by time-dependent E × B drift in the excitation of drift-cyclotron-type instability modes.
Moving on to figure 7, where variations in the electrons’ axial current density and mobility with modulation frequency and amplitude are presented, it is notably seen that modulation at
$f_{M}=40\,\text{MHz}$
and
$A_{M}=0.25 E_{y}$
results in the largest reductions relative to the baseline condition – approximately 38 % in electron current and approximately 30 % in electron mobility. While noting the difference in the simulations’ set-up, these reductions are consistent with the trend seen in the purely azimuthal studies (figure 4), where the greatest reduction in electron transport – by approximately 30 % – also occurred at this modulation frequency.

Figure 7. Electrons’ axial current density (
$J$
) and electrons’ axial mobility (
$\mu$
) from the radial-azimuthal simulations with various modulation frequencies and amplitudes. The colour bars represent the spatiotemporal mean value [over
$5{-}15\,\mu\text{s}$
(35–45
$\mu\text{s}$
absolute time) and entire domain] averaged over three simulation repetitions. Panels (a) and (c) display the absolute values of
$J$
and
$\mu$
, and panels (b) and (d) show the normalised change of
$J$
and
$\mu$
in the modulated simulations with respect to the respective quantities in the baseline case. The error bars indicate one standard deviation among the three simulation repetitions in each case. The solid and the dotted black lines represent the average and one standard deviation corresponding to three simulation repetitions for the baseline case.
Modulation at the ion plasma frequency (5.8 MHz) also reduces electrons’ current density and mobility by up to approximately 24 % at modulation amplitudes of
$A_{M}=0.1E_{y}$
and
$0.25 E_{y}$
.
In contrast to both above-mentioned observations, modulation at the electron cyclotron frequency (560 MHz) with a low amplitude (
$A_{F}=0.05E_{y}$
) significantly decreases electron current and mobility by 37 % and 24 %, respectively. Yet, as the modulation amplitude increases, the electron axial transport is drastically enhanced, rising by more than two orders of magnitude at
$A_{M}=0.25E_{y}$
and
$0.5 E_{y}$
. This level of enhancement in electron transport can be attributed to significant heating of particles under these modulation conditions, leading to a breakdown of electrons’ magnetic confinement due to elevated electron temperatures (substantially elongated electron cyclotron orbit along the axial direction).
To support the previous explanation, Appendix A presents the overall impact of electrostatic modulation – and, consequently, the altered effect of instabilities – on the particles’ velocity distribution functions and the time-averaged plasma profiles. The plots in Appendix A notably show that modulation at 560 MHz with
$A_{M}\geq 0.1E_{y}$
results in considerable heating of the electrons and ions, and a substantial broadening of their velocity distributions along both the radial and azimuthal directions.
5. Discussions
5.1. Deeper analysis of wave interactions
Bispectral analysis (Kim & Powers Reference Kim and Powers1979) is a standard method to identify the nonlinear interactions among spectral components. The bispectrum quantifies phase coupling between three frequency components that satisfy the resonance condition
$f_{1}+f_{2}=f_{3}$
, revealing the presence of coherent three-wave interactions. To isolate the strength of this coupling independent of signal amplitude, the bispectrum is normalised to obtain the bicoherence, which ranges from 0 (no phase correlation) to 1 (perfect phase locking). Regions of high bicoherence indicate frequency triads (Newell, Nazarenko & Biven Reference Newell, Nazarenko and Biven2001) that are phase-coherently coupled, providing insight into the underlying nonlinear energy transfer mechanisms.
In this section, we followed the approach described by Ritz et al. (Reference Ritz1988) to carry out bicoherence analysis on our simulation data. The code implementation available on GitHub (Przybocki Reference Przybocki2025) was used.
The bicoherence analysis was performed on spatiotemporal azimuthal electric field fluctuations data from the 1-D azimuthal simulation cases acquired at a sampling frequency (
$f_{s}$
) of 10 GHz. The total duration analysed was 2–10
$\mu \text{s}$
, yielding 80 000 time samples per spatial location. Spatial data were taken from two arbitrary located points in the simulation domain. The results were found to be largely insensitive to the specific choice of locations or their relative separation.
The bispectral computation used a Hann window with 50 % overlap and an FFT length of 2000, resulting in a frequency resolution of
$\Delta f={f_{s}}/({n_{fft}})=5$
MHz. The bicoherence was then computed by averaging over 80 overlapping blocks of data. These parameters were chosen to balance spectral resolution with statistical representativeness. However, the limited frequency resolution imposed by the FFT length means that interactions involving low-frequency components are not well resolved in the bicoherence analysis.
The salient results of the analysis performed are described in the following.
In the unmodulated baseline case, the bicoherence map in figure 8 shows that coherent energy transfer is confined to discrete, localised bright spots arranged on a grid-like pattern, spaced approximately 152 MHz apart. These correspond to interactions of a 152-MHz wave with itself and its primary harmonics, generating higher-order harmonics. The nature of this wave is not known as it does not coincide with the characteristic frequencies of the system. However, important to our discussion is that it provides a comparative reference for the bicoherence maps from the modulated cases.

Figure 8. Bicoherence map of azimuthal electric field fluctuations from the baseline (unmodulated) 1-D azimuthal simulation.
In the presence of modulation with different frequencies, the bicoherence maps in figure 9 consistently show a grid-like pattern of bright spots spaced approximately at the modulation frequencies, particularly in cases with
$f_{M}=20, 40, 70$
and
$150$
MHz.

Figure 9. Bicoherence maps of azimuthal electric field fluctuations from the 1-D azimuthal simulations with various modulation frequencies.
To explain and understand these patterns, we refer to figure 10, which shows zoomed-in views of the frequency spectra from the FFT of the azimuthal electric field signal previously presented in figure 2. From the plots in figure 10, we observe that the application of excitation to plasma does not manifest as a single sharp spectral line at the modulation frequency. Instead, two sidebands appear equally spaced around the modulation frequency (although not necessarily with comparable amplitudes), while the amplitude at the exact excitation frequency is notably suppressed. This spectral splitting can be physically interpreted as the result of three-wave coupling between the modulating wave (with frequency
$f_{M}$
) and a newly generated lower-frequency mode (with frequency
$f_{LF}$
). The latter mode refers to the spectral peak near 2 MHz observed in the FFT spectra (figure 2), which appears across most modulation cases. Such modulation leads to the redistribution of spectral power from the central frequency (
$f_{M}$
) into sidebands located at
$f_{M}\pm f_{LF}$
. This behaviour is noticeable in the FFT plots of figure 10 for modulation frequencies
$f_{M}\geq 20$
MHz
$.$
From the full FFT plots in figure 2, the low-frequency mode excited and coupled with the modulating wave is found to be around the frequency (
$f_{LF}$
) of
$2$
MHz.

Figure 10. Zoomed-in views on the plots in figure 2 for modulation frequencies in the range of 20–800 MHz; average temporal FFTs of the azimuthal electric field (
$E_{z}$
) signal from the 1-D azimuthal simulations over all spatial locations and eight simulation repetitions.
The FFT plots of figure 2 also showed for most modulation cases – particularly with
$f_{M}\geq 20$
MHz – the excitation of a series of peaks at frequencies corresponding to integer multiples of the excited sidebands (
$n(f_{M}\pm f_{LF}),$
with
$n=1,2, \ldots$
) around
$f_{M}$
, suggesting the excitation of harmonics of the sidebands of modulation frequency as well.
In line with these interpretations from the FFT plots, the bright spots in the bicoherence maps in figure 9, especially visible for
$f_{M}=20, 40, 70$
and
$150$
MHz, indicate strong phase coupling between the modulation frequency and either of the sideband components, satisfying the three-wave nonlinear interactions and leading to the excitation of higher-order harmonics. This confirms that the observed high-frequency FFT peaks in figure 2 beyond
$f_{M}$
represent higher-order harmonics of the sidebands of the modulation frequency rather than being independent spectral modes.
For lower modulation frequencies of
$f_{M}=1, 4.5$
and
$5.8$
MHz, the spectral response in figure 2 exhibited a broad-band rise in the mid-frequency content. This coincides with extended bicoherence features (diagonal bright ridge around
$f_{1}\approx f_{2}$
for
$f_{M}=1$
MHz and extended bright areas for
$f_{M}=4.5, 5.8$
MHz).
These patterns in the bicoherence maps indicate the presence of coherent three-wave energy transfer over an extended frequency range, where frequency components within the mid-band are phase-locked and engage in nonlinear interactions. These couplings redistribute energy into higher frequencies through sum interactions (
${f_{3}}=f_{1}+f_{2}$
) and into lower frequencies through difference interactions
$(f_{3}=|f_{1}-f_{2}| )$
over extended frequencies. This explains the observed spectral broadening in the FFTs of figure 2 at low-frequency modulation.
Similar spectral behaviours were observed from the FFT plots of the 2-D simulations (figure 5), especially evident in the
$f_{M}=40\,\mathrm{MHz}$
case, where the excitation was again seen not to produce a single peak, but instead result in the appearance of a series of sideband peaks around the modulation frequency. The FFT spectra in figure 5 also exhibited harmonics of the modulation sidebands, suggesting similar underlying nonlinear coupling mechanisms consistent with the 1-D azimuthal cases analysed and discussed here.
5.2. Closer assessment of variations in instability-induced electron transport
As discussed in § 2.1, instability-induced axial electron transport in E × B discharges is governed by the correlation between density and azimuthal electric-field fluctuations
$\langle \tilde{n}_{e}\tilde{E}_{z}\rangle$
. By changing the spectral characteristics and influencing both the amplitude of dominant azimuthal instabilities as well as the relative phase between azimuthal fluctuations in the electron number density
$\tilde{n}_{e}$
and azimuthal electric field
$\tilde{E}_{z}$
, electrostatic modulation can modify this correlation and thus the effective axial current and mobility.
To demonstrate this, figure 11 shows the normalised (with respect to averages of
$\tilde{n}_{e}$
and
$| \tilde{E}_{z}|$
) correlation term
$\langle \tilde{n}_{e}\tilde{E}_{z}\rangle$
along with the instability-induced axial electron mobility – obtained from (2.3) and averaged over time and across simulation repetitions – for various excitation frequencies. The plots shown are for the 1-D azimuthal simulations as demonstrative examples.

Figure 11. Variation of (a) instability-induced electrons mobility (
$\mu _{inst}$
) and (b) the normalised fluctuations correlation term
$\langle \tilde{n}_{e}\tilde{E}_{z}\rangle$
, from the 1-D azimuthal simulations with various modulation frequencies. The colour bars represent the spatiotemporal mean value (over
$2{-}10\, \mu\text{s}$
and entire domain) over eight simulation repetitions. The error bars indicate one standard deviation among the eight simulation repetitions in each case. The solid and the dotted black lines respectively represent the average and one standard deviation of
$\mu _{inst}$
and normalised
$\langle \tilde{n}_{e}\tilde{E}_{z}\rangle$
corresponding to eight simulation repetitions for the baseline case.
We note two observations from panels (a) and (b) in figure 11: (i) in the transport suppression regime (low-to mid-frequency modulation), phase alignment between
$\tilde{n}_{e}\,$
and
$\tilde{E}_{z}$
diminishes, lowering the induced axial electron mobility; (ii) in the transport enhancement regime (higher frequency modulation, especially near electron cyclotron frequency), the correlation between
$\tilde{n}_{e}\,$
and
$\tilde{E}_{z}$
increases, enhancing transport even if the original ECDI peak was seen to weaken from the FFT plots in figure 2.
5.3. Further examination of modulation regimes for effective transport suppression
In § 4.2, we identified in the 2-D radial-azimuthal configuration that a modulation frequency of
$f_{M}=40\,\text{MHz}$
with the amplitude of
$A_{M}=0.25E_{y}$
is the most effective to reduce electron transport among the studied modulation scenarios. Building on this finding, we now investigate how fine variations around this specific frequency and amplitude can further impact the instability spectra and the consequent electron transport. This analysis serves as a sensitivity study for fine-tuning frequency and amplitude, aiming to determine whether minor adjustments might lead to an even greater reduction in electron transport – an essential step in optimising these parameters for practical applications. To achieve this, the modulation frequency is varied from 10 to 100 MHz in 10 MHz increments, while maintaining a constant amplitude of
$0.25 E_{y}$
. Similarly, the modulation amplitude is varied incrementally from
$0.05 E_{y}$
to
$0.5 E_{y}$
, with the frequency fixed at 40 MHz.
The consequent variations in the frequency and wavenumber spectra with the modulation frequency and amplitude are presented in figures 12 and 13. Across all parameter studies, the peak associated with the ECDI is consistently weakened. Notably though, a modulation frequency around 40 MHz combined with an amplitude near
$0.25E_{y}$
still represents the most substantial impact on reducing the ECDI amplitude. This reduction in ECDI amplitude directly contributes to a significant decrease in electron transport, highlighting the effectiveness of this specific frequency-amplitude configuration in controlling this instability mode and enhancing transport suppression.

Figure 12. Variation of the frequency spectrum relative to the baseline case from the radial-azimuthal simulations with various modulation frequencies (
$10{-}100\,\text{MHz}$
) and amplitudes (
$0.05{-}0.5\, E_{y}$
). In the left-column plots, the modulation amplitude is fixed at
$0.25E_{y}$
and in the right-column plots, the modulation frequency is fixed at 40 MHz. In each case, the spectrum represents the average of the temporal FFT of azimuthal electric field (
$E_{z}$
) signal over all azimuthal and radial positions and over three simulation repetitions. The characteristic frequencies displayed in grey dashed lines from left to right correspond to theoretical ion acoustic frequency (
$f_{IA}$
), ion plasma frequency (
$f_{pi}$
), Hall circulation frequency (
$f_{E}=E_{y,0}/B_{x}L_{z}$
), and first and second harmonics of electron cyclotron frequency (
$f_{ce},2f_{ce}$
).

Figure 13. Variation of the azimuthal wavenumber (
$k_{z}$
) spectrum relative to the baseline case from the radial-azimuthal simulations with various modulation frequencies (
$10{-}100\, \text{MHz}$
) and amplitudes (
$0.05{-}0.5 E_{y}$
). In the left-column plots, the modulation amplitude is fixed at
$0.25E_{y}$
and in the right-column plots, the modulation frequency is fixed at 40 MHz. In each case, the spectrum represents the spatiotemporally averaged [over
$5{-}15\, \mu\text{s}$
(35–45
$\mu\text{s}$
absolute time) and all radial positions] spatial FFT of azimuthal electric field (
$E_{z}$
) signal averaged over three simulation repetitions. The yellow and blue lines indicate the theoretical
$k_{z}$
of the ECDI and the MTSI based on (2.1) and (2.2), respectively.
The average electrons’ axial current and mobility plots in figures 14 and 15 further suggest that modulation at 40 MHz with amplitude
$0.1E_{y}{-}0.3E_{y}$
represents the optimal range of configuration for the minimisation of electron cross-field transport in the studied 2-D plasma configuration.

Figure 14. Electrons’ axial current density (
$J$
) and electrons’ axial mobility (
$\mu$
) from the radial-azimuthal simulations with various modulation frequencies (
$10{-}100\, \text{MHz}$
). The modulation amplitude is fixed at
$0.25E_{y}$
.The colour bars represent the spatiotemporal mean value [over
$5{-}15\, \mu\text{s}$
(35–45
$\mu\text{s}$
absolute time) and entire domain] averaged over three simulation repetitions. Panels (a) and (c) display the absolute values of
$J$
and
$\mu$
, and panels (b) and (d) show the normalised change of
$J$
and
$\mu$
in the modulated simulations with respect to the respective quantities in the baseline case. The error bars indicate one standard deviation among the three simulation repetitions in each case. The solid and the dotted black lines represent the average and one standard deviation corresponding to three simulation repetitions for the baseline case.

Figure 15. Electrons’ axial current density (
$J$
) and electrons’ axial mobility (
$\mu$
) from the radial-azimuthal simulations with various modulation amplitude (
$0.05{-}0.5 E_{y}$
). The modulation frequency is fixed at 40 MHz. The colour bars represent the spatiotemporal mean value [over
$5{-}15\, \mu\text{s}$
(35–45
$\mu\text{s}$
absolute time) and entire domain] averaged over three simulation repetitions. Panels (a) and (c) display the absolute values of
$J$
and
$\mu$
, and panels (b) and (d) show the normalised change of
$J$
and
$\mu$
in the modulated simulations with respect to the respective quantities in the baseline case. The error bars indicate one standard deviation among the three simulation repetitions in each case. The solid and the dotted black lines represent the average and one standard deviation corresponding to three simulation repetitions for the baseline case.
6. Conclusion
This study demonstrated that externally applied electrostatic modulation can substantially influence instability-induced electron transport in E × B plasma discharges. Through one-and two-dimensional kinetic simulations, we explored how the frequency and amplitude of an imposed axial electric field affect the spectral energy distribution of plasma instabilities and the resulting cross-field transport.
In the 1-D azimuthal simulations, modulation at frequencies below 150 MHz, particularly near 40 MHz, consistently suppressed the ECDI, reducing its spectral amplitude and shifting energy towards longer wavelengths and lower frequencies. This spectral redistribution correlated with a notable decrease in axial electron transport, with the strongest reduction – up to 30 % – observed at 40 MHz. In contrast, modulation near or above the electron cyclotron frequency intensified high-frequency oscillations and led to a recovery or increase in transport.
The 2-D radial-azimuthal simulations exhibited similar trends at corresponding modulation frequencies. Nonetheless, modulation at 560 MHz produced more pronounced spectral broadening and stronger transport enhancement at high amplitudes, reaching increases of more than two orders of magnitude in axial current. At moderate amplitudes, the same modulation suppressed transport, highlighting the non-monotonic and amplitude-sensitive nature of the response.
The spectral analyses carried out via FFT and bicoherence revealed that electrostatic modulation reshapes the energy pathways among interacting wave modes. Sideband structures around the modulation frequencies, arising from nonlinear interactions between the modulation wave and a low-frequency mode excited near 2 MHz, were observed in the FFT spectra from the 1-D simulations for modulation frequencies
$\geq$
20 MHz, and in the 2-D cases, particularly at 40 MHz. Grid-like bicoherence patterns confirmed that these interactions also generate higher-order harmonics of the modulation-frequency sidebands. At lower modulation frequencies (1–6 MHz), broadband bicoherence features for the 1-D azimuthal data demonstrated extended phase coupling across a wide range of frequencies, consistent with coherent three-wave energy exchange and the corresponding spectral broadening seen in the FFT spectra.
Analysis of the variations in the normalised cross-correlation term
$\langle \tilde{n}_{e}\tilde{E}_{z}\rangle$
between electron-density and azimuthal-field fluctuations with modulation frequency – performed for the 1-D azimuthal cases as a demonstrative example – showed that, in the regimes where transport was suppressed (low-to mid-frequency modulation), the phase alignment between
$\tilde{n}_{e}$
and
$\tilde{E}_{z}$
weakened, reducing axial mobility. Conversely, at higher modulation frequencies, particularly near the electron-cyclotron frequency, stronger phase correlation enhanced transport even when the amplitude of the dominant instability peak was seen to diminish.
Together, these findings establish that external modulation can steer instability-induced transport in E × B plasmas by altering nonlinear mode coupling and the phase relationship between key fluctuating quantities. Within the mechanism-assessment scope of this work performed using physically idealised E × B plasma configurations, the results highlight the possibility of tailoring modulation parameters to regulate cross-field mobility and provide a foundation for future experimental validation, with the aim of controlling E × B plasma devices.
Acknowledgements
Editor Edward Thomas, Jr thanks the referees for their advice in evaluating this article.
Funding
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The simulation data that support the findings of this study are available from the corresponding author upon reasonable request.
Appendix: Supplementary results
A. Variations in plasma profiles and particles’ distribution functions in the radial-azimuthal configuration
The time-averaged radial profiles of plasma properties in the radial-azimuthal plasma configuration under different modulation conditions are presented in figure 16. Furthermore, the time-averaged particles’ velocity distribution functions along the radial and azimuthal directions from the corresponding simulations are provided in figure 17. These figures illustrate the influences of modulation with various frequencies and amplitudes, and the altered instability spectra, on the plasma properties and particles velocity spreads.

Figure 16. Time-averaged [over
$5{-}15\,\mu\text{s}$
(35–45
$\mu\text{s}$
absolute time)] radial profiles of the plasma properties from the radial-azimuthal simulations with various modulation frequencies and amplitudes averaged over all azimuthal locations and over three simulation repetitions. The rows from top to bottom represent ion number density (
$n_{i}$
), radial electron temperature (
$T_{ex}$
), azimuthal electron temperature (
$T_{ez}$
) and the ratio of the radial to azimuthal electron temperature (
$T_{ex}/T_{ez}$
). Note: in plots of
$n_{i}$
and
$T_{ez}$
for
$f_{M}=560\, \text{MHz}$
, the scale of the
$y$
-axis is different than the corresponding plots for the other modulation frequencies.

Figure 17. Time-averaged [over
$5{-}15\, \mu\text{s}$
(35–45
$\mu\text{s}$
absolute time)] electrons’ and ions’ velocity distribution functions (EVDF and IVDF, respectively) along the radial and azimuthal directions from the radial-azimuthal simulations with various modulation frequencies and amplitudes. The distribution functions correspond to the electrons or ions within the entire domain and are averaged over three simulation repetitions. Note: in the plots for
$f_{M}=560\, \text{MHz}$
, the scale of the
$x$
-axis is different than the corresponding plots for the other modulation frequencies.
Notably, for modulation at 560 MHz with
$A_{M}\geq 0.1E_{y}$
, remarkable heating of both electrons and ions in the radial and azimuthal directions is observed. This effect is evident in the temperature profiles (figure 16) and is manifested as the pronounced broadening of the velocity distribution functions (figure 17).
























































































