# Fluid and hybrid simulations of the ionization instabilities in Hall
thruster

O. Chapurin<EMAIL_ADDRESS>Department of Physics and Engineering
Physics, University of Saskatchewan, Saskatoon SK S7N 5E2, Canada  A. I.
Smolyakov  Department of Physics and Engineering Physics, University of
Saskatchewan, Saskatoon SK S7N 5E2, Canada  G. Hagelaar  J.B. Boeuf  LAPLACE,
Université de Toulouse, CNRS, INPT, UPS, 118 Route de Narbonne, 31062
Toulouse, France  Y. Raitses  Princeton Plasma Physics Laboratory, Princeton,
New Jersey 08540, USA

###### Abstract

Low-frequency axial oscillations (5-50 kHz) stand out as a pervasive feature
observed in many types of Hall thrusters. While it is widely recognized that
the ionization effects play the central role in this mode, as manifested via
the large scale oscillations of neutral and plasma density, the exact
mechanism(s) of the instabilities remain unclear. To gain further insights
into the physics of the breathing mode and evaluate the role of kinetic
effects, a one-dimensional time-dependent full nonlinear low-frequency model
describing neutral atoms, ions, and electrons, is developed in full fluid
formulation and compared to the hybrid model in which the ions and neutrals
are kinetic. Both models are quasineutral and share the same electron fluid
equations that include the electron diffusion, mobility across the magnetic
field, and the electron energy evolution. The ionization models are also
similar in both approaches. The predictions of fluid and hybrid simulations
are compared for different test cases. Two main regimes are identified in both
models: one with pure low-frequency behaviour and the other one, where the
low-frequency oscillations coexist with higher frequency oscillations (with
the characteristic time scale of the ion channel flyby time, 100-200 kHz). The
other test case demonstrate the effect of a finite temperature of injected
atoms which is shown to have a substantial effect on the oscillation
amplitude.

##  I Introduction

Hall thrusters are successfully used for electric propulsion in space, e.g.
for satellite orbit keeping, and becoming an enabling technology of choice for
long-term missions, such as trips to Mars. Despite the relatively long history
of practical use (since 1972morozov2000fundamentals ), the crucial physical
aspects of their operation are poorly understood. In the absence of predictive
modeling capabilities, scaling of these devices for large (e.g., for long-term
missions) and for low (for microsatellites) power is very difficult and
expensive. The quantitative understanding of the physics of these devices
remains an important task.

Plasma discharge in Hall thrusters (HT) is supported by the electrons drifting
in the closed (periodic) azimuthal $\mathrm{E}\times\mathrm{B}$ direction,
while the thrust is created by ions (effectively unmagnetized due to large
gyroradius) accelerated by the electric field in the axial direction. One of
the characteristics of Hall thrusters is the presence of turbulence and
structures (azimuthal and axial) that affect their operation. Studies of
nonlinear phenomena in these plasmas are not only of great practical
importance but also address fundamental problems of plasma physics and plasma
turbulence. In particular, the turbulent electron transport in such devices is
orders of magnitude larger than the classical collisional transport (across
the magnetic field) predicts. Inhomogeneous plasmas with
$\mathrm{E}\times\mathrm{B}$ electron drift typically prone to various drift
instabilities, both due to the fluid smolyakov2016fluid  and kinetic
mechanisms lafleur2016theory ; janhunen2018nonlinear ; janhunen2018evolution ;
villafana2021rz ; jimenez2021twod , which drive high cross-field electron
currents (for more details see Ref. kaganovich2020physics, ).

Among the plethora of wave phenomena in a Hall thruster device, low-frequency
oscillations propagating in axial direction stand as one of the most common
and observed in most types of Hall thrusters morozov2000fundamentals . They
appear as the axial discharge current oscillations with frequencies of 5-50
kHzesipchuk1974plasma . A strong periodic depletion of atoms in the ionization
region is observed during the oscillations, suggesting the ionization nature.
In the literature they are known as breathing modes (due to slow periodic
plasma bursts out of the channel exhaust). Analytical studies of these
phenomena are difficult due to the importance of nonlinear effects and the
global nature of solutions, thus numerical methods have to be used.
Qualitatively, the oscillation period of breathing modes depends on the travel
time of neutral particles to the ionization region, e.g. for the
characteristic $1\text{\,}\mathrm{c}\mathrm{m}$ and atom velocity
$150\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$ gives
$15\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$. Overall, we understand the
phenomenology of the oscillations but cannot accurately predict their
existence and amplitude. Generally accepted phenomenological description of
these oscillations is described as the following sequence: decrease of the
discharge current $\,\to\,$ decrease of ionization $\,\to\,$ increase of the
neutral density in the exhaust region $\,\to\,$ increase of the electron
conductivity in that region $\,\to\,$ increase of the current and ionization
$\,\to\,$ neutral depletion $\,\to\,$ decrease of the current and so on.
According to this picture, the oscillation frequency is related to the time
necessary for the neutrals to refill the ionization region.

The 0D predator-prey model proposed earlier fife1997numerical ;
barral2008origin ; hara2014perturbation  is appealing because of its
simplicity but fails to identify the conditions for the instability. Moreover,
more accurate treatments show that the basic two-component (plasma-neutral)
system with uniform ion and neutral velocities is stable hara2014perturbation
; chapurin2021mechanism ; LafleurJAP2021 . A simple model was proposed that
the ion back-flow region which occurs near the anode as a result of large
contribution of the electron diffusion current (due to the density gradient)
and quasineutrality constrain provides a critical excitation mechanism for the
breathing mode chapurin2021mechanism . Linear resistively unstable modes
chable2005numerical  and fluctuations of electron temperature and power
absorption hara2014perturbation ; LafleurJAP2021  were also investigated as
possible triggers of the breathing modes.

In general, several physical mechanisms affect the breathing mode excitation
and characteristics: electron momentum and energy losses to the wall,
anomalous cross-field transport and heating, the ion backflow, and
recombination at the anode. These mechanisms are inter-related, depend in a
complex way on the magnetic field configuration, and are not easily
quantifiable. Numerical models that include many of these effects were
proposed makowski2001review ; hagelaar2004modelling ; boeuf1998low ;
barral2001fluid ; barral2008origin ; barral2009low . However, some calibration
and adjustment of the parameters are required to satisfactorily reproduce the
breathing modes characteristics observed experimentally giannetti2021numerical
. Therefore further insights on key physical processes are required to expand
the predicting powers of such models, especially to new parameters range and
new operational regimes.

While many numerical models for breathing modes based on fully fluid
formulations, time-dependent, the hybrid modeling was also undertaken using
the kinetic description for ions and neutrals boeuf1998low ; MorozovPPR2000h ;
chable2005numerical ; barral2001fluid ; ShashkovPoP2017 ; gavrikov2021hybrid .
The extent to which the ion and neutrals kinetic effects influence the
breathing mode excitation and characteristics remain a mute point of many
studies. The goal of this paper is to analyze the role of ions and neutrals
kinetic effects under the same physics of the electron dynamics which is
treated with the fluid theory. We use the axial one-dimensional full fluid and
hybrid models and compare their results.

The basic fluid model describes ion and atoms with first two fluid moments
(conservation of mass and momentum), and electrons are considered in drift-
diffusion approximation with a full electron energy balance. In the hybrid
model heavy species (ions, atoms) are kinetic (via particle-in-cell method)
and electrons are fluid (modelled in exactly the same way in both approaches).
Both models include plasma recombination at the anode and neutral dynamics
with ionization due to electron-neutral impact. Plasma discharge is supported
by the ionization process driven by the axial current due to the applied
potential across the domain. For the fluid simulations BOUT++ computational
frameworkdudson2009bout++  is used. The hybrid code was developed in the
LAPLACE laboratory, France hagelaar2002two ; hagelaar2003role ;
hagelaar2004modelling . The parameters of the simulations are chosen according
to the “Fluid/Hybrid” test case in the LANDMARK (Low temperAture magNetizeD
plasMA benchmaRKs) benchmarking projectlandmark . Previously, a comparison
between fluid and hybrid models for the axial direction of Hall thruster
configuration was presented in Ref. barral2001fluid, , however, this model did
not include electron pressure gradients, thus omitting the effects of electron
diffusion resulting in the formation of the presheath region near the anode
and ions transition through the ion-sound barrier. Neither it included the
full electron energy balance.

One of the important finding of the present paper is the identification of two
distinct regimes of breathing oscillations; the result which was confirmed
with both fluid and hybrid models during this benchmark. We show that the
regime with higher electron energy losses exhibits the low-frequency mode
(${\sim}$14\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$) that coexists with the
higher frequency ion “transient-time” oscillations
(${\sim}$150\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$$)esipchuk1974plasma . In
the second regime, with low electron energy losses, purely breathing
oscillations are observed, the so-called solo regime. We believe that
different mechanisms are involved in these regimes.

For the first regime, we identify the higher frequency oscillations as the
excitation of the resistive modes (convective instability with the
characteristic ion fly-by frequency)litvak2001resistive ; chable2005numerical
; fernandez2008growth ; koshkarov2018current . Such resistive-type modes
appear in simple models without ionization or electron diffusion. The main
feature of resistive modes is a strong dependence of growth rate and frequency
on the electron mobility (resistivity)fernandez2008growth . Similar features
are shown in this work, while that the low frequency mode (breathing mode) has
weak or no dependency on the electron mobility (see Appendix A).

In fact, the frequency of resistive modes can vary significantly, ranging
0.1-10 MHz, and may become close to the breathing modes at the lower end of
its spectrum. Some axial thruster models (with ionization but without electron
diffusion)morozov2000fundamentals ; MorozovPPR2000f  claimed that these modes
might be responsible for the breathing modes observed in Hall thrusters. For
clarity, here we will call breathing modes only those associated with atom
depletion, whose frequencies scale accordingly to atom fly-by time and
ionization processes.

In our earlier work, we have proposed a reduced model (only ion and atom
dynamics included) for the second, solo regime of the breathing mode. In this
regime, the instability is triggered by the ion backflow (negative ion
velocity) in the near-anode presheath region chapurin2021mechanism . It was
demonstrated that such configuration is prone to low-frequency oscillations,
where the ion backflow region is necessary. Recently, it has been shown in
Ref. gavrikov2021some,  with a more rigorous formal analysis that the sign-
alternating ion velocity profile with a positive slope (i.e. negative ion
velocity near anode and positive towards the exit) indeed is a necessary
condition to excite the oscillations. It is also pointed out that the problem
cannot be reduced to the 0D predator-prey model.

Finally, we demonstrate the effect of atom temperature in the solo regime. We
find that a finite energy spread of injected atoms strongly suppress the
oscillations amplitude compared to the cases of the injection with the same
velocity (called below as monokinetic injection). It is found that even a
small spread in atom temperature for the solo regime notably lowers the
amplitude of breathing mode. For the first regime, with the presence of
resistive modes, the atom temperature effect is negligible. We also show the
role of ion heating (due to the resistive modes) and selective ionization of
neutrals.

The paper is organized as follows. In Section II the detailed description of
both fluid and hybrid models is given. Section III defines the main features
of two distinct regimes of low-frequency oscillations and presents results for
three test cases with a detailed comparison between models.

##  II Fluid and hybrid models of low-frequency dynamics

Detailed description of the full fluid model and the hybrid model is presented
in this section. Two models share the same electron fluid equations (drift-
diffusion approximation and energy evolution). The models are considered in
the electrostatic and quasineutral approximation, with the three species:
neutral atoms, ions, and electrons. The ionization effects included via the
electron-atom collisions, serving as a mechanism for supporting plasma
discharge. Atom losses are only due to ionization effects, radial atom losses
were not included. The models also includes the self-consistent electric
field, the anode plasma recombination, the electron pressure effects, and the
electron heat flux across the magnetic field.

The simulated length of $5\text{\,}\mathrm{c}\mathrm{m}$ is assumed in the
axial direction of a Hall thruster ($x$-direction), with the channel exit in
the middle where the radial magnetic field has its maximum, Fig. 1. The
profile of the magnetic field amplitude given by
$B=B_{0}\exp{\left[-\left(x-x_{0}\right)^{2}/2\delta_{B}^{2}\right]}$, where
$x_{0}=$2.5\text{\,}\mathrm{c}\mathrm{m}$$ is the channel exit location, and
$\delta_{B}$ is the characteristic width coefficient for the magnetic field
profile, which are set
$\delta_{B,\text{in}}=$1.1\text{\,}\mathrm{c}\mathrm{m}$$,
$\delta_{B,\text{out}}=$1.8\text{\,}\mathrm{c}\mathrm{m}$$, respectively for
the inner and outer regionslandmark .

Figure 1: The magnetic field profile used in simulations, with the channel
exit located $2.5\text{\,}\mathrm{c}\mathrm{m}$ from anode (dashed line).

###  II.1 Fluid model

First, a short description of each species dynamics is given and then the full
system of time-dependent equations is formulated. For the neutral atoms, the
constant flow velocity $V_{a}$ along the channel is considered, and the
continuity equation with the source term is used to describe their dynamics:

$\frac{\partial n_{a}}{\partial t}+V_{a}\frac{\partial n_{a}}{\partial
x}=-\beta n_{a}n_{e},$ (1)

where $n_{a}$ is the atom number density, $n_{e}$ is the electron number
density, $\beta(\varepsilon)$ is the ionization rate coefficient that depends
on the electron energy $\varepsilon=(3/2)T_{e}$, where $T_{e}$ is the electron
temperature. The ionization due to electron-atom impact produces a pair of ion
and electron with a loss of neutral atom, hence the sink term $-\beta
n_{i}n_{e}$ (same and opposite sign source terms are included in the ion and
electron continuity equations). The ionization rate
$\beta(\varepsilon)=\langle\sigma(v)v\rangle$ is obtained via
BOLSIGhagelaar2005solving  by averaging over Maxwellian EEDF with the cross-
sections from the SIGLO databasesiglo .

The ion species are unmagnetized (gyroradius is much larger than the thruster
dimensions for a typical magnitude of a magnetic field in the thruster) and
described with the conservation of number density and momentum equations:

$\displaystyle\frac{\partial n_{i}}{\partial t}+\frac{\partial}{\partial
x}\left(n_{i}V_{i}\right)=\beta n_{a}n_{e},$ (2) $\displaystyle\frac{\partial
V_{i}}{\partial t}+V_{i}\frac{\partial V_{i}}{\partial
x}=\frac{e}{m_{i}}E+\beta n_{a}\left(V_{a}-V_{i}\right),$ (3)

where $n_{i}$ is the ion density, $V_{i}$ is the ion flow velocity, $E$ is the
axial electric field, $e$ is the elementary charge, $m_{i}$ is the ion mass
(Xenon, $131.293\text{\,}\mathrm{a}\mathrm{m}\mathrm{u}$). The ion pressure
term and generalized viscosity tensor are neglected in this model (ions are
ballistic with temperatures much lower than that of electron component).

The magnetized electron species are described with the first three fluid
moment equations: (electron inertia is neglected)

$\displaystyle\frac{\partial n_{e}}{\partial t}+\frac{\partial}{\partial
x}\left(n_{e}V_{ex}\right)=\beta n_{n}n_{e},$ (4) $\displaystyle
0=-\frac{e}{m_{e}}\mathbf{E}-\frac{eB}{m_{e}}\left(\mathbf{V}_{e\perp}\times\mathbf{\hat{z}}\right)-\frac{1}{n_{e}m_{e}}\frac{\partial\left(n_{e}T_{e}\right)}{\partial
x}-\nu_{m}\mathbf{V}_{e\perp},$ (5)
$\displaystyle\frac{3}{2}\frac{\partial}{\partial
t}\left(nT_{e}\right)+\frac{5}{2}\frac{\partial}{\partial
x}\left(n_{e}V_{ex}T_{e}\right)+\frac{\partial q_{e}}{\partial
x}=-n_{e}V_{ex}\frac{\partial\phi}{\partial
x}-n_{e}n_{a}\mathrm{K}-n\mathrm{W},$ (6)

where $n_{e}$ is the electron density,
$\mathbf{V}_{e\perp}=(V_{ex},V_{e\theta})$ is the electron flow velocity
perpendicular to the magnetic field ($x$ and $\theta$ are axial and azimuthal
coordinates), $m_{e}$ is electron mass, $B$ is the external magnetic field,
$\nu_{m}$ is the total electron momentum exchange frequency, $\mathrm{W}$ is
the anomalous energy loss coefficient, $\mathrm{K}$ is the collisional energy
loss coefficient, generated by BOLSIGhagelaar2005solving  as a table-valued
function, and $q_{e}$ is the electron heat flux. Phenomenological anomalous
electron energy loss coefficient (e.g., described the sheath radial energy
losses) $\mathrm{W}$ is introduced boeuf1998low  as

$\mathrm{W}=\nu_{\varepsilon}\varepsilon\exp{\left(-U/\varepsilon\right)},$
(7)

where $\varepsilon=3T_{e}/2$, $U=$20\text{\,}\mathrm{e}\mathrm{V}$$,
$\nu_{\varepsilon}$ is the anomalous losses frequency coefficient. The heat
flux across the magnetic field is

$q_{e}=-\frac{5}{2}\mu_{e}nT_{e}\frac{\partial T_{e}}{\partial x}.$ (8)

The electron momentum conservation equation (5) is simplified assuming no
pressure gradients nor equilibrium electric fields in other than axial
direction, then the axial electron velocity (denoted further as $V_{e}$) can
be expressed as

$V_{e}=-\mu_{e}E-\frac{\mu_{e}}{n_{e}}\frac{\partial(nT_{e})}{\partial x},$
(9)

where the electron mobility $\mu_{e}$ is the well-known classical electron
mobility across the magnetic field:

$\mu_{e}=\frac{e}{m_{e}\nu_{m}}\frac{1}{1+\omega_{ce}^{2}/\nu_{m}^{2}},$ (10)

where $\omega_{ce}=eB/m_{e}$ is the electron cyclotron frequency. Eq. (9)
commonly called the drift-diffusion equation. The model of electron transport
as based on the assumption of the following total electron momentum exchange
collision frequency:

$\nu_{m}=\nu_{en}+\nu_{walls}+\nu_{B},$ (11)

where the electron-neutral collision frequency $\nu_{en}$, electron-wall
collision frequency $\nu_{walls}$, and anomalous Bohm frequency $\nu_{B}$ are
given with

$\displaystyle\nu_{en}=k_{m}n_{a},$ (12)
$\displaystyle\nu_{walls}=\alpha${10}^{7}\text{\,}\mathrm{[}\mathrm{s}^{-1}\mathrm{]}$,$
(13) $\displaystyle\nu_{B}=\left(\beta_{a}/16\right)eB/m_{e}.$ (14)

where
$k_{m}=$2.5\text{\times}{10}^{-13}\text{\,}\mathrm{m}^{-3}\mathrm{s}^{-1}$$,
$\alpha$ and $\beta_{a}$ are free parameters. For the electron mobility model
different parameters are used inside and outside the channel (denoted
additionally as in, out): the near wall conductivity contribution
$\alpha_{\text{in}}=1,\ \alpha_{\text{out}}=0$, the anomalous contribution is
set to $\beta_{a,\text{in}}=0.1$, $\beta_{a,\text{out}}=1$.

Here we seek the low-frequency and bulk plasma modes, thus electron inertia is
neglected (as shown above) and further the full plasma quasineutrality is
assumed. One can quantifymakowski2001review  the quasineutral approximation
using the Poisson equation $\varepsilon_{0}\partial E/\partial
x=e(n_{i}-n_{e})$ ($\varepsilon_{0}$ is the permittivity of free space). With
the typical values of electric field
$E\approx${10}^{4}\text{\,}\mathrm{V}\mathrm{/}\mathrm{m}$$ and the size of
acceleration zone $1\text{\,}\mathrm{c}\mathrm{m}$, average plasma density
$n_{0}=${10}^{17}\text{\,}\mathrm{m}^{-3}$$, the difference
$(n_{i}-n_{e})/n_{0}\approx$5\text{\times}{10}^{-4}$$. Instead of the Poisson
equation, the electric field is found from the electron momentum equation as
shown below. Note that while the quasineutrality neglects a potential drop on
the Debye sheath near the anode, it still allows the presheath region to form
if the electron pressure is includedCohenZurPoP2002 ; AhedoPoP2005 ;
dorf2003anode ; dorf2005experimental , as in our electron model. The presheath
is the region where the electric field is induced to accelerate ions towards a
plasma boundary to compensate the electron current due to the pressure
gradient.

Thus, the full system of time-dependent fluid equations to be solved include
Eqs. (1,2,3,6), along with the drift-diffusion form for the electron velocity
Eq. (9). Full quasineutrality $n=n_{i}=n_{e}$ is enforced and the self-
consistent electric field is found via the electron drift-diffusion equation
(9), given by

$E=\frac{J_{T}}{en\mu_{e}}-\frac{V_{i}}{\mu_{e}}-\frac{1}{n}\frac{\partial
nT_{e}}{\partial x},$ (15)

where the total current density $J_{T}=en\left(V_{i}-V_{e}\right)$. Here
$J_{T}$ is constant in space (divergenceless current), which can be seen by
combining continuity equations for ions and electrons with quasineutral
assumption. We will use the integral approach, consisting of the evaluation
the total current density $J_{T}$ via the constraint $\int_{0}^{L}Edx=U_{0}$
(then it is substituted to Eq. (15) to evaluate the electric field), which
yields to

$J_{T}=\dfrac{U_{0}+\bigintss_{0}^{L}\left(\dfrac{V_{i}}{\mu_{e}}+\dfrac{1}{n}\dfrac{\partial
p_{e}}{\partial x}\right)dx}{\bigintss_{0}^{L}\dfrac{dx}{en\mu_{e}}},$ (16)

where $L$ is the system length, $U_{0}$ is the applied voltage (without sheath
voltage).

The fluid model is solved with the following boundary conditions. A constant
mass flow rate $\dot{m}$ and the full recombination of plasma that flows to
the anode determines the value of $n_{a}$ at the anode boundary,

$n_{a}(0)=\frac{\dot{m}}{m_{i}AV_{a}}-\frac{nV_{i}(0)}{V_{a}},$ (17)

where $A$ is the anode surface area. The value of the ion velocity is imposed
at the anode $V_{i}(0)=-b_{v}\sqrt{T_{e}/m_{i}}$ with the parameter
$b_{v}=0\text{--}1$, the Bohm velocity factor that can be varied. Both anode
and cathode electron temperature are fixed with
$T_{e}(0)=T_{e}(L)=$2\text{\,}\mathrm{e}\mathrm{V}$$.

As noted in Ref. makowski2001review, , plasma acceleration in the
configuration of the axial direction of a Hall thruster shows similarities to
the flow in de Laval nozzle. Indeed, the whole acceleration region can be
split into subsonic $V_{i}<c_{s}$ and supersonic $V_{i}>c_{s}$ regions, where
$c_{s}=\sqrt{T_{e}/m_{i}}$ is the ion sound speed. While in de Laval nozzle
the transition through sonic point happens at the region with the smallest
cross-section of the channel (due to extrema condition and regularity
requirement), for the bounded plasma configuration the position is determined
via the nonlinear relationship between plasma parameters (and their first
derivative) at the sonic point and a value of the total
currentsmolyakov2019stationary . Note that the total current $J_{T}$, given by
Eq. (16), is a function of $U_{0}$ and the integral dependence on all main
plasma parameters, thus the problem is inherently nonlocal, which has no
analogy with the standard Laval nozzle. Another difference is that in the
axial direction of Hall thrusters the presheath region can induce the backward
ion flow in a large portion of the thruster channel (see Case 2 below).

###  II.2 Hybrid model

The hybrid model has the same electron equations as in the fluid model, while
ions and neutrals are modeled via particle-in-cell (PIC) method
hagelaar2002two ; hagelaar2003role ; hagelaar2004modelling . The plasma
recombination effect is also included via the relationship (17), but the ion
velocity at the anode is not forced to satisfy the Bohm velocity. The
ionization is included via the electron-atom impact with the Monte Carlo
simulations (via the null collision methodbirdsall1991particle ;
hagelaar2008modelling ) with the known ionization rate $\beta(\varepsilon)$,
obtained in the same way as described above the fluid model. Neutral atoms are
injected with the constant flow rate $\dot{m}$ either with the constant
velocity $V_{a}$ (monokinetic), thus $f(v_{x})=\delta(v_{x}-V_{a})$, or with
the half-Maxwellian velocity distribution function at the left wall (anode),

$f(v_{x})=\frac{2v_{x}}{v_{Ta}^{2}}\exp({-\frac{v_{x}^{2}}{v_{Ta}^{2}}}),\ \
v_{x}>0,$ (18)

where $v_{Ta}^{2}=2k_{B}T_{a}/m_{a}$ is the atom thermal speed
($m_{a}=m_{i}$), $T_{a}$ is the atom temperature (in $\mathrm{K}$). For the
half-Maxwellian injection (18) the average flow velocity is
$v_{Ta}/\sqrt{\pi}$. Ions are assumed only singly charged and unmagnetized,
thus they are only accelerated by an electric field. Ions are produced
accordingly to the ionization rate coefficient, i.e. self-consistently with
the electron temperature evolution and local atom density. Both atoms and ions
are lost at boundaries. In this model the ion velocity is not forced to the
Bohm velocity, like in the fluid ion model. With the quasineutral approach
($n_{e}=n_{i}$ is forced at every time step), plasma density is evaluated from
the ion particle distribution and thereafter used for the electron temperature
(6) and the electric field (9) calculations. Formally, the evolution of the
distribution function for ions $f_{i}(x,v_{x},t)$ and atoms $f_{a}(x,v_{x},t)$
is described with the Boltzmann equation for each specie:

$\displaystyle\frac{\partial f_{i}}{\partial t}+v_{ix}\frac{\partial
f_{i}}{\partial x}+\frac{e}{m_{i}}E\frac{\partial f_{i}}{\partial
v}=S(x,v_{x}),$ (19) $\displaystyle\frac{\partial f_{a}}{\partial
t}+v_{ax}\frac{\partial f_{a}}{\partial x}=-S(x,v_{x}),$ (20)

where $S(x,v_{x})$ is the collisional source term due to the ionization. In
the case of the monokinetic target species (atoms), the ionization leads to
the ion creation with the atoms velocity $V_{a}$ (constant) and the source
term can be expressedboeuf1998low  as $S(x,v_{x})=\beta
n_{e}n_{a}\delta(v_{x}-V_{a})$. For the simulations with finite atom
temperature, newly created ions assigned velocities by sampling from the
isotropic Maxwellian distribution with the temperature $T_{a}$ with the
standard sampling techniquescartwright2000loading . Solutions to Eqs. (19,20)
effectively obtained by solving the motion equations for the corresponding
specie (the method of characteristics via PIC method). Eqs. (19,20) and the
electron fluid equations (6, 9) form the complete set of equations solved in
the hybrid model.

##  III Simulation results and comparison

The fluid and the hybrid models described above were studied for the three
test cases (denoted as Cases 1-3). Cases 1-2 will demonstrate two distinct
regimes of low-frequency oscillations, and Case 3 shows the effect of atom
temperature. Note that both Cases 1-2 use monokinetic atoms with the velocity
$v_{a}=$150\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$. Cases 1-2 are chosen
with the following observation: the larger values of
$\nu_{\varepsilon,\text{in}}$ allow the high frequency (of ion fly-by time)
modes to appear and coexist with the low-frequency modes, represented in Case
1. The simulations with lower value of $\nu_{\varepsilon,\text{in}}$ reveal
only the large amplitude low-frequency oscillations, we call it the solo
regime. In Case 1 the anomalous energy loss coefficient
$\nu_{\varepsilon,\text{in}}=$0.95\text{\times}{10}^{7}\text{\,}\mathrm{s}^{-1}$$,
and in Case 2
$\nu_{\varepsilon,\text{in}}=$0.4\text{\times}{10}^{7}\text{\,}\mathrm{s}^{-1}$$.
This is the only parameter distinguishing Cases 1 and 2. Note that for all
cases reported in this paper we keep
$\nu_{\varepsilon,\text{out}}=${10}^{7}\text{\,}\mathrm{s}^{-1}$$. Besides the
different time-dependent behaviour, these regimes show a notable difference in
the time-averaged axial profiles of the ion velocity and the electron energy,
see Figs. 2a,2b.

(a)

(b)

Figure 2: Averaged in time profiles of ion velocity (a) and electron
temperature (b) for Case 1
($\nu_{\varepsilon}=$0.4\text{\times}{10}^{7}\text{\,}\mathrm{s}^{-1}$$) and
Case 2
($\nu_{\varepsilon}=$0.95\text{\times}{10}^{7}\text{\,}\mathrm{s}^{-1}$$). The
result obtained with the hybrid model.

For both Cases the ion velocity profile is similar near the exit and beyond,
but in the near-anode region the ion backflow region (where ions are moving
towards the anode) is much shorter for Case 1 (Fig. 2a). The ion backflow
region is associated with the presheath formation near the anode (with
negative electric field). The size of the presheath region as a function of
$\nu_{\varepsilon,\text{in}}$ is shown in Fig. 4, where the transition between
regimes with large and short backflow region happens near the value
$\nu_{\varepsilon,\text{in}}=$0.75\text{\times}{10}^{7}\text{\,}\mathrm{s}^{-1}$$.
The electron temperature for Case 2 mostly concentrated near the channel exit
(Fig. 2b) with very low values in the near-anode region, and for Case 1 its
spread more uniformly and to the near-anode region, which might be due to the
higher electron flow velocity near the anode, see Fig. 3 for the time-averaged
electron velocity components from Eq. (9). Note substantially higher electron
current due to the pressure gradient in Case 1 for
$x<$0.4\text{\,}\mathrm{c}\mathrm{m}$$, which also results in a larger total
electron flow velocity. Indeed, the average gradient parameter
$L_{n}^{-1}=\partial_{x}n_{i}/n_{i}$ is about five times larger for Case 1.
Recall, that the presheath region in formed due to diffusive electron current
by inducing ion current (total current is conserved), generating negative
electric in this region. All features presented above clearly distinguish Case
1 and Case 2.

Figure 3: Electron flow velocity components, response to the electric field
and due to due for Case 1 and 2 (result obtained with the hybrid model).
Figure 4: Extent of the presheath zone, defined as the region with negative
electric field near the anode, as a function of anomalous electron energy
losses coefficient $\nu_{\varepsilon}$ (result obtained with the hybrid
model).

Finally, Case 3 demonstrates the effect of finite atom temperature, where
atoms injected to the system with the half-Maxwellian distribution, Eq. (18).
The main effect of atom temperature is found to be a significant reduction of
breathing mode amplitude to those observed in Case 2. Besides finite atom
temperature, all parameters for Case 3 are exactly the same as in Case 2, and
the atom temperature is set to $T_{a}=$500\text{\,}\mathrm{K}$$ (average
injected flow velocity $142\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$, close to
monokinetic $142\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$ used in Cases 1-2).

###  III.1 Case 1: Low electron energy losses; the co-existence of low and
high frequency modes

This case exhibits both low and higher frequency oscillations in the fluid and
the hybrid models. The hybrid model shows a smaller total current amplitude,
see Figs. 5a, 5b. The time-averaged total currents are close,
$8.2\text{\,}\mathrm{A}$ in the hybrid model, and $8.3\text{\,}\mathrm{A}$ in
the fluid model. The time-averaged ratio of the ion current (at the plume
exit, $x=$5\text{\,}\mathrm{c}\mathrm{m}$$) to the total current is 45% in the
hybrid model and 48% in the fluid.

(a)

(b)

Figure 5: Amplitudes of the total, ion, and electron currents in fluid model
(a) and hybrid model (b). Ion and electron currents are evaluated at
$x=$5\text{\,}\mathrm{c}\mathrm{m}$$. Small averaging window of length
$2.3\text{\,}\mu\mathrm{s}$ is applied to the ion and electron currents (to
filter out high-frequency noise).

Spectral power of the total current also shows some differences in both low-
and high-frequency range, Figs. 6a,6b. The main low-frequency mode in the
fluid model is $11.4\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$, while it is
$14.4\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ in the hybrid model. The total
current signal in the hybrid model contains more noise (statistical noise due
to the use of macroparticles), but the high-frequency component is clearly
seen at around $125\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$. In the fluid
model the high-frequency mode is larger, centered at about
$175\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$.

(a)

(b)

Figure 6: Spectral density of the total current yield in fluid model (a) and
hybrid model (b).

Besides the currents, a more rigorous comparison between the two models is
shown in terms of time-averaged axial profiles of various plasma quantities.
Due to the oscillatory nature of these solutions, the averaging time window
was chosen as the ten periods of the corresponding main low-frequency mode in
each simulation. The profiles depicted in Figs. 7a-7d. The main discrepancy
lies in the peak plasma density in the ionization (source) region at about
$1.3\text{\,}\mathrm{c}\mathrm{m}$ from the anode (higher value in the hybrid
model). Also, the ion velocity in the plume
($x>$2.5\text{\,}\mathrm{c}\mathrm{m}$$) is slightly higher in the hybrid
model.

(a)

(b)

(c)

(d)

Figure 7: Comparison of time-averaged axial macroscopic profiles resulted from
fluid and hybrid models of neutral density (a), plasma density (b), ion flow
velocity (c), and electron temperature (d).

The ion phase space (hybrid model) is shown in Fig. 8a. The ion velocity
distribution function (IVDF) is highly inhomogeneous inside the channel,
suggesting that the higher fluid moments may play a role (recall, that the
closure in the ion fluid model is given by nullified pressure term).
Typically, ion pressure effects are neglected due to a low ion temperature;
ion formation is due to the ionization process as they carry the low atom
temperature even accelerating by the axial electric field.

(a)

(b)

Figure 8: Instantaneous image of ion distribution function (in space of axial
coordinate and axial velocity) in the hybrid model (a). Ion temperature
spatial profile (time-averaged) evaluated from ion kinetic representation in
the hybrid model (b).

To check the validity of this assumption for this case the fluid moments were
calculated from the kinetic particle representation in the hybrid model. It
allows to test the ion momentum balance equation (to identify the role of the
ion pressure term) in a more complete form

$\frac{\partial V_{i}}{\partial t}+V_{i}\frac{\partial V_{i}}{\partial
x}=\frac{e}{m_{i}}E-\frac{1}{n_{i}}\frac{\partial p_{i}}{\partial x}+\beta
n_{a}\left(V_{a}-V_{i}\right),$ (21)

where the fluid moments, such as ion density $n_{i}$, ion flow velocity
$V_{i}$, and ion pressure $p_{i}$ are evaluated from the ion distribution
function in the following way:

$\displaystyle n_{i}=\int f_{i}dv_{ix},$ (22) $\displaystyle
V_{i}=\frac{1}{n_{i}}m_{i}\int v_{x}f_{i}dv_{ix},$ (23) $\displaystyle
p_{i}=m_{i}\int v_{i}^{\prime}v_{i}^{\prime}f_{i}dv_{ix},$ (24)

$v_{i}^{\prime}=v_{i}-V_{i}$ is the random component of particle velocity, and
$V_{i}$ is the average (flow) velocity. The time-averaged profile of the ion
temperature $T_{i}=p_{i}/n_{i}$ is shown in Fig. 8b, revealing higher than
typical values of the ion temperature, with the average over the whole domain
of $1.7\text{\,}\mathrm{e}\mathrm{V}$. For the momentum balance test, each
term in Eq. (21) was evaluated as a function of time and space and then
averaged in time over a few periods of the main low-frequency mode. Fig. 9a
shows the difference of left- and right-hand sides of equation (21) that
ideally must add up to zero. The plotted terms were normalized to the value
$V^{2}/L$ where the ballistic ion velocity is $V^{2}=2eU_{0}/m_{i}$ with the
potential difference $U_{0}=$300\text{\,}\mathrm{V}$$ over the system length
$L$. It is seen that the ion pressure term notably improves the overall ion
momentum balance, suggesting the fluid model for this configuration should not
ignore the ion pressure.

(a)

(b)

Figure 9: Difference between left and right hands sides of the ion momentum
balance equation (ion-mom) evaluated from the kinetic description of the
hybrid model (a). Separate terms of the same equation averaged in time
(denoted with angle brackets) over few periods of the low-frequency component
(b).

At this moment, it is interesting to further inspect each term in the Eq.
(21). The unsteady term $\partial V_{i}/\partial t$ is negligible due to time
averaging, while the rest four terms have comparable values inside the
channel, see Fig 9b. In the near-anode region, the ion pressure is negligible
($T_{i}\approx$0.1\text{\,}\mathrm{e}\mathrm{V}$$), and no ionization, so the
ions accelerate towards the anode in a weak negative electric field
ballistically. Then in the ionization (source) region, we see all terms are
comparable. The ion pressure and the collisional drag compensate the ballistic
acceleration, so the ion inertia remains low. Finally, to the right in the
acceleration zone ($x>$1.5\text{\,}\mathrm{c}\mathrm{m}$$) the inertial and
ballistic terms start to dominate (which continues outside the channel). It
can be seen that the ion pressure term changes the sign at
$x\approx$1.6\text{\,}\mathrm{c}\mathrm{m}$$ (due to ion density profile, Fig.
7b) and contributes to the ion acceleration. At the same time, the collisional
drag continues to slow down ions.

Based on the results above, the ion pressure force term was added to the fluid
model with the temperature kept constant for simplicity. We understand that
for the self-consistent treatment, the ion energy evolution shall be included
or, at least, equation of state. We used
$T_{i}=$1.2\text{\,}\mathrm{e}\mathrm{V}$$ and the pressure $p_{i}=nT_{i}$
(average in Fig 6b is 1.7 eV, but fluid might result in stationary result).
This results in a better agreement between two models, the total current
amplitude decrease and the main low-frequency mode increase, Figs. 10a,10b.
The low-frequency mode increased to
$13.9\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$, and the high-frequency peak is
shifted to a lower value of $153\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$. The
ratio of the ion current to the total current shows improved agreement, 45%,
the same as in the hybrid model. However, the time-averaged total current
value of $8.7\text{\,}\mathrm{A}$ became somehow larger than in the hybrid
model (8.2 A). Also, the higher frequency component is closer to the hybrid
result. As for the peak plasma density discrepancy (Fig. 7b), it remained
similar.

(a)

(b)

Figure 10: Comparison of total current in the fluid model and hybrid model (a)
and the spectral power density of the total current in the fluid model (b).
The result obtained is for the fluid model model with the ion pressure term
included.

###  III.2 Case 2: High electron energy losses; the solo regime of the low
frequency mode

This case is subject to the low-frequency oscillations only, with the only
difference to Case 1 in the value of the anomalous electron energy loss
coefficient, which is
$\nu_{\varepsilon,\text{in}}=${10}^{7}\text{\,}\mathrm{s}^{-1}$$. As in the
previous case, the total current amplitude is higher in the fluid model ,
Figs. 11a, 11b, but the main oscillation frequency in two models is similar,
Figs. 12a, 12b.

(a)

(b)

Figure 11: Amplitudes of total, ion, and electron currents in fluid model (a)
and hybrid model (b) for Case 2. Ion and electron currents are evaluated at
$x=$5\text{\,}\mathrm{c}\mathrm{m}$$.

(a)

(b)

Figure 12: Spectral density of the total current for fluid model (a) and
hybrid model (b).

Time-averaged profiles well agree in two models, with the only notable
discrepancy in the plasma density profile, which is ${\sim}30\%$ lower inside
the channel region for the fluid model, see Figs. 13a-13d.

(a)

(b)

(c)

(d)

Figure 13: Comparison of spatial distribution of time-averaged macroscopic
profiles between fluid and hybrid models for Case 2.

Unlike in Case 1, the IVDF for this case (Fig. 14a) reveals that the ion
population remains cold (everywhere except in the near-anode region), and the
ion momentum balance must be well satisfied without the ion pressure term.
This is seen in Fig. 14b, where the ion pressure force remains low everywhere
in the channel. In fact, in this case, the ballistic ion acceleration is more
pronounced, dominating everywhere except the stall point ($V_{i}\approx 0$).
The average ion temperature, calculated similarly as in Case 1, is
$0.3\text{\,}\mathrm{e}\mathrm{V}$ and not exceeding
$1\text{\,}\mathrm{e}\mathrm{V}$ in the domain (about 5 times smaller than in
Case 1).

The main stages of the breathing mode dynamics are illustrated with plasma and
atom densities, Fig. 15. After atoms reach the ionization zone and undergo
ionization, plasma density increases and quickly depleted
(${\sim}$1\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$) to the left
(due to the backflow region with negative velocity) and the right of
ionization zone. Then all of the ion flux that reached the anode recombines
and form the peak in neutral density at the anode, increasing the number of
atoms advect to the ionization zone; this process repeats. In our previous
work, Ref. chapurin2021mechanism, , this setup was studied in detail, and it
was shown that the ion backflow region plays a crucial role in these
oscillations and that they can exist without the recombination mechanism.

(a)

(b)

Figure 14: Instantaneous image of ion distribution function (in space of axial
coordinate and axial velocity) in the hybrid model (a). Separate terms of the
ion momentum balance equation (22) evaluated from the ion kinetic description
(hybrid model), averaged in time (denoted with angle brackets) over few
periods of the low-frequency component (b). (a)  Figure 15: Neutral and plasma
density evolution during one oscillation period. Dashed line separates the
region with negative (to the left) and positive (to the right) ion velocity.
Figure 16: Minimum and maximum total current values during oscillations for
various ion velocities at the anode, expressed as fractions of the Bohm
velocity. Note that oscillations are absent for $V_{i}<-0.65c_{s}$.

It is found that the oscillation dynamics in Case 2 highly depend on the
recombination at the anode, modeled with Eq. (17), and the number of
recombined atoms proportional to the oscillation current amplitude. Note that
in Case 1, recombination plays a minor role and does not affect the presented
results. Turning off the recombination in Case 2 nullifies the oscillations
and the stationary solutions obtained (in both models). We suggest that the
main difference in this case between the fluid and the hybrid model, in this
case, lies in the differences with the ion velocity boundary condition.
Normally, in the fluid quasineutral models, this boundary condition is fixed
to the Bohm velocity, as in our fluid model. When the Bohm boundary condition
is scaled with the factor $b_{v}=0\textrm{--}1$ (ion velocity effectively
decreased at the anode) in the fluid model, the oscillations amplitude also
decreases, see Fig. 16. It clearly shows that the plasma recombination
provides the additional feedback in this configuration, and the main
difference with the hybrid model lies in the ion boundary conditions causing
larger oscillation amplitude in the fluid model. In the hybrid model (in our
setup), ions are not forced to satisfy the Bohm condition, and the flow
velocity is established self-consistently. We modeled this behaviour with the
free boundary condition for the ion velocity at the anode in the fluid model
(forcing $\partial_{x}^{2}V_{i}(0)=0$). In the presence of ion backflow (like
we see in Case 1,2) in the quasineutral approximation, the ion velocity at the
anode is defined self-consistently by the characteristic flowing from
$v_{i}=0$ inside the channel; thus, a fixed boundary is not required.

The obtained result in the fluid model with the free boundary for the ion
velocity at the anode reveals less violent oscillations and generally better
agrees with the hybrid model results. The total current oscillation amplitude
is lower, Fig. 17, closer to the hybrid mode result. Time-averaged atom and
plasma density profiles agreement is improved, Figs. 18a, 18b.

Figure 17: Comparison of the total current obtained in both models where the
free boundary condition used in the fluid model for ion velocity at the anode.

(a)

(b)

Figure 18: Atom (a) and ion (b) axial profiles, averaged in time, compared in
both models where the free boundary condition used in the fluid model.

###  III.3 Case 3: Effects of finite temperature of neutral atoms

Unlike Case 2, the oscillation amplitude in this case is much smaller in the
hybrid model, Fig. 20a, with the amplitude of ${\sim}$1\text{\,}\mathrm{A}$$
(in comparison to Case 2 with ${\sim}$13\text{\,}\mathrm{A}$$). In the
previous cases, injected atoms in the hybrid model were kept monokinetic (zero
thermal energy). Here, the finite atom temperature is included while all other
parameters are kept as in Case 2. Atoms are injected at the anode with the
half-Maxwellian flux distribution, Eq. (18) with a temperature
$T_{a}=$500\text{\,}\mathrm{K}$$. The average injection velocity for the half-
Maxwellian is
$V_{0}=v_{Ta}/\sqrt{\pi}=$142\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$, which
is close to $150\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$ used in the previous
monokinetic runs. It was noticed that the atom flow velocity in the hybrid
model exhibits a clear spatial dependence, “accelerating” along the channel.
This effect is known as selective ionization, observed both in experiments
meezan2001anomalous ; hargus2002interior  and simulationshara2012one . It is
clear that atoms in the fluid model with the simple advection equation (1)
cannot capture this effect; hence the atom momentum balance equation was
included in the fluid model:

$\frac{\partial\left(n_{a}V_{a}\right)}{\partial t}+\frac{\partial}{\partial
x}\left(n_{a}V_{a}^{2}\right)=-\beta
n_{a}n_{i}V_{a}-\frac{T_{a}}{m_{i}}\frac{\partial\left(n_{a}\right)}{\partial
x},$ (25)

where the closure is given with the constant temperature
$T_{a}=$500\text{\,}\mathrm{K}$$. The illustration that the momentum balance
given by Eq. (25) has a sufficient number of terms, the atom fluid moments
were calculated from the kinetic representation of atoms in the hybrid model,
and the balance is compared, see Fig. 19. Thus, Eq. (25) along with the
continuity equation

$\frac{\partial n_{a}}{\partial
t}+\frac{\partial\left(n_{a}V_{a}\right)}{\partial x}=-\beta n_{a}n_{i}$ (26)

solved in the fluid model, with $T_{a}=$500\text{\,}\mathrm{K}$$ and the fixed
atom flow velocity at the anode
$V_{a}(0)=$142\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$ (Dirichlet condition),
and the same recombination boundary condition given by Eq. (17). This allowed
to recover qualitatively the atom flow velocity behaviour, Fig. 20b, but in
the same time the fluid model resulted in a completely stationary solution,
Fig. 20a.

(a)  Figure 19: Atom momentum balance equation terms (a). Note that the number
of macroparticles (and atom density) significantly decrease to the right along
the channel elevating the noise.

(a)

(b)

Figure 20: The total current in fluid and hybrid simulations (a); comparison
of time-averaged spatial profile of atom flow velocity for both models (b).
Figure 21: Power spectrum of the total current in the hybrid simulation. Note
that there is increase in frequency (from 10.2 kHz to 18.2 kHz) in compare to
Case 2 (with monokinetic atoms).

To identify better the role of the atom distribution in the hybrid model, we
configured the shifted Maxwellian distributions for atoms with various
temperatures denoted as $T_{a,s}$, and the fixed shift
$150\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$. Thus, the value of
$T_{a,s}=$0\text{\,}\mathrm{K}$$ corresponds to Case 2. We found that even a
low spread for atom velocities, e.g. corresponding to
$T_{a,s}\approx$50\text{\,}\mathrm{K}$$, heavily damps the oscillation
amplitude to the order $1\text{\,}\mathrm{A}$, similar to what we obtained
with the half-Maxwellian and $T_{a}=$500\text{\,}\mathrm{K}$$.

Interestingly, we were not able to find higher amplitude low-frequency
oscillations for this case in the hybrid model (simply scanning parameters
like $\nu_{\varepsilon}$ or $\beta_{a}$). However, slightly varying the
position of the maximum magnetic field profile, higher amplitude breathing
modes were observed even with thermal atoms. This important behaviour needs
further attention which is out of the scope of this paper and is left for
future studies.

##  IV Summary

Generally, plasma models for Hall thruster configuration have to deal with
long time scales (defined by slow neutral dynamics, requiring at least
${\sim}10^{6}$ electron plasma periods) and large spatial scales (above
${\sim}10^{3}\leavevmode\nobreak\ \lambda_{De}$). A scale separation by many
orders is a typical problem in plasma modeling. Full kinetic descriptions for
the problem involving low-frequency dynamics in the axial direction may not be
practical due to computational cost and potential numerical problems. Reduced
plasma models, such as the fluid model formulated via the conservation laws
derived from higher dimensional kinetic formulation, are less computationally
expensive and capture main physical phenomena in many situations. Fluid
formulations typically allow easier analysis and physical interpretation of
the results. A good compromise is achieved with hybrid models. One species is
modeled with fluid equations, and the other is kinetic (e.g. fluid electrons
with neglected inertia are common for such low-frequency models). Typically,
full kinetic simulations for plasma thrusters require some speed-up
techniques, such as artificial increase of dielectric permittivity, an
increase of electron to ion mass ratio, geometrical downscaling szabo2001fully
; taccogna2005plasma ; liu2010pic ; charoy2020numerical , all to overcome
numerical constraints and resolve low-frequency dynamics. Full kinetic
(particle-in-cell) 2D axial-azimuthal model (with the artificial increase of
dielectric permittivity) reproduces the low-frequency ionization oscillations
with no assumptions made on anomalous electron transport coefficients
(electron current develops self-consistently via azimuthal drift
instabilities), presented in Ref. coche2014two, . Note that the radial
coordinate may introduce important modifications to the azimuthal drift mode
(excitation of large scale MTSI modes) with an influence on the anomalous
electron transportjanhunen2018evolution ; villafana2021rz ; jimenez2021twod .
Thus, ideally, the detailed model of electron transport will be reproduced in
a 3D simulation.

Low-frequency plasma dynamics in the axial direction of a Hall thruster are
studied with full fluid and hybrid (kinetic/fluid) models. The model
parameters are taken close to those of the LANDMARK benchmarking
projectlandmark . We identified and distinguished two low-frequency
oscillations regimes: one is subject to the low-frequency oscillations only,
and another with both low and higher (time of ion flyby oscillations)
frequency components. Another distinct feature is the extent of the ion
backflow region: in the case with the pure low-frequency component, the ion
backflow on average is about half of the channel length, while for the regime
where two modes coexist, it is much shorter, Fig. 2a. These regimes are
observed in fluid and hybrid models and presented as Cases 1,2 in the
numerical experiments. Finally, Case 3 illustrates the effects of finite atom
temperature.

The main difference between the full fluid and hybrid model results, as
expected, is caused by the different approaches in heavy species modeling. The
kinetic method for heavy particles in the hybrid model automatically includes
effects of the non-equilibrium distribution function (i.e. higher fluid
moments effects, only limited by statistical of the noise of the PIC
approach). However, as seen from the results, a simple fluid model for both
ion and atom components (first two moments) is well sufficient to reproduce
the main results and capture breathing mode frequency and the resistive mode
presence.

It is found that in Case 1 (regime with co-existence of low and higher
frequency modes), the finite value of ion pressure played a notable role in
oscillatory behaviour. Ion temperature is often assumed negligible in Hall
thrusters operation (similar to a neutral temperature, up to about
$1000\text{\,}\mathrm{K}$ or $0.1\text{\,}\mathrm{e}\mathrm{V}$), thus it was
not included in the primary fluid model. Using kinetic representation from the
hybrid model, we show that the ion momentum balance without the ion pressure
was notably violated. Thus the ion pressure (with constant temperature) term
was added into the fluid model. It resulted in better conformity between the
two models. Nevertheless, both fluid and hybrid models predicted the existence
of two modes in this configuration even without the ion pressure in the fluid
model, suggesting that these oscillations are the so-called resistive modes
chable2005numerical ; koshkarov2018current . We believe these modes are
observed in the models without electron pressure, many of which were done in
the primary studies and modeling of axial dynamics of Hall
thrustersmorozov2000fundamentals ; chable2005numerical . We conjecture that
the resistive modes of higher frequencies play an important role in the
excitation of the low-frequency modes.

For Case 2, which shows only the low-frequency oscillations, we already
performed extensive studies, Ref. chapurin2021mechanism, , where we identified
that the mechanism of these oscillations lies in the ion backflow region
(presheath) and that they can be additionally reinforced with the plasma
recombination. This led to the observation that the difference in the ion
boundary conditions at the anode played a crucial role in our benchmark. Bohm
velocity for ions in the fluid model generated a larger atom yield (due to
recombination, Eq. (17)) from the anode in comparison to the hybrid model
where this velocity was unconstrained and found to be smaller on average (thus
leading to smaller atom yield). Modifying the ion velocity boundary condition
in the fluid model (replacing the Bohm condition with a free boundary
condition) resulted in better conformity between the two models. Similar to
the previous case, we conclude that the main physical behaviour was identified
in both models even without this modification.

Finally, in Case 3, with finite atom temperature in the hybrid model but the
same values of other parameters as in Case 2, it is shown that the advection
equation in the fluid model, Eq. (1), with the constant flow velocity is not
sufficient to describe atom dynamics. An effect of selective ionization of
neutral particles is observed in this case, where average macroscopic velocity
increased more than twice along the channel. Thus, the fluid model was
supplemented with the atom momentum balance equation (25), which resulted in a
better agreement between the two models, reproducing the selective ionization
effect. Though it resulted in a completely stationary solution, the hybrid
model preserved the breathing oscillations (with similar frequency to Case 2)
but with much smaller amplitude.

In summary, benchmarking of the fluid and hybrid models show very close
agreement for averaged plasma parameters profiles. Both models reveal the
existence of the two different regimes of the low-frequency oscillations in
Hall thrusters. While qualitatively, the two regimes are identifiable in both
models, there are some quantitative differences in the frequencies and the
amplitude of the oscillations. These differences are attributed to the ion
finite thermal (pressure) effects which were not originally included in the
fluid model. Account of the finite ion pressure improves the agreement. The
finite temperature (energy) spread of the neutral atoms provides a strong
stabilizing effect on the oscillations. These results highlight limitations of
the fluid models that have to be considered in future modeling of practical
devices.

##  Appendix A Resistive modes

Additionally to Case 1 with the coexistence of both breathing and resistive
modes with relatively small amplitudes, we present cases with more prominent
resistive modes and show their effect on ion heating. Larger value of electron
anomalous energy losses,
$\nu_{\varepsilon,\text{in}}=$1.2\text{\times}{10}^{7}\text{\,}\mathrm{s}^{-1}$$,
leads to a solution where restive mode dominate, see Fig. 23a. It reveals
larger amplitude and clearly distinct higher-frequency oscillations (168 kHz),
with a small low-frequency modulation (14.4 kHz), Fig. 23b its spectrum. As it
was demonstrated for Case 1, ion temperature effects were not negligible, and
ions were heated to a few electron volts. We noticed that the average ion
temperature is higher with presence of the resistive modes as shown in Fig.
22, where the resistive modes appear for values
$\nu_{\varepsilon,\text{in}}=$0.9\text{\,}\mathrm{s}^{-1}$$ and higher.

Figure 22: Ion average thermal energy as a function of the anomalous electron
energy loss coefficient. As we move into regime with resistive modes present
(studied in Case 1), ion heating increases and ion pressure effects may become
important.

(a)

(b)

Figure 23: Currents (a) and the total current spectra (b).

Now, we show that the frequency of the resistive mode scales with the average
electron collision frequency $\bar{\nu}_{m}$, Fig. 24a, and thus higher
electron mobility ${\sim}1/\bar{\nu}_{m}$ leads to higher frequencies. At the
same time, the ion heating effect is stronger for the higher frequency of
resistive modes, Fig. 24b. It is interesting to note that the breathing mode
frequency stays approximately the same (Fig. 24a), along with the size of the
presheath region.

(a)

(b)

Figure 24: Frequency of breathing and resistive modes as a function of
averaged (over time and space inside the channel) total electron momentum
exchange frequency $\bar{\nu}_{m}$ (a). Average ion temperature (in time and
whole domain) as function of $\bar{\nu}_{m}$ (b). The parameter varied in this
study was the anomalous Bohm coefficient $\beta_{a,in}$ inside the channel,
directly affecting the total electron frequency, Eq. (11).

###### Acknowledgements.

This work is supported in part by NSERC Canada, and US Air Force Office of
Scientific Research FA9550-15-1-0226.

## Data availability

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

## References

  * [1] A.I. Morozov and V.V. Savelyev.  Fundamentals of stationary plasma thruster theory.  In Reviews of plasma physics, pages 203–391. Springer, 2000.
  * [2] A.I. Smolyakov, O. Chapurin, W. Frias, O. Koshkarov, I. Romadanov, T. Tang, M. Umansky, Y. Raitses, I.D. Kaganovich, and V.P. Lakhin.  Fluid theory and simulations of instabilities, turbulent transport and coherent structures in partially-magnetized plasmas of discharges.  Plasma Physics and Controlled Fusion, 59(1):014041, 2016.
  * [3] T. Lafleur, S. D. Baalrud, and P. Chabert.  Theory for the anomalous electron transport in Hall effect thrusters. i. insights from particle-in-cell simulations.  Physics of Plasmas, 23(5):053502, 2016.
  * [4] S. Janhunen, A. Smolyakov, O. Chapurin, D. Sydorenko, I. Kaganovich, and Y. Raitses.  Nonlinear structures and anomalous transport in partially magnetized E$\times$B plasmas.  Physics of Plasmas, 25(1):011608, 2018.
  * [5] S. Janhunen, A. Smolyakov, D. Sydorenko, M. Jimenez, I. Kaganovich, and Y. Raitses.  Evolution of the electron cyclotron drift instability in two-dimensions.  Physics of Plasmas, 25(8):082308, 2018.
  * [6] W. Villafana, F. Petronio, A. C. Denig, M. J. Jimenez, D. Eremin, L. Garrigues, F. Taccogna, A. Alvarez-Laguna, J. P. Boeuf, A. Bourdon, P. Chabert, T. Charoy, B. Cuenot, K. Hara, F. Pechereau, A. Smolyakov, D. Sydorenko, A. Tavant, and O. Vermorel.  2D radial-azimuthal particle-in-cell benchmark for E$\times$B discharges.  Plasma Sources Science and Technology, 30(7):075002, jul 2021.
  * [7] M. Jimenez.  2D3V Particle-in-cell simulations of Electron Cyclotron Drift Instability and anomalous electron transport in E$\times$B plasmas.  Master’s thesis, University of Saskatchewan, 2021.
  * [8] I.D. Kaganovich, A. Smolyakov, Y. Raitses, E. Ahedo, I.G. Mikellides, B. Jorns, F. Taccogna, R. Gueroult, S. Tsikata, A. Bourdon, J.-P. Boeuf, M. Keidar, A. T. Powis, M. Merino, M. Cappelli, K. Hara, J. A. Carlsson, N. J. Fisch, P. Chabert, I. Schweigert, T. Lafleur, K. Matyash, A.V. Khrabrov, R. W. Boswell, and A. Fruchtman.  Physics of E$\times$B discharges relevant to plasma propulsion and similar technologies.  Physics of Plasmas, 27(12):120601, 2020.
  * [9] Y.B. Esipchuk, A.I. Morozov, G.N. Tilinin, and A.V. Trofimov.  Plasma oscillations in closed-drift accelerators with an extended acceleration zone.  Soviet Physics Technical Physics, 18:928, 1974.
  * [10] J. Fife, M. Martinez-Sanchez, J. Szabo, J. Fife, M. Martinez-Sanchez, and J. Szabo.  A numerical study of low-frequency discharge oscillations in Hall thrusters, page 3052.  2012\.
  * [11] S. Barral and E. Ahedo.  On the origin of low frequency oscillations in Hall thrusters.  AIP Conference Proceedings, 993(1):439–442, 2008.
  * [12] Kentaro Hara, Michael J Sekerak, Iain D Boyd, and Alec D Gallimore.  Perturbation analysis of ionization oscillations in Hall effect thrusters.  Physics of Plasmas, 21(12):122103, 2014.
  * [13] O. Chapurin, A.I. Smolyakov, G. Hagelaar, and Y. Raitses.  On the mechanism of ionization oscillations in Hall thrusters.  Journal of Applied Physics, 129(23):233307, 2021.
  * [14] T. Lafleur, P. Chabert, and A. Bourdon.  The origin of the breathing mode in hall thrusters and its stabilization.  Journal of Applied Physics, 130(5):053305, 2021.
  * [15] S. Chable and F. Rogier.  Numerical investigation and modeling of stationary plasma thruster low frequency oscillations.  Physics of plasmas, 12(3):033504, 2005.
  * [16] K. Makowski, Z. Peradzynski, S. Barral, and M.A. Dudeck.  Review of the plasma fluid models in stationary plasma thrusters.  High Temperature Material Processes: An International Quarterly of High-Technology Plasma Processes, 5(2), 2001.
  * [17] G.J.M. Hagelaar, J. Bareilles, L. Garrigues, and J.-P. Boeuf.  Modelling of stationary plasma thrusters.  Contributions to Plasma Physics, 44(5-6):529–535, 2004.
  * [18] J.-P. Boeuf and L. Garrigues.  Low frequency oscillations in a stationary plasma thruster.  Journal of Applied Physics, 84(7):3541Hall3554, 1998.
  * [19] S. Barral, Z. Peradzynski, K. Makowski, and M.A. Dudeck.  Fluid model of Hall thruster – comparison with hybrid model.  High Temperature Material Processes: An International Quarterly of High-Technology Plasma Processes, 5(2), 2001.
  * [20] S. Barral and E. Ahedo.  Low-frequency model of breathing oscillations in Hall discharges.  Physical Review E, 79(4):046401, 2009.
  * [21] V. Giannetti, M.M. Saravia, L. Leporini, S. Camarri, and T. Andreussi.  Numerical and experimental investigation of longitudinal oscillations in Hall thrusters.  Aerospace, 8(6):148, 2021.
  * [22] A. I. Morozov and V. V. Savel’ev.  One-dimensional hybrid model of a stationary plasma thruster.  Plasma Physics Reports, 26(10):875–880, 2000.
  * [23] Andrey Shashkov, Alexander Lovtsov, and Dmitry Tomilin.  A one-dimensional with three-dimensional velocity space hybrid-pic model of the discharge plasma in a hall thruster.  Physics of Plasmas, 24(4):043501, 2017.
  * [24] Gavrikov M.B. and Taiurskii A.A.  Hybrid model of a stationary plasma thruster (in Russian).  Preprints Keldysh Institute of Applied Mathematics, (0):35–47, 2021\.
  * [25] B.D. Dudson, M.V. Umansky, X.Q. Xu, P.B. Snyder, and H.R. Wilson.  Bout++: A framework for parallel plasma fluid simulations.  Computer Physics Communications, 180(9):1467Hall1480, 2009.
  * [26] G.J.M. Hagelaar, J. Bareilles, L.B.J.P. Garrigues, and J.-P. Boeuf.  Two-dimensional model of a stationary plasma thruster.  Journal of Applied Physics, 91(9):5592–5598, 2002.
  * [27] G.J.M. Hagelaar, J. Bareilles, L. Garrigues, and J.-P. Boeuf.  Role of anomalous electron transport in a stationary plasma thruster simulation.  Journal of Applied Physics, 93(1):67–75, 2003.
  * [28] Low temperAture magNetizeD plasMA benchmaRKs.  https://www.landmark-plasma.com/.
  * [29] A. A. Litvak and Nathaniel J. Fisch.  Resistive instabilities in Hall current plasma discharge.  Physics of Plasmas, 8(2):648–651, 2001.
  * [30] E. Fernandez, M. K. Scharfe, C. A. Thomas, N. Gascon, and M. A. Cappelli.  Growth of resistive instabilities in E$\times$B plasma discharge simulations.  Physics of Plasmas, 15(1):012102, 2008.
  * [31] O. Koshkarov, A.I. Smolyakov, I.V. Romadanov, O. Chapurin, M.V. Umansky, Y. Raitses, and I.D. Kaganovich.  Current flow instability and nonlinear structures in dissipative two-fluid plasmas.  Physics of Plasmas, 25(1):011604, 2018.
  * [32] A. I. Morozov and V. V. Savel’ev.  One-dimensional hydrodynamic model of the atom and ion dynamics in a stationary plasma thruster.  Plasma Physics Reports, 26(3):219–224, 2000.
  * [33] Gavrikov M.B. and Taiurskii A.A.  Some mathematical questions of plasma ionization (in Russian).  Preprints Keldysh Institute of Applied Mathematics, (94):1–48, 2021\.
  * [34] G.J.M. Hagelaar and L.C. Pitchford.  Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models.  Plasma Sources Science and Technology, 14(4):722, 2005.
  * [35] SIGLO database.  www.lxcat.net/SIGLO.  [Retrived: June-2013].
  * [36] A. Cohen-Zur, A. Fruchtman, J. Ashkenazy, and A. Gany.  Analysis of the steady-state axial flow in the Hall thruster.  Physics of Plasmas, 9(10):4363–4374, 2002.
  * [37] E. Ahedo and J. Rus.  Vanishing of the negative anode sheath in a Hall thruster.  Journal of Applied Physics, 98(4):043306, 2005.
  * [38] L. Dorf, V. Semenov, and Y. Raitses.  Anode sheath in Hall thrusters.  Applied Physics Letters, 83(13):2551–2553, 2003.
  * [39] L. Dorf, Y. Raitses, and N. J. Fisch.  Experimental studies of anode sheath phenomena in a Hall thruster discharge.  Journal of Applied Physics, 97(10):103309, 2005.
  * [40] A. Smolyakov, O. Chapurin, I. Romadanov, Y. Raitses, I. Kaganovich, G. Hagelaar, and J.-P. Boeuf.  Stationary profiles and axial mode oscillations in Hall thrusters.  In AIAA Propulsion and Energy 2019 Forum, pages 4080–4080, 2019\.
  * [41] C.K. Birdsall.  Particle-in-cell charged-particle simulations, plus Monte Carlo collisions with neutral atoms, PIC-MCC.  IEEE Transactions on plasma science, 19(2):65–85, 1991.
  * [42] G.J.M. Hagelaar.  Modelling methods for low-temperature plasmas.  PhD thesis, Université Toulouse III Paul Sabatier (UT3 Paul Sabatier), 2008.
  * [43] K.L. Cartwright, J.P. Verboncoeur, and C.K. Birdsall.  Loading and injection of Maxwellian distributions in particle simulations.  Journal of Computational Physics, 162(2):483–513, 2000.
  * [44] N.B. Meezan, W.A. Hargus Jr, and M.A. Cappelli.  Anomalous electron mobility in a coaxial Hall discharge plasma.  Physical Review E, 63(2):026410, 2001.
  * [45] W.A. Hargus and M.A. Cappelli.  Interior and exterior laser-induced fluorescence and plasma measurements within a Hall thruster.  Journal of Propulsion and Power, 18(1):159–168, 2002.
  * [46] K. Hara, I.D. Boyd, and V.I. Kolobov.  One-dimensional hybrid-direct kinetic simulation of the discharge plasma in a Hall thruster.  Physics of Plasmas, 19(11):113508, 2012.
  * [47] J.J. Szabo.  Fully kinetic numerical modeling of a plasma thruster.  PhD thesis, Massachusetts Institute of Technology, 2001.
  * [48] F. Taccogna, S. Longo, M. Capitelli, and R. Schneider.  Plasma flow in a Hall thruster.  Physics of Plasmas, 12(4):043502, 2005.
  * [49] H. Liu, B. Wu, D. Yu, Y. Cao, and P. Duan.  Particle-in-cell simulation of a Hall thruster.  Journal of Physics D: Applied Physics, 43(16):165202, apr 2010.
  * [50] T. Charoy.  Numerical study of electron transport in Hall thrusters.  PhD thesis, Institut polytechnique de Paris, 2020.
  * [51] P. Coche and L. Garrigues.  A two-dimensional (azimuthal-axial) particle-in-cell model of a Hall thruster.  Physics of Plasmas, 21(2):023503, 2014.