Manipulating plasma turbulence in cross-field plasma sources using unsteady electrostatic forcing

Also known as electron–cyclotron drift instability or beam cyclotron instability, the EDI is thought to be a primary factor in driving anomalous transport in HETs [15, 16]. Initially investigated in the context of collisionless longitudinal shocks in space plasmas in the 1960s and 1970s [6], the theory behind EDI formation is still an active topic of research. An EDI forms when electron Bernstein modes at high frequencies at multiples of the electron cyclotron frequency Doppler shift towards lower frequencies and the ion acoustic mode frequency range. Once within a similar frequency range, the two modes merge to form the EDI [27].

The dispersion relation, linking wavenumber k and frequency ω, for the EDI has been derived in both three dimensions and one dimension extensively by Ducrocq et al from the perturbed Vlasov equation [28]. Equation (2) shows the three-dimensional dispersion relationship

$$\begin{align} &1 + k^{2}\lambda_\mathrm{De}^{2} + g\left(\frac{\omega - k_{x}\text{v}_{\textbf{E}}}{\Omega_\mathrm{ce}}, (k_{y}^{2} + k_{z}^{2}),k_{y}^{2}r_\mathrm{ce}^{2}\right) \nonumber\\ &\quad- \frac{k^2\lambda_\mathrm{De}^2\omega_\mathrm{pi}^2}{(\omega - k_{z}\text{v}_\mathrm{b,i})^2} = 0, \end{align} \tag{ 2 }$$

where g is the Gordeev function, $\textbf{v}_\textbf{E}$ is the electron azimuthal drift velocity and $\textbf{v}_\mathrm{b,i}$ is the primary axial ion beam velocity [29]. The wave number vector in three dimensions is $\mathbf{k} = [k_x,k_y,k_z]^\mathrm{T}$. $\lambda_\mathrm{De}$, $\Omega_\mathrm{ce}$, $r_\mathrm{ce}$ and $\omega_\mathrm{pi}$ denote the Debye length, electron cyclotron frequency, the electron Larmor radius at thermal velocity $\text{v}_{\text{the}}$ and the ion plasma frequency, respectively.

Based on Ducrocq et al's analysis of equation (2) in the axial-azimuthal plane, it has been shown that the stability transitions for the EDI occur when $\omega - k_{x}\text{v}_\textbf{E} = n\Omega_\mathrm{ce}$. For the EDI ω is small in comparison to $\Omega_\mathrm{ce}$, resulting in resonance peaks when $k_{x}\text{v}_\textbf{E} \approx n\Omega_\mathrm{ce}$, $n = \pm\,1, 2, 3\ldots$. The resultant unstable wavenumbers, the cyclotron harmonics, are denoted nk₀ where $k_0 = \Omega_\mathrm{ce}/\text{v}_\mathbf{E}$ [30]. In the limit where the wavenumber parallel to magnetic field is sufficiently large, such as for an HET with an appropriately sized radial dimension, the relation reduces to an ion acoustic-like dispersion relation with approximate analytical solution [12, 31, 32]

$$\begin{align} \omega_R = \textbf{k}\cdot\textbf{v}_\mathrm{b,i}\pm\frac{kc_\mathrm{s}}{\sqrt{1+k^2\lambda_\mathrm{De}^2}} \end{align} \tag{ 3a }$$

$$\begin{align} \gamma \approx \pm\sqrt{\frac{{\pi}m_\mathrm{e} }{8m_\mathrm{i}}}\frac{\mathbf{k}\cdot\mathbf{v}_\textbf{E}}{\left(1+k^2\lambda_\mathrm{De}^2\right)^{3/2}} \end{align} \tag{ 3b }$$

where ω_R and γ are the angular wave frequency and growth rate of the wave, respectively, $c_\mathrm{s}$ is the ion sound speed and $m_\mathrm{i}$ and $m_\mathrm{e}$ are the ion and electron mass, respectively. Assuming the instability is primarily in the azimuthal dimension, $\textbf{k} \approx [k_{x}\, 0\, 0]^\mathrm{T}$ and therefore $k \approx k_{x}$. Setting $\partial\gamma/\partial k = 0$ finds k, ω_R and γ at maximum growth rate

$$\begin{align} [k]_\mathrm{max} = \frac{\sqrt{2}}{2\lambda_\mathrm{De}} \end{align} \tag{ 4a }$$

$$\begin{align} \left[\omega_R\right]_\mathrm{max} = \mathbf{k}\cdot\mathbf{v}_\mathrm{b,i}\pm\frac{\sqrt{3}\omega_\mathrm{pi}}{3} \end{align} \tag{ 4b }$$

$$\begin{align} \left[\gamma\right]_\mathrm{max} = \sqrt{\frac{\pi}{27}}\frac{\omega_\mathrm{pi}\text{v}_{\mathbf{E}}}{\text{v}_{\text{the}}} .\end{align} \tag{ 4c }$$

The wave phase and group velocities are defined by $\text{v}_\phi = \omega_R/k$ and $\text{v}_{g} = \partial\omega_R/\partial \mathbf{k}$, respectively. From equation (3a ) we find for the EDI

$$\begin{align} \text{v}_\phi(k) = \hat{\mathbf{k}}\cdot\textbf{v}_\mathrm{b,i}\pm\frac{c_\mathrm{s}}{\sqrt{1+k^2\lambda_\mathrm{De}^2}} \end{align} \tag{ 5 }$$

$$\begin{align} \!\!\!\quad\text{v}_g(k) = \mathbf{v}_\mathrm{b,i}\pm\hat{\mathbf{k}}\frac{c_\mathrm{s}}{\left(1+k^2\lambda_\mathrm{De}^2\right)^{3/2}} \end{align} \tag{ 6 }$$

where $\hat{\mathbf{k}} = \mathbf{k}/k$. Evaluating the EDI phase velocity at maximum growthrate we obtain $\text{v}_{\phi} = \sqrt{3/2}c_\mathrm{s}$.

The electron cross-field mobility (also known as cross-field transport) describes the electron movement in the axial dimension out of the highly magnetized region of the thruster and is parameterized by $\mu = u_{z}/E_{0}$. The cross-field electron mobility is roughly analogous to the axial electron current. Electron cross-field mobility is predicted from the electron momentum conservation equation [33]

$$\begin{align} m_{\alpha} \frac{\mathrm{D} \mathbf{u}_{\alpha}}{\mathrm{D} t}& =q_{\alpha}\left(\mathbf{E}+\mathbf{u}_{\alpha} \times \mathbf{B}\right)-\frac{1}{n_{\alpha}} \nabla p_{\alpha} \nonumber\\ &-m_{\alpha} \sum_\beta\nu_{\alpha \beta}\left(\mathbf{u}_{\alpha}-\mathbf{u}_{\beta}\right) \end{align} \tag{ 7 }$$

where $\mathrm{D}\textbf{u}_{\alpha}/\mathrm{D}t$ denotes the total derivative and $\nu_{\alpha\beta}$ the collision frequency of species α with species β. The flow velocity, electric and magnetic field strength, particle charge, density and pressure are denoted by u, E, B, q, n and p, respectively. Ignoring the pressure gradient and inertia terms in the axial and azimuthal directions, with a uniform magnetic field in the radial direction, B₀, the momentum conservation equation becomes

$$\begin{align} 0 = qn_\mathrm{e}E_{z}- qn_\mathrm{e}u_{x}B_{y} - m\nu_\mathrm{en}n_\mathrm{e}u_{z} \end{align} \tag{ 8a }$$

$$\begin{align} 0 = qn_\mathrm{e}E_{x}- qn_\mathrm{e}u_{z}B_{y} - m\nu_\mathrm{en}n_\mathrm{e}u_{x} .\end{align} \tag{ 8b }$$

Taking into account the presence of gradients and instabilities observed in the azimuthal direction of the thruster [34] and following the analysis of Lafleur [35], we calculate the effective cross-field transport parameter, taking average field quantities from the electron momentum conservation equations giving equation (9) where $\nu_\mathrm{en}$ is the electron-neutral collision frequency

$$\begin{align} \mu_{\text{effective}} = \frac{\frac{|q|}{m\nu_\mathrm{en}}}{1 + \frac{\Omega_\mathrm{ce}^{2}}{\nu_\mathrm{en}^2}}\left[1-\frac{\Omega_\mathrm{ce} < n_\mathrm{e}E_{x}>}{\nu_\mathrm{en}n_\mathrm{e}E_{z}}\right] .\end{align} \tag{ 9 }$$

For an even more detailed derivation of the cross-field mobility term, taking into account pressure gradient and inertia terms in the axial and azimuthal directions, see the work of Lafleur, Baalrud and Chabert [12, 31].

The underlying physics that govern the turbulent behaviour observed within HETs and the effects of unsteady forcing are outlined and modelled using a reduced-order one-dimensional three-velocity (1D-3 V) azimuthal fully kinetic electrostatic particle-in-cell (PIC) simulation code using Imperial College Plasma Propulsion Laboratory's (IPPL) PlasmaSim written in Julia [36]. Previous authors have concluded that an azimuthal 1D-3 V model is sufficient for capturing the plasma dynamics responsible for both the EDI development and the resultant anomalous electron mobility that is under investigation [37]. A simplified model of the HET reduces the complexity of the interactions between the many different plasma modes that exist in a real thruster, making analysing the links between the two properties under investigation far more tractable along with having a significantly reduced computational cost.

PIC simulation is a kinetic method for plasma simulation, representing electrons and ions as groups of 'super-particles' moving according to the particle equations of motion on a grid of cells used for calculating the statistical thermodynamic and electromagnetic properties within a plasma [38]. PlasmaSim has been validated for use in three 1D-3 V cases for the axial, radial and azimuthal directions of an HET, respectively, against other plasma models and experimental data [39]. The branch used in this investigation is derived from the azimuthal 1D-3 V electrostatic method popularized by Birdsall and Langdon [38] using the hard sphere model of Bird to model particle collisions [40]. Unless otherwise stated, it should be assumed that the PIC model is identical to that outlined by Birdsall and Langdon. The directions x, y and z denote the azimuthal, radial and axial directions of the HET, respectively, as shown in figure 2.

Zoom In Zoom Out Reset image size

Given the technical complexity of producing a functioning PIC code, it is necessary to rigorously validate the results of any new code with respect to the available test cases and benchmarks provided by the community [7]. Furthermore, in order to identify the effects of applying unsteady forcing to the plasma, a suitable baseline simulation configuration to compare to had to be established. LANDMARK Test Case 1 [37] was identified as a suitable baseline for both validation and comparison purposes.

Figure 2(a) shows a diagram of the physical locations of the simulated dimension (dashed line), along which Poisson's equation is satisfied, and the particle-tracking domain (highlighted surface) in PlasmaSim on a reference SPT-100 thruster, reproduced with permission from Tejeda et al [24]. Figure 2(b) shows the two-dimensionalized version of the simulated dimension and particle tracking domain 'unwrapped' from the annular channel along with particle boundary conditions.

To accurately reproduce the baseline test case, the following boundary conditions were implemented for the particles and electric field. Wherever possible, the boundary conditions were chosen to replicate those used by Lafleur et al [35] to allow for a comparison of the results to those in the literature. Initially, electrons, neutrals and ions are all generated with a randomized Maxwellian thermalized velocity at a given temperature, $T_{\mathrm{e}0}$, $T_{\mathrm{n}0}$ and $T_{\mathrm{i}0}$, respectively, and randomized position within the particle-tracking region. Sizes of the super-particles are set according to specified number densities, n₀ for electrons and ions and $n_{n0}$ for neutrals. The azimuthal boundary condition is periodic, with Poisson's equation solved numerically using the finite difference method with Dirchelet boundary conditions of $\phi_{x = 0} = \phi_{x = X_\mathrm{max}} = 0$ where particles passing outside the limits of $x\gt X_\mathrm{max}$ or x < 0 are re-introduced at the same velocity on the other side of the domain. In order to simplify the modelling of the radial dimension of the HET, another periodic boundary condition was used which again reintroduced particles that left the extents of the simulation domain, $y\gt Y_\mathrm{max}$ or y < 0, on the other side of the domain. In the axial dimension, electrons passing out of the domain at the anode end (z = 0) were eliminated before being reintroduced at a random azimuthal location at the cathode end ($z = Z_\mathrm{max}$) with a thermalized velocity of temperature $T_{\mathrm{e}0}$. The axial component of the velocity vector is enforced to be into the particle tracking domain to prevent hot electrons being eliminated straight after reinitialization. Similarly, for the ions, upon reaching the cathode end they are eliminated and reintroduced in a similar fashion at the anode end with temperature $T_{\mathrm{i}0}$. Neutral particles are specularly reflected at the axial particle tracking boundary. The static electric and magnetic field strengths, E₀ and B₀, are set in the axial and radial directions, respectively.

Elastic collisions between charged particles and neutrals were captured in the model using the direct-simulation Monte Carlo method with expected particle collision frequencies between species α and β obtained using the hard sphere model of Bird, where the collision frequency $\nu_{\alpha\beta} = n_\alpha \overline{\sigma_{\alpha\beta}\text{v}}$ and the collision cross-section $\sigma_{\alpha\beta} = \pi(r_{\alpha}+r_{\beta})^2$ where r is the van der Walls radius of the colliding particle [40]. Ionization reactions between electrons and neutral particles were not captured as part of the model.

Given the statistical nature of PIC simulation, Julia's [36] inbuilt random number generator using the Xoshiro256++ algorithm is used to generate independent trial runs for each scenario using a different random number generator seed before averaging to remove the effects of statistical noise. Limitations on the simulation parameters for the spatial discretization of the modelled dimension ($\Delta x$) are given by the requirement that $\Delta x$ must be less than $\lambda_\mathrm{De}$ to resolve the Debye length of the plasma and not simply capture collective effects. The time step of the simulation is required to resolve both the electron plasma ($\Delta t_\mathrm{pe}$) and cyclotron frequency ($\Delta t_\mathrm{ce}$) from the Nyquist criterion, and satisfy the Courant–Freidrichs–Lewey condition ($\Delta t_\mathrm{CFL}$) which are given by below [38]

$$\begin{align} \Delta t_\mathrm{pe} < 2\pi/\omega_\mathrm{pe} \end{align} \tag{ 10a }$$

$$\begin{align} \Delta t_\mathrm{ce} < 2\pi/\Omega_\mathrm{ce} \end{align} \tag{ 10b }$$

$$\begin{align} \Delta t_\mathrm{CFL} < \Delta x/|\text{v}_\mathrm{the}| \end{align} \tag{ 10c }$$

$$\begin{align} \Delta x < \lambda_\mathrm{de} .\end{align} \tag{ 10d }$$

An additional safety factor of 25 is then applied to the minimum time step.

The electron cross-field mobility parameter can be calculated numerically from the simulation particles by evaluating the average axial velocity of the electrons according to equation (11) where N is the number of electron super-particles at each time step.

$$\begin{align} \mu_{\text{PIC}} = \frac{\sum_{j=1}^{N}\text{v}_{jz}}{NE_{0}} .\end{align} \tag{ 11 }$$

The key publications laying out the conditions and expected results of the baseline configuration are given in 'Nonlinear structures and anomalous transport in partially magnetized E × B plasmas' by Janhunen et al [30] and 'Theory for the anomalous electron transport in Hall effect thrusters I & II' by Lafleur et al [32, 35], which were used in the validation of the results obtained by PlasmaSim. The baseline plasma properties used for the initialization of the simulation are given in table 1. Where possible, the values of the plasma properties were selected to be identical to those used by Lafleur et al [35] to validate the results generated by PlasmaSim against the results available in the literature. The neutral number density was arbitrarily selected so that the neutral gas pressure would be sufficiently low to observe the instability [35].

Table 1. Baseline operating and numerical parameters used in the PIC simulations.

Plasma property	Value
$T_\mathrm{e0}$ (eV)	2.0
$T_\mathrm{i0}$ (eV)	0.1
$T_\mathrm{n0}$ (K)	300
B0 (T)]	0.02
E0 (V cm−1)	200
$n_\mathrm{n0}$ (m−3)	$1.0\times 10^{-18}$
n0 (m−3)	$1.0\times 10^{-17}$

Simulation parameters such as the time step, data capture frequency and spatial discretization were tuned to satisfy the conditions given by equations (10a )–(10d ) to reduce numerical inaccuracy of the model and achieve an accurate depiction of the EDI and cross-field mobility development. The final simulation parameters settled on are shown in table 2. Again, where possible, the baseline simulation parameters were chosen to closely resemble those used by Lafleur et al [35]. Results from a mesh-refinement study confirmed that a sufficient number of cells (NG) for resolving the instability was 175 in the x dimension with a time step (t_step) of $1\times10^{-11}$ s. A data capture rate (t_write) of $1\times10^{-10}$ s was selected by Nyquist's criterion to capture interesting plasma behaviour at frequencies up to and including the electron plasma frequency whilst reducing the size of the simulation output files. The extents of the simulation domain X_max, Y_max and Z_max were chosen as 0.005 m, 0.0075 m and 0.01 m, respectively. One hundred super particles per cell for each species ($N/NG$) were selected from the literature as an appropriate number for initialization of the simulation.

Table 2. Simulation parameters.

Simulation parameter	Value
tstep (s)	$1.0\times10^{-11}$
twrite (s)	$1.0\times10^{-10}$
NG	175
$N/NG$	100
Xmax (m)	0.005
Ymax (m)	0.0075
Zmax (m)	0.01

The primary method to affect unsteady forcing was by varying the axial electric field, which was thought could provide a mechanism for manipulating azimuthal instabilities and hence the cross-field mobility of the plasma. It has previously been observed that forcing the axial electric field can lead to changes in the temporal spectral energy distribution of the azimuthal electric field; however, specific effects of forcing over a wide band of forcing frequencies have yet to be fully characterized.

In a real HET, axial electric field forcing can be accomplished by varying the potential of the anode or the cathode to achieve the desired electric field. This has previously been suggested by Tejeda et al [24] but not yet been tested experimentally. The axial electric field in an HET is non-uniform, as shown in figure 1, but, given that the 1D-3 V azimuthal version of PlasmaSim does not resolve Poisson's equation in this direction, the axial electric field profile is instead simplified to a uniform field varying in time with field strength given by $E_{z} = E_{0} + A_\text{EF}\sin{\omega_\text{EF} t}$, where E₀ is the average static electric field strength and $A_\text{EF}\sin{\omega_\text{EF} t}$ is the superimposed oscillation. A_EF is the amplitude and $\omega_\text{EF}$ is the frequency of the forcing of the electric field, t is the simulation time.

The plasma distributions obtained in PlasmaSim for the parameters specified for the baseline test case are shown in figure 3. Both the electron density and the azimuthal electric field show wave-like behaviour with high wavenumber oscillations over the first 500 ns giving way to low-wavenumber, progressively more chaotic, high-amplitude oscillations as the instability saturates towards 4000 ns.

Zoom In Zoom Out Reset image size

The wavenumber and frequency of the wave-like behaviour can be calculated from the electric field distribution by counting peaks along lines XX and YY of figure 4 respectively after saturation. The wave phase velocity can be evaluated from $\Delta x/\Delta t$. This method gives a wave phase velocity of $5.4\times10^{3}$ m s⁻¹ , wave number of $4.1\times10^{3}$ rad m⁻¹ and a frequency of $2.1\times10^{7}$ rad s⁻¹. Given the electron temperature at 4500 ns is approximately $21\times T_{\mathrm{e}0}$ the Debye length of the plasma would have increased to $1.5\times10^{-4}$ m. Using equations (4a ), (4b ) and (6), we calculate the expected wavenumber at maximum growth rate $[k]_\mathrm{max}$ at 4500 ns of $4.6\times10^{3}$ rad m⁻¹, frequency $[\omega_R]_\mathrm{max}$ of $2.1\times10^{7}$ rad s⁻¹ and phase velocity $\text{v}_\phi([k]_\mathrm{max})$ of $5.5\times10^{3}$ m s⁻¹ for the EDI. The results show good agreement with those obtained from the analytical dispersion relationship for the wavenumber of maximum growth rate. Phase velocity, wavenumber and frequency show a 2.8%, 11.7% and 0.5% error from the theoretical values, respectively. This was deemed to be within acceptable bounds given the assumptions of the theoretical model of an unbounded collisionless plasma and suggests PlasmaSim is accurately modelling the expected physics of the EDI. k, ω and $\text{v}_\phi$ of the resultant wave behaviour appears to be closely related to the wavenumber of maximum growth rate $[k]_\mathrm{max}$.

Zoom In Zoom Out Reset image size

It should be noted that in the numerical model both the electron temperature and electron number density oscillate between approximately $30\times T_\mathrm{e0}$ and $10\times T_\mathrm{e0}$, and $1.8\times10^{17}$ m⁻³ and $0.6\times10^{17}$ m⁻³, respectively, after the instability saturates. The values of the phase velocity, wavenumber and frequency calculated from equations (4a ), (4b ) and (6) do not account for these oscillations, using instead the instantaneous values for electron temperature and number density at 4500 ns of $21\times T_\mathrm{e0}$ and $1.0\times10^{17}$ m⁻³, respectively. This is a likely explanation for the remaining discrepancy between the numerically calculated wavenumber and the theoretical value as the theoretical value is sensitive to these parameters.

The dispersion relationship of the EDI was calculated numerically by taking a two-dimensional fast Fourier transform of the azimuthal electric field from PlasmaSim in MATLAB [41] to convert it to a wavenumber–frequency diagram, as shown in figure 5. The numerical and analytical dispersion relationships show good agreement, with the wavenumber of the largest growth rate in figure 5 being close to the analytically calculated value of $[k]_\mathrm{max}$ and $[\omega_R]_\mathrm{max}$ for the EDI from equations (4a ) and (4b ). These results were compared to Lafleur et al's [12] results finding the same dispersion relationship despite Lafleur et al using a 2D-3 V axial–azimuthal PIC model. This gives confidence in the one-dimensional model retaining enough complexity to model the EDI sufficiently despite the dimensionality reduction.

Zoom In Zoom Out Reset image size

The electron cross-field mobility calculated from equation (11) and averaged over the time period after instability saturation (6 µs–10 µs) was found as 5.38 m² V⁻¹ s⁻¹. To reduce the effects of statistical variation, the result was averaged over five independent simulation runs with a standard deviation of 0.94 m² V⁻¹ s⁻¹, agreeing with the cross-field mobility parameter obtained by Lafleur [35] for similar simulation parameters.

Initial investigations focused on applying forcing at frequencies at or near the fundamental plasma frequencies, the electron cyclotron frequency $\Omega_\mathrm{ce}$, the electron plasma frequency $\omega_\mathrm{pe}$ and the ion plasma frequency $\omega_\mathrm{pi}$. Three forcing amplitudes A_EF of 10 V cm⁻¹, 50 V cm⁻¹ and 100 V cm⁻¹ were tested. The resulting cross-field mobility increased significantly from the baseline value by up to a factor of 20 times when forcing at the electron cyclotron frequency with the largest increase at a forcing amplitude of 100 V cm⁻¹ and the smallest increase at an amplitude of 10 V cm⁻¹. Only a small increase in cross-field mobility was noted when forcing at the electron plasma frequency, as shown in figure 6; however, this was not judged to be statistically significant. A slight reduction in cross-field mobility was observed forcing at the ion plasma frequency; however, this was only statistically significant at a forcing amplitude of 50 V cm⁻¹.

Zoom In Zoom Out Reset image size

The temporal spectral amplitude distribution for the azimuthal electric field was compared with the baseline distribution for each forcing frequency. To reduce the effects of statistical noise, the amplitude spectrum for the baseline test case was averaged over five independent simulation runs to obtain the expected distribution whilst the amplitude spectrum for the forced case was averaged over three independent runs. Significant divergences from the baseline distribution were identified for forcing frequencies matching the electron cyclotron frequency and the ion plasma frequency, as can be seen from figures 7 and 8, respectively. Forcing at the electron cyclotron frequency resulted in a reduction in spectral amplitude around the frequency of maximum growth rate of the EDI $[\omega_R]_{\mathrm{max}}$ whilst all other frequencies saw an increase in amplitude, particularly for frequencies exceeding the 'so-called' Hall circulation frequency $\omega_{\textrm {E}} = 2\pi\text{v}_{\textrm {E}}/X_\text{max}$ where $\text{v}_{\textrm {E}}$ is the drift velocity $E_0/B_0$. Conversely, forcing at the ion plasma frequency leads to an increase in spectral amplitude for frequencies between $\omega_\mathrm{pi}$ and ω_E with a reduction in the spectral amplitude for frequencies greater than ω_E or lower than $\omega_\mathrm{pi}$, the amplitude of the peaks associated with the electron cyclotron frequency and its resonances appeared to be damped when forcing at this frequency was applied.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To attempt to fully characterize the effects of unsteady forcing, a wider range of frequencies from the MHz range up to the GHz range were investigated. Simulations were run with forcing frequencies $\omega_\mathrm{EF}$ from $1\times10^{6}$ Hz to $1\times10^{10}$ Hz in logarithmic steps, i.e 10⁶ Hz, 10^6.1 Hz, 10^6.2 Hz etc, to encompass the whole range. Again, three forcing amplitudes $A_\mathrm{EF}$ of 10 V cm⁻¹, 50 V cm⁻¹ and 100 V cm⁻¹ were tested to quantify the effects of varying forcing amplitude. The resultant power spectrum data and associated cross-field mobility were averaged over three independent simulation runs similarly to above.

Approximate forcing effects can be grouped roughly according to the qualitatively observed effects in different frequency bands. To aid discussion, the 'low' frequency band shall be frequencies below the ion plasma frequency $\omega_\mathrm{pi}$, 'mid-range' frequencies shall be between $\omega_\mathrm{pi}$ and the Hall circulation frequency $\omega_\text{E}$ and 'high' frequencies shall be frequencies above the Hall circulation frequency. Forcing frequencies $\omega_\mathrm{EF}$ in the low frequency range, below the ion plasma frequency, typically resulted in spectral amplitude reduction in the low frequency range, as demonstrated in figure 9. Furthermore, above $\omega_\mathrm{pi}$, the spectral amplitude showed a characteristic reduction, particularly around the electron cyclotron frequency resonances and the Hall circulation frequency. No obvious resonant behaviour was observed at multiples of the forcing frequency in this range.

Zoom In Zoom Out Reset image size

Mid-range forcing for $\omega_\mathrm{EF}$ between the ion plasma frequency $\omega_\mathrm{pi}$ and the Hall circulation frequency $\omega_\text{E}$ resulted in interesting behaviour whereby large peaks in amplitude centred on the forcing frequency $\omega_\text{EF}$ and its resonances appear in the spectral amplitude distribution for the azimuthal electric field. An example of this behaviour is shown in figure 10, demonstrating a coupling between the axial and azimuthal electric field behaviour. Mid-range forcing also appears to reduce the spectral amplitude in both the high- and low- frequency ranges. As the forcing frequency approaches $\omega_\text{E}$, the amplitude of the electron cyclotron peak and its resonances become increasingly damped, potentially indicating a coupling between the mid-range frequency modes and the high frequency electron cyclotron modes.

Zoom In Zoom Out Reset image size

Figure 11 shows a characteristic response of the plasma when forced within the high-frequency range, demonstrating the observed phenomena of 'frequency lock-in'. When the plasma is forced at a frequency close to the electron cyclotron frequency, the peak typically located at $\Omega_\mathrm{ce}$ is observed to shift slightly away from its usual position towards the forcing frequency or lock-in frequency. Similar behaviour is observed at the resonances of the electron cyclotron frequency, sometimes with a significant amplitude reduction or peak sharpening at the other resonances. When $\omega_\mathrm{EF}$ is not located near an electron cyclotron frequency resonance, no significant change in the spectral amplitude distribution is observed. Forcing close to the Hall circulation frequency $\omega_\text{E}$ and its resonances appears to amplify the other Hall circulation resonances.

Zoom In Zoom Out Reset image size

Generally speaking, increasing the amplitude of the forcing function appears to increase the change in the observed spectral amplitude. Forcing at 10 V cm⁻¹ does not appear to significantly reduce the low-frequency spectral amplitude compared to forcing at 50 V cm⁻¹ when $\omega_\text{EF}$ is in the low-frequency range. When $\omega_\text{EF}$ is in the mid- to high-frequency region, the response for 10 V cm⁻¹ amplitude is similar to that observed for the higher forcing amplitudes tested.

The cross-field mobility calculated for the plasma at each forcing condition over a range of frequencies is shown in figures 12–14 for forcing amplitudes of 10 V cm⁻¹, 50 V cm⁻¹ and 100 V cm⁻¹, respectively, and compared to the expected value calculated from the unforced baseline case. Taking into account the statistical variation in the result output by the PIC algorithm, the figures show both the mean value and the region of one standard deviation from the mean. The mean baseline cross-field mobility is represented by a horizontal red line with the standard deviation of this shown by two dashed horizontal red lines. Statistically significant results that showed a deviation of the cross-field mobility from the baseline value were identified using a two-sample, two-tailed t-test at the 5% confidence level and are circled in red.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As expected from the prior investigation of plasma fundamental frequencies, for all three forcing amplitudes tested the cross-field mobility peaks when the forcing frequency is close to $\Omega_\mathrm{ce}$ with cross-field mobility increasing with amplitude. There is a reduction in cross-field mobility by a factor of up to a half relative to the baseline test case for forcing frequencies close to the frequency of the fundamental cyclotron harmonic $\omega_R(k_0)$ and its resonances; this is most pronounced for a forcing amplitude of 50 V cm⁻¹. At the same forcing amplitude, the cross-field mobility appears to be reduced somewhat for forcing frequencies in the range between the ion plasma frequency and the Hall circulation frequency, and no such effect is distinguishable for the other trial forcing amplitudes. Forcing at 10 V cm⁻¹ appears to result in the smallest deviation in cross-field mobility to the baseline case for low- and mid-range frequencies.

The reduction in cross-field mobility in the low frequency range is centred upon the frequency of the EDI's unstable cyclotron harmonics $\omega_R(k_0),\,\omega_R(2k_0),\,\omega_R(3k_0),\,\ldots,\,\omega_R(nk_0).$ The frequencies of these harmonics are calculated from the dispersion relationship for the EDI from the equation of angular wave frequency (3a ). Substituting in the unstable cyclotron wavenumbers nk₀ to this equation gives the wave frequency as

$$\begin{align} \omega_R(nk_0) = \pm\frac{nk_0c_\mathrm{s}}{\sqrt{1+(nk_0)^2\lambda_\mathrm{De}}} \end{align} \tag{ 12 }$$

where $n = 1,2,3\ldots$ and $\lambda_\mathrm{De}$ is a function of the electron temperature. Forcing at the frequency $\omega_R(nk_0)$ was observed to reduce the temporal spectral amplitude of frequencies other than those corresponding to an unstable wavemode in the low frequency region, potentially as a result of energy pathways amplifying the unstable EDI mode or as a result of frequency lock-in type behaviours. Figure 15 demonstrates this behaviour for forcing at the frequency of the fundamental cyclotron harmonic $\omega_R(k_0)$ for the temporal spectral distribution. Observing the region of the spectrum between $1\times10^{7}$ rad s⁻¹ and $1\times10^{8}$ rad s⁻¹, the spectral amplitude for the forced plasma shows peaks at $\omega_R(k_0)$, $\omega_R(2k_0)$ and $\omega_R(3k_0)$ with a reduction in spectral amplitude everywhere else compared to the baseline case. Similar results were observed for cases where forcing is applied at the frequencies of two other unstable cyclotron harmonics $2k_0$ and $3k_0$.

Zoom In Zoom Out Reset image size

Comparing the development of the spatial spectral amplitude distribution for the azimuthal potential field of the baseline case shown in figure 16 with that forced at the frequency of the fundamental cyclotron harmonic $\omega_R(k_0)$ shown in figure 17, it can be seen that forcing induces a reduction in spectral amplitude for wavemodes that are not excited, $2k_0$, $3k_0,\ldots$, whilst the excited wavemode k₀ sees an increase in amplitude. Counterintuitively, we see a reduction not an increase in cross-field mobility by forcing at the frequency of the unstable modes of the EDI. Referring back to equation (9), it is known that the cross-field mobility is dependent on the cross-correlation term $\lt n_eE_x\gt$. Calculating this term for the baseline case and comparing it with the calculated term for the forced case at the frequency of the fundamental cyclotron harmonic at 50 V cm⁻¹ shows a reduction of over 50% for the forced case. One possible theory for the reduction could be due to the forcing causing E_x and $n_\mathrm{e}$ to oscillate out of phase, reducing $\lt n_{e}E_x\gt$ and the resultant cross-field mobility.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It was observed from the spatial Fourier transform of the potential field that $|\phi_x(k)|^2\propto k^{-3}$, as illustrated in figure 16. This establishes an energy–wavenumber relationship seemingly similar to those predicted by Kolmogorov cascade theory whereby a turbulent field's spectral amplitude is proportional to k⁻ⁿ where n is a constant to be determined [42, 43]. The $|\phi_x(k)|^2\propto k^{-3}$ relationship has been observed to hold for instabilities in tokamak fusion plasmas studied by Qi [44], although more work is required to determine if there is any applicable similarity to tokamak instabilities here that could also be exploited for HETs and on the significance of the $n = -3$ relationship overall.

Forcing close to the electron cyclotron frequency led to a large increase in cross-field mobility. It is theorized that this could be a result of disruption to the circular motion of the electrons normally oscillating at that frequency in the E × B field, causing the electrons to break confinement. Whilst promoting the electron mobility classically would seem to be a detriment to HET systems, due to the reduction in current utilization efficiency, it could also be a pathway to reducing the anode voltage of such a system, enabling the possibility of a direct drive architecture.