# On the mechanism of ionization oscillations in Hall thrusters

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  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 ionization oscillations involving plasma and neutral density
(breathing modes) are the most violent perturbations in Hall thrusters for
electric propulsion. Because of its simplicity, the zero-dimensional (0-D)
predator-prey model of two nonlinearly coupled ordinary differential equations
for plasma and neutral density has been often used for the characterization of
such oscillations and scaling estimates. We investigate the properties of its
continuum analog, the one-dimensional (1-D) system of two nonlinearly coupled
equations in partial derivatives (PDE) for plasma and neutral density. This is
a more general model, of which the standard 0-D predator-prey model is a
special limit case. We show that the 1-D model is stable and does not show any
oscillations for the boundary conditions relevant to Hall thruster and the
uniform ion velocity. We then propose a reduced 1-D model based on two coupled
PDE for plasma and neutral densities that is unstable and exhibit oscillations
if the ion velocity profile with the near the anode back-flow (toward the
anode) region is used. Comparisons of the reduced model with the predictions
of the full model that takes into account the self-consistent plasma response
show that the main properties of the breathing mode are well captured. In
particular, it is shown that the frequency of the breathing mode oscillations
is weakly dependent on the final ion velocity but shows a strong correlation
with the width of the ion back-flow region.

##  I Introduction

Hall thrusters are electric propulsion plasma devices with cross-field
$E\times B$ configuration, where electrons are effectively trapped by the
magnetic field, and heavy ions are accelerated by the electric field
generating thrust. Despite the long history and a large number of successful
missions, as well as the relatively simple design of these devices, many
physical phenomena are not well understood. Large amplitude discharge current
oscillations in the axial direction of Hall thrusters, or the so-called
breathing modes, are among the most strong perturbations that affect the
operation, e.g. with current oscillations reaching 100% and even extinguishing
the dischargeZhurin  _et al._ (1993). However note that some modifications of
Hall thrusters exhibit no or very little low-frequency oscillationsRaitses and
Fisch (2001); Wei  _et al._ (2012a). The physics of the breathing mode has
been studied in many papers but the full understanding is far from complete.
Simulations of low-frequency axial dynamics of Hall thrusters were done
earlier with various approximations, including fluidMorozov and Savelyev
(2000); Hagelaar  _et al._ (2002); Barral and Ahedo (2009), and hybridMorozov
and Savelyev (2000); Boeuf and Garrigues (1998); Fife (1998); Chable and
Rogier (2005); Hara, Boyd, and Kolobov (2012) (fluid electrons and kinetic
ions and atoms) models. Ionization is one of the key physics elements of the
breathing mode oscillations. In addition to the basic coupling of plasma
density and neutral gas due to the ionization, many models include self-
consistent electric field calculations along the channel, and some models
include electron pressure evolution and detailed electron energy balance.

At the same time, a simple zero-dimensional (0-D) model was proposed Fife  _et
al._ (2012); Barral and Ahedo (2008) to explain the mechanism of breathing
mode, in which neutral and ion (plasma) densities are coupled in the form of
predator-prey model (Lotka-Volterra equations). The predator-prey model is a
system of two nonlinearly coupled ordinary differential equations that
describes periodic oscillations in a number of physics and biology systems.
This model is claimed to predict the frequencies close to the observed in
experiments. It was, however, noted that for the constant ion and neutral
velocities, the model is stable and does not predict any instability Hara  _et
al._ (2014) so the conditions for the excitation of the ionization
oscillations are not clear. Additional effects and modifications Hara  _et
al._ (2014); Wang, Wei, and Yu (2011) were proposed for the predator-prey
model to make it unstable but it remains unclear to what degree the 0-D models
are capable to predict the breathing mode oscillations (discussed in section
II).

Nevertheless, simple reduced models are of considerable interest. First, such
models are useful as scaling tools, e.g. for the design of power supply
sources. Another motivation stems from the considerable complexity of the full
models. In the full models with many parameters that are unknown
experimentally, it is rather difficult to predict the conditions for the
instability and therefore difficult to translate the predictions of the
theoretical models into practical recommendations. The search for physical
mechanism (or mechanisms) and critical controlling parameters continues Dale
and Jorns (2019a). Reduced models focusing on selected specific phenomena may
allow easier experimental validation and eventually lead to a better
understanding of the physics of the instability, its conditions, and the
development of the full predictive model(s) Dale and Jorns (2019b, c).

In this work, we analyze the properties of the continuum one-dimensional (1-D)
model, based on two coupled partial differential equations (PDE) for the
evolution of neutral and plasma densities. The 0-D predator-prey is a
reduction of the 1-D model. If anything, the more complete 1-D model should
exhibit the same physics (and e.g. predict similar characteristics of the
oscillations) if the reduction from the 1-D to the 0-D model is appropriate
for the relevant physics, and the oscillations in the 0-D model are not the
artefact introduced by the simplifications. We show, however, that the 1-D
model is not only stable with boundary conditions of a typical Hall thruster
and the standard assumptions of the uniform velocity (as it is used in 0-D
model) but also does not show the oscillations observed in the 0-D model. On
the contrary, any oscillations, even introduced externally, are damped in time
and space. We then propose a modified 1D model with a non-uniform ion velocity
profile and with a back-flow (toward the anode) the region near the anode. We
show that such a model is unstable and exhibit ionization oscillations with
characteristics roughly consistent with the predictions of the full fluid
modelSmolyakov  _et al._ . Rather good agreement of the results of the reduced
(simplified) and full models suggests that the reduced model captures well
essential physics relevant to the breathing mode and may point to the critical
conditions for the instability.

The paper is organized as follows. In section II, various modifications of the
0-D model are discussed. In section III, the continuum 1-D model is
introduced. The stationary solutions are obtained here, their stability and
response to external perturbations are discussed for the uniform profile of
the ion velocity. Section IV presents a reduced 1-D continuum model with the
ion back-flow region near the anode and shows that this model exhibit unstable
nonlinear oscillations. In section V, the predictions of the reduced model are
compared against the full self-consistent simulations.

##  II Zero-dimensional (0-D) predator-prey models

In this section we discuss the various modifications and properties of the
zero-dimensional (0-D) predator-prey models for the ionization oscillations.
Such a simplest model includes basic time and space dynamics of plasma and
neutrals and is based on the continuity equations for the ion density $n_{i}$
(assuming full plasma quasineutrality), and neutral atom density $n_{a}$:

$\displaystyle\frac{\partial n_{a}}{\partial t}+v_{a}\frac{\partial
n_{a}}{\partial x}=-\beta n_{a}n_{i},$ (1) $\displaystyle\frac{\partial
n_{i}}{\partial t}+\frac{\partial}{\partial x}\left(n_{i}v_{i}\right)=\beta
n_{a}n_{i},$ (2)

where $v_{i},v_{a}$ are ion and atom flow velocities, respectively, and
$\beta$ is the ionization rate coefficient. We assume that $v_{a}$ is a
constant, but $v_{i}$ and $\beta$, in general, could be non-uniform in space.
A basic assumption leading to a simple zero-dimensional predator-prey model
Fife  _et al._ (2012) is that the ionization zone and acceleration zone is
replaced by the transition layer of the width $L$. The integration of Eqs.
(1)-(2) over this region results in the system of ordinary differential
equations (ODE):

$\displaystyle\frac{\partial}{\partial t}\left\langle
n_{a}\right\rangle+\frac{1}{L}\left.\left(n_{a}v_{a}\right)\right|_{0}^{L}=-\beta\left\langle
n_{a}n_{i}\right\rangle,$ (3) $\displaystyle\frac{\partial}{\partial
t}\left\langle
n_{i}\right\rangle+\frac{1}{L}\left.\left(n_{i}v_{i}\right)\right|_{0}^{L}=\beta\left\langle
n_{a}n_{i}\right\rangle,$ (4)

where averages over the transition layer are defined as
$\left\langle\left(...\right)\right\rangle=L^{-1}\int_{0}^{L}\left(...\right)dx$
and $\beta$ is assumed constant here. A set of approximations are made to
obtain the original predator-prey modelFife  _et al._ (2012). First, the ion
and neutral fluxes at the boundaries of the transition layer are assumed in
the form: $\left.\left(n_{i}v_{i}\right)\right|_{0}=0$,
$\left.\left(n_{i}v_{i}\right)\right|_{L}=nv_{i}$ and
$\left.\left(n_{a}v_{a}\right)\right|_{0}=n_{a}v_{a}$ and
$\left.\left(n_{a}v_{a}\right)\right|_{L}=0$, where $v_{i}$ is the final ion
velocity and $n(t)$ is plasma density at the exit of the transition layer, at
$x=L$, $n(t)\equiv n_{i}(t,L)$, and $N(t)$ is neutral density at $x=0$,
$N(t)\equiv n_{a}(t,0)$. Thus, these boundary conditions imply full
ionization, $n_{a}(L)=0$, zero flux of the ions from the left
boundary,$\left.\left(n_{i}v_{i}\right)\right|_{0}=0$, and full acceleration
the transition layer, $\left.\left(n_{i}v_{i}\right)\right|_{L}=nv_{i}$. Then,
the boundary values for the plasma and neutral densities are used to
approximate the averages over the transition layer: $\left\langle
n_{a}n_{i}\right\rangle\simeq nN$, $\left\langle n_{a}\right\rangle\simeq N$,
and $\left\langle n_{i}\right\rangle\simeq n$. These steps lead to the 0-D ODE
system Fife  _et al._ (2012)

$\displaystyle\frac{dN}{dt}-\frac{1}{L}Nv_{a}=-\beta Nn,$ (5)
$\displaystyle\frac{dn}{dt}+\frac{1}{L}nv_{i}=\beta Nn.$ (6)

These equations have the equilibrium solution with $n_{eq}=v_{a}/L\beta$ and
$N_{eq}=v_{i}/L\beta$. Considering the perturbations
$\left(\widetilde{n}(t),\widetilde{N}(t)\right)\sim\exp\left(-i\omega
t\right)$ near the equilibrium, $n=$ $n_{eq}+\widetilde{n}(t)$, $N=$
$N_{eq}+\widetilde{N}(t)$ one obtains stable oscillations with

$\omega^{2}=\beta^{2}n_{eq}N_{eq}=v_{i}v_{a}/L^{2}$ (7)

It was later noted Barral and Ahedo (2008) that oscillations (time dependence)
in $N$ are inconsistent with the constant mass rate boundary condition at
$x=0$, $\left.\left(n_{a}v_{a}\right)\right|_{0}=Nv_{a}=\dot{M}/Am_{a}=const,$
which means that the $N$ should be constant for the constant neutral flow
velocity; $A$ is the cross-section area. A modification of the 0-D model was
suggested in Ref. Hara  _et al._ , 2014 where the value of neutral density at
the left (entrance) side was fixed constant,
$\left.\left(n_{a}v_{a}\right)\right|_{0}=N_{0}v_{a}$, and the value at the
exit was assumed time dependent,
$\left.\left(n_{a}v_{a}\right)\right|_{L}=N\left(t\right)v_{a},$ so the
neutral balance equation takes the form

$\frac{dN}{dt}+\frac{1}{L}\left(N-N_{0}\right)v_{a}=-\beta Nn.$ (8)

The plasma balance equation was also modified with an additional term due to
the sheath losses at the lateral walls

$\frac{dn}{dt}+\frac{1}{L}nv_{i}+\frac{1}{d}nc_{s}=\beta Nn,$ (9)

where sheath losses are estimated based on the Bohm condition for the ion
velocity, where $d$ is the radial channel width and
$c_{s}=\sqrt{T_{e}/m_{i}}$. With these modifications, the equilibrium value of
the neutral density in equations (8-9) is only corrected by the sheath losses,

$N_{eq}=\frac{1}{\beta}\left(\frac{v_{i}}{L}+\frac{c_{s}}{d}\right),$ (10)

while the equilibrium value of plasma density can be very different from that
in equation (5) and takes the form

$n_{eq}=\frac{v_{a}}{\beta L}\left(\frac{N_{0}}{N_{eq}}-1\right).$ (11)

Note that in this model, the value of the neutral density at $x=0$ should be
sufficiently large $N_{0}>N_{eq}=\left(v_{i}/L+c_{s}/d\right)/\beta$.
Considering perturbations near $N_{eq}$ and $n_{eq}$ one obtains damped
oscillations Hara  _et al._ (2014)

$\omega=-\frac{i}{2}\frac{v_{a}}{L}\frac{N_{0}}{N_{eq}}\pm\sqrt{\beta^{2}n_{eq}N_{eq}-\frac{1}{4}\left(\frac{v_{a}}{L}\right)^{2}\left(\frac{N_{0}}{N_{eq}}\right)^{2}}.$
(12)

The essential problem of the predator-prey model is that it does not predict
instability so no condition for the oscillations can be determined. It was
argued that the electron dynamics should be important and more complex models
were proposed such as using the electron (instead of ion) continuity equation
with the drift-diffusion approximation for the electron velocity Yamamoto  _et
al._ (2002); Yamamoto, Komurasaki, and Arakawa (2005); Yamamoto  _et al._
(2005), ion flows Wang, Wei, and Yu (2011), electron energy evolution Hara
_et al._ (2014); Hara and Kubota (2017), and two-zone model Dale and Jorns
(2019b, c). Ref. Wei  _et al._ , 2015 provides an overview of various models
and involved mechanisms.

##  III Continuum (1-D) predator-prey model

Ideally, the reduction from the continuous 1-D model to the 0-D models as was
discussed in the previous section should preserve the essential properties and
features of the more complete model and do not introduce any fundamental
changes to the 1-D model. Therefore, is of interest to study if the 1-D model
has the properties predicted by the 0-D models. Here we discuss properties of
the stationary and time-dependent solutions of the one-dimensional
(continuous) model consisting of the continuity equations (1) and (2) for
atoms and ions with constant flow velocities $v_{a}$ and $v_{i}$, and the
ionization rate $\beta$. These are the same assumptions with which the
original predator-prey model was derived. The stationary problem is formulated
as two coupled equations:

$\displaystyle v_{a}\frac{\partial n_{a}}{\partial x}=-\beta n_{a}n_{i},$ (13)
$\displaystyle v_{i}\frac{\partial n_{i}}{\partial x}=\beta n_{a}n_{i}.$ (14)

The exact solution to this system is found in the form:

$\displaystyle
n_{a,st}(\xi)=n_{a0}\left(1+\frac{n_{0}v_{i}}{n_{a0}v_{a}}\right)\dfrac{1}{\left(n_{0}v_{i}\right)/\left(n_{a0}v_{a}\right)\exp{\left(\xi\right)}+1},$
(15) $\displaystyle
n_{i,st}(\xi)=n_{0}\left(1+\frac{n_{a0}v_{a}}{n_{0}v_{i}}\right)\dfrac{\exp{\left(\xi\right)}}{\exp{\left(\xi\right)}+\left(n_{a0}v_{a}\right)/\left(n_{0}v_{i}\right)},$
(16)

where $n_{0},n_{a0}$ are values of ion density and atom density, respectively,
at the left boundary, and normalized length is $\xi=x/l_{0}$, with
$l_{0}=v_{i}v_{a}/\beta\left(n_{0}v_{i}+n_{a0}v_{a}\right)$. It can be noted
that the solutions of the stationary problem depend only on two parameters,
$n_{0},n_{a0}$. Typical stationary solutions for $v_{i}/v_{a}=10$ are depicted
in Figs. 1a,1b by solid lines. Note that the position of the crossing point,
where the ion and neutral densities are equal, and the localization of the
ionization source $\beta nn_{a}$ depend on the value $n_{0}$: they move to the
right with decreasing $n_{0}$, and goes to infinity for $n_{0}\to 0$.

(a) $n_{0}/n_{a0}=0.05,\ l_{0}=$1.31\text{\,}\mathrm{c}\mathrm{m}$$

(b) $n_{0}/n_{a0}=0.0001,\ l_{0}=$1.96\text{\,}\mathrm{c}\mathrm{m}$$

Figure 1: Stationary solutions given by Eqs. (13) and (14) are shown by solid
lines, with higher (a) and lower (b) ion density at the left boundary $n_{0}$.
The time dependent solutions of Eqs. (1), and (2), converging to the
stationary solution are shown for the case (a) at three consecutive time steps
$t_{1}<t_{2}<t_{3}$, dot dashed (initial condition), dashed, and dotted line,
respectively.

Numerical studies of the time-dependent equations (1) and (2) were performed
with fixed boundary values at the left $n_{0},n_{a0}$ and free boundary
conditions (spatial second derivative is zero) on the right end. This analysis
shows that all perturbations are damped and converge to the stationary
solutions given by equations (15) and (16), as shown in Fig. 1a at different
times. Thus, the 1-D predator-prey model with a constant spatial profile of
the ion velocity does not exhibit oscillatory or unstable behaviour, unlike
zero-dimensional predator-prey models. The 1-D continuum model presented here
in some sense is similar to the asymptotic low-frequency model derived in the
limit $n/N_{a}\ll 1$ Barral and Ahedo (2008, 2009), which also predicts the
neutral oscillations with the frequencies of the order of given by equation
(7). Contrary to Barral and Ahedo (2008, 2009), our model does not show
neutrally stable oscillations but only damped modes.

One can envisage that the boundary value of plasma density could be perturbed,
e.g. by some sheath instability. Such perturbations may also arise in the
experiments with segmented anode Raitses  _et al._ (2006), external
modulations of the anode voltage Romadanov, Raitses, and Smolyakov (2018), and
two-stage thruster configurations Fisch  _et al._ (2001). Therefore, it is
interesting to investigate how this nonlinear system responds to external
perturbations of the ion density and whether such external perturbations may
grow in amplitude while propagating from the left boundary. We will impose
harmonic external modulations of the ion density at the left boundary in the
form $n_{i}(0)=n_{i0}\left(1+r\sin{2\pi f_{ext}t}\right)$, where $f_{ext}$ is
the modulation frequency, with $r=0.15$. We consider the 1-D continuum model,
Eqs. (1,2) with $v_{i}=10v_{a}$, $n_{i0}/n_{a0}=0.05$,
$l_{0}=$1.31\text{\,}\mathrm{c}\mathrm{m}$$, Fig. 1a, with the simulation
length $15l_{0}$. Since velocities of neutrals and ions are different, there
are two natural modes that can be excited during the external modulations,
with the wavelengths given with $\lambda_{a}=v_{a}/f_{ext}$ and
$\lambda_{i}=v_{i}/f_{ext}$. We consider four values of the driving
frequencies $f_{ext}l_{0}/v_{a}=(0.1,0.5,1,10)$. For these values, the
corresponding natural wavelengths of the neutral and ion characteristic
wavelengths are $\lambda_{a}=(10,2,1,0.1)l_{0}$ and
$\lambda_{i}=(100,20,10,1)l_{0}$, respectively.

(a) $f_{ext}l_{0}/v_{a}=0.1$

(b) $f_{ext}l_{0}/v_{a}=0.5$

(c) $f_{ext}l_{0}/v_{a}=1$

(d) $f_{ext}l_{0}/v_{a}=10$

Figure 2: Instantaneous profiles of perturbed ion $\widetilde{n}_{i}$ (blue)
and neutral $\widetilde{n}_{a}$ (red) densities for the various driving
frequencies $f_{ext}l_{0}/v_{a}=(0.1,0.5,1,10)$, shown, respectively, in (a),
(b), (c), (d). Three consecutive time snapshots $t_{1}<t_{2}<t_{3}$ are shown
for each variable, depicted with the dot dashed line, dashed line, and solid
line, respectively; they divide one oscillation period in 3 equal intervals.
Crossing point (of stationary ion and atom densities profiles) location is
shown with the black dashed line.

The resulting response of the system of Eqs. (1,2) under the external
modulations of the ion density is presented in Fig. 2 for the perturbed
variables, defined as $\widetilde{n}=n-n_{st}$, where $n_{st}$ is the
corresponding stationary solution. It can be noted that in all cases ion
density perturbations grow in absolute values before reaching the crossing
point and leaving the ionization region, whereas the amplitude of neutral
density perturbations decays, which is expected in the ionization region. Note
that while the absolute value of the perturbations of the ion density seems
increasing, the relative value $\widetilde{n}/n_{st}$ does not grow. Thus, in
the presence of external ion density perturbations in the near-anode region,
the system exhibit externally driven oscillations that are advected along the
channel but show a limited increase of the amplitude. However, the atom
response is slightly more complicated and besides the natural mode can exhibit
the nonlinear response (on ion natural wavelength). It is also affected by the
atom stationary profile (Fig. 1a), as atom perturbations decay rapidly after
the crossing point. It is seen that for the low frequency perturbations the
natural atom mode is excited (Fig. 2a), while for the higher frequencies the
ion characteristic wavelengths become more and more dominant in the atom
response (seen gradually through Figs. 2b-2d). We remind that the amplitude of
the ion density perturbations was fixed in these cases. The amplitude of the
neutral density perturbations is maximal for the lower frequencies. We have
also investigated the case when the neutral density was modulated at the left
boundary. Such a system shows similar behaviour and does not show the unstable
modes.

##  IV Ionization oscillations in the 1-D continuous model with the ion back-
flow in the near anode region

As it was discussed in the previous section, the reduced fluid model,
described with Eqs. (1, 2), with spatially uniform profiles of $v_{i}$,
$v_{a}$, and $\beta$ is stable. Here we will show that the same system with
spatial dependence of $v_{i}$, when $v_{i}$ is negative near the anode, the
so-called ion back-flow, becomes unstable and exhibit self-consistent
oscillations. The ion back-flow region naturally forms in the presheath, a
quasineutral region near the anode with the negative electric field due to the
electron diffusionCohen-Zur  _et al._ (2002); Ahedo and Rus (2005); Dorf,
Semenov, and Raitses (2003); Dorf, Raitses, and Fisch (2005). The positive
density gradient near the anode creates a strong electron flow to the anode.
To maintain ambipolarity (the total current remains uniform), the ion back-
flow current occurs to compensate the increase of the electron current. The
current ambipolarity is required to keep plasma quasineutral. This process is
somewhat similar to the formation of the electric field in the presheath
region near the plasma boundary where the electric field is induced to
accelerate ions toward the boundary (up to the Bohm velocity) to balance
electron and ion losses.

In this configuration the only fixed boundary condition is a value of the
neutral density on the left, as their characteristics travel strictly to the
right (atom velocity is a positive constant). Ion density boundary condition
is free on the left, $\partial^{2}_{x}n_{i}(0)=0$. First, we consider the
simple configuration with the velocity profile $v_{i}$ as a strictly linear
function of the position. The ionization rate $\beta$ is taken uniform. One
might argue the constant $\beta$ is an unrealistic assumption for a typical
Hall thruster where electron energy near the anode is low due to higher
mobility, shaping the $\beta$ profile; we consider these effects in Section V
in the comparison with the full fluid model. Recombination of plasma at the
anode is not included here: the ion flux converted to the neutral flux at the
anode is not added to the neutral flow. It can also be included, increasing
the amplitude of oscillations, but it is not required for the existence of the
oscillations. An example of such oscillations shown for a typical Hall
thruster parameters (the stationary plasma thruster “SPT 100”Morozov and
Savelyev (2000)), with the channel length $3\text{\,}\mathrm{c}\mathrm{m}$,
ionization rate
$\beta=$9\text{\times}{10}^{-14}\text{\,}\mathrm{m}^{3}\mathrm{s}^{-1}$$,
$v_{a}=$150\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$. Atom density at the
anode $n_{a0}$ is evaluated from $n_{a0}v_{a}=\dot{m}/Am_{a}$, assuming
$\dot{m}=$5\text{\,}\mathrm{m}\mathrm{g}\mathrm{s}^{-1}$$, Xenon atom mass,
and $A=$12.75\pi\text{\,}\mathrm{c}\mathrm{m}^{2}$$. For the simplicity, we
assumed the linear ion flow velocity profile
$v_{i}=\left(2.5x-1.5\right)\leavevmode\nobreak\
$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$ ($x$ is in
[$\mathrm{c}\mathrm{m}$]), with the ion back-flow extent
$L_{b}=$0.6\text{\,}\mathrm{c}\mathrm{m}$$. The resulting time and space
evolution of ion and atom densities is shown in Fig. 3. The observed frequency
is $12.8\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$, in the range of the values
observed in experiments and simulations for breathing mode.

(a)

(b)

Figure 3: Evolution in space and time of atom (a) and ion (b) density,
normalized to $n_{a0}$.

Now, for the same setup with ion velocity profile in the form
$v_{i}=\left(2.5x-1.5\right)\leavevmode\nobreak\
$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$ ($x$ is in
[$\mathrm{c}\mathrm{m}$]), we will show effects of the amplitude of the
ionization rate $\beta$, which is still assumed uniform here for simplicity.
The ionization rate coefficient $\beta$ needs to be sufficiently large to
support plasma discharge. Interestingly, in addition to a threshold $\beta$
value for the discharge to exist, there is also a threshold between stable and
oscillatory plasma behaviour. In the oscillatory regime, there is a value of
$\beta$ when the amplitude of the oscillations is the largest (note that here
we assume that $\beta$ is constant along the channel). It shouldn’t be of
surprise, as when the values of $\beta$ are too high, the ionization takes
place immediately near the left end and so both density profiles have maximum
values at the left end. The dependency of the oscillation amplitude on $\beta$
is shown in Fig. 4a, where the peak-to-peak amplitude of the averaged in space
density $\langle n_{i}\rangle_{x}(t)$ is shown.

The frequency dependence on $\beta$ is presented in Fig. 4b. If the ion back-
flow region is very short, i.e. $v_{i}(0)\approx 0$, the ionization mostly
will occur at the anode for the same reasons, and there would not be enough
spatial separation between high and low atom density points. Dependencies,
similar to those shown in Figs. 4a, 4b, are also observed for different values
of the ion back-low extent $L_{b}$. They follow the same scaling as presented
below in Figs. 5b, 6.

(a)

(b)

Figure 4: Peak-to-peak values of averaged in space ion density $\langle
n_{i}\rangle_{x}(t)$ as a function of ionization rate (a). The main frequency
component of observed oscillations as a function of the ionization rate (b).

The effect of the ion back-flow region is investigated further for our reduced
model by taking the ion velocity profiles with a variable ion back-flow
extent. The Fig. 5a shows the velocity profiles taken as parabolas specified
with the three points,
$v_{i}(0)=$-1.5\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$,
$v_{i}(L)=$20\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$, and
$v_{i}(L_{b})=$0\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$, where
$L_{b}$ is the length of the ion back-flow region, the distance from the anode
to the transition point where the ion velocity reverse its sign. The value of
$L_{b}$ was varied, with the system length set to
$3\text{\,}\mathrm{c}\mathrm{m}$, the ionization rate
$\beta=$5\text{\times}{10}^{-14}\text{\,}\mathrm{m}^{3}\mathrm{/}\mathrm{s}$$,
atom flow velocity $v_{a}=$150\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$.
Simulations using these velocity profiles reveal that that oscillation
frequency scales with the atom fly-by frequency $v_{a}/L_{b}$ in the back-flow
region, Fig. 5b. The oscillations amplitude decreases with decreasing $L_{b}$
and the oscillations disappear for $L_{b}\to 0$, as shown in Fig. 6. Also, no
oscillations were observed for monotonically increasing ion velocity profiles
with $v_{i}(0)\geq 0$.

(a)

(b)

Figure 5: Velocity profiles (a) for a variable ion back-flow region $L_{b}$
while the exhaust velocity is fixed to
$20\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$
($300\text{\,}\mathrm{e}\mathrm{V}$ is a typical ion energy in the exhaust
plume of “SPT100”Morozov and Savelyev (2000)). The corresponding frequency (b)
as a function of $L_{b}$ in compare with neutral flyby frequency of
$v_{a}/L_{b}$ in the backflow region. (a)  Figure 6: Peak-to-peak values of
averaged in space ion density $\langle n_{i}\rangle_{x}(t)$ as a function of
the ion back-flow region $L_{b}$. Oscillations disappear as
$L_{b}=$1\text{\,}\mathrm{m}\mathrm{m}$$ goes to zero.

Zero-dimensional predator-prey models predict that the frequency depends on
the exhaust ion velocity. Here, we show that for our model this effect is very
small and, mostly, the back-flow region “defines” the frequency. Fig. 7a shows
the family of the ion velocity profiles, constructed in a similar manner (by a
parabola) but keeping the back-flow length the same
($1\text{\,}\mathrm{c}\mathrm{m}$) and varying
$v_{i}(L)=(14,18,22,26)\leavevmode\nobreak\
$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$. The resulting frequency for all
these profiles was found to be close to 10 kHz, Fig. 7b, predicted by the
previous example in Fig. 5b with the $L_{b}=$1\text{\,}\mathrm{c}\mathrm{m}$$.

(a)

(b)

Figure 7: Velocity profiles (a) and corresponding frequency (b) as a function
of $v_{i}(L)$, for the variable ion exhaust velocity $v_{i}(L)$ while keeping
the back-flow region $L_{b}$ constant.

The mechanism of the oscillations can be considered in three main stages,
starting from a stationary state in Fig. 1a, but with the ion velocity back-
flow region. In the first stage we observe enhanced ionization process in the
region of plasma stagnation ($v_{i}\approx 0$). The ionization process in this
stage may be approximated as $\dot{n}_{i}=\beta n_{i}n_{a}$,
$\dot{n}_{a}=-\beta n_{i}n_{a}$ (where $n_{i},n_{a}$ are some average plasma
and atom densities, respectively, in that region). In this stage, effects of
atom inflow and ion outflow are relatively small here. The ionization proceeds
with the exponential rise plasma density $n_{i}\sim\exp{\left(\beta
n_{a}t\right)}$ and depletion of atom density. In the second stage, the newly
created plasma peak splits and partially moves to the anode (via back-flow)
wiping down the rest of atoms near the anode. Finally, in the third stage, the
system enters the period of “regeneration” with both plasma and atom densities
low, while atoms slowly fill the region near the anode. The atom density
remains higher than the ion density, and the ion density balance is
predominantly follows the equation $\dot{n}_{i}=\beta
n_{i}n_{a}-\gamma_{loss}n_{i}$, where $\gamma_{loss}n_{i}$ represents the ion
advection losses via the back-flow. In this quasistationary state, there is no
exponential rise of ionization but slow inflow of atoms (to the right) filling
this region toward the stagnation region. The third stage is also the slowest
and defines the breathing mode period, hence the backflow lengths defines the
period (Fig. 5b). Additionally, if the recombination at the anode is included,
it increases a number of atoms injected to the system for the stage three, and
it increases the amplitude of the oscillations.

##  V Reduced vs full self-consistent model

In this section we compare the reduced and full models for the breathing
modes, and demonstrate that our simple reduced model, Eqs. (1,2), reproduces
well the results from the full self-consistent time-dependent model Smolyakov
_et al._ . First, we give a brief description of the full fluid self-
consistent model for low-frequency ionization oscillations.

###  V.1 Full self-consistent model

The full model of low-frequency axial plasma dynamics in Hall thruster is
considered in electrostatic and quasineutral approximation for three species:
neutrals, ions, and electrons. Our model is fully fluid both for the ion,
neutral, and electron components. The electron equations in general are
similar to those used in Refs. Boeuf and Garrigues (1998); Barral and Ahedo
(2009). The full fluid model consists of time-dependent PDE equations for
neutral atom density $n_{a}$, plasma density $n$ (ion and electron), ion flow
velocity $v_{i}$, and electron temperature $T_{e}$:

$\displaystyle\frac{\partial n_{a}}{\partial t}+v_{a}\frac{\partial
n_{a}}{\partial x}=-\beta n_{a}n,$ (17) $\displaystyle\frac{\partial
n}{\partial t}+\frac{\partial}{\partial x}\left(nv_{i}\right)=\beta n_{a}n,$
(18) $\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),$ (19)
$\displaystyle\frac{3}{2}\frac{\partial}{\partial
t}\left(nT_{e}\right)+\frac{5}{2}\frac{\partial}{\partial
x}\left(nv_{ex}T_{e}\right)+\frac{\partial q_{e}}{\partial
x}=nv_{ex}E-nn_{a}\mathrm{K}-n\mathrm{W},$ (20)

where the electric field $E$ is obtained from the electron momentum balance
equation

$v_{ex}=-\mu_{e}E-\frac{\mu_{e}}{n_{e}}\frac{\partial(nT_{e})}{\partial x},$
(21)

obtained with neglected electron inertia, where $\mu_{e}$ is the electron
mobility across the magnetic field (described below). Quasineutrality
condition, along with the constraint $\int_{L}Edx=U_{0}$ ($U_{0}$ is the
applied potential) leads to the total current
$J_{T}=en\left(v_{i}-v_{e}\right)$

$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}}},$ (22)

which is spatially uniform but may oscillate in time. Other quantities in the
system (17-20) are the constant atom flow velocity $v_{a}$, the ionization
rate coefficient $\beta$ (obtained with BOLSIG+Hagelaar and Pitchford (2005)
for Maxwellian EEDF using SIGLO databasesig ), the elementary charge $e$, the
electron mass $m_{e}$, the ion (Xenon) mass
$m_{i}=$131.293\text{\,}\mathrm{a}\mathrm{m}\mathrm{u}$$, the anomalous energy
loss coefficientBoeuf and Garrigues (1998) $\mathrm{W}$, the collisional
energy loss coefficient $\mathrm{K}$ (also via BOLSIG+), and the electron heat
flux given by

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

The electron transport across magnetic field is described in the form of the
magnetized mobility

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

where $\nu_{m}$ the total electron momentum exchange collision frequency,
$\omega_{ce}=eB/m_{e}$ is the electron cyclotron frequency. In this model
$\nu_{m}$ is represented in the form as was adopted in Ref. Hagelaar  _et al._
(2004),

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

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},$ (26)
$\displaystyle\nu_{walls}=\alpha${10}^{7}\text{\,}\mathrm{[}\mathrm{s}^{-1}\mathrm{]}$,$
(27) $\displaystyle\nu_{B}=\left(\beta_{a}/16\right)eB/m_{e}.$ (28)

where
$k_{m}=$2.5\text{\times}{10}^{13}\text{\,}\mathrm{m}^{-3}\mathrm{s}^{-1}$$,
$\alpha$ and $\beta_{a}$ are adjusting constants. The profile of external
magnetic field $B$ is shown if Fig. 8a, with the channel’s exit in the peak of
magnetic field intensity, given by
$B=B_{0}\exp{\left(-(x-l)^{2}/2\delta_{B}^{2}\right)}$, where
$l=$2.5\text{\,}\mathrm{c}\mathrm{m}$$ is the channel length and $\delta_{B}$
defines a magnetic field profile width. This model use different parameters
inside ($x<l$) and outside ($x\geq l$) the channel Hagelaar  _et al._ (2003,
2004); Hagelaar (2007), the near wall conductivity contribution (27)
$\alpha_{in}=0.2,\ \alpha_{out}=0$, and the anomalous contribution (28) is set
to $\beta_{a,in}=0.1$, $\beta_{a,out}=1$. The anomalous electron energy loss
coefficientBoeuf and Garrigues (1998) $\mathrm{W}$ is modeled as

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

where $\varepsilon=3/2T_{e}$, $U=$20\text{\,}\mathrm{V}$$, and
$\nu_{\varepsilon}$ is electron energy anomalous loss coefficient. A constant
mass flow rate $\dot{m}$ determines the value of $n_{a}$ at the boundary
together with the recombination of plasma that flows to the anode, hence the
boundary condition:

$n_{a}(0)=\frac{\dot{m}}{m_{i}Av_{a}}-\frac{nv_{i}(0)}{v_{a}},$ (30)

where $A$ is the anode surface area of a thruster. Bohm type condition for ion
velocity can be imposed at the anode $v_{i}(0)=-b_{v}\sqrt{T_{e}/m_{i}}$,
where $b_{v}=0\text{--}1$ is the Bohm velocity factor which can be varied.
Both anode and cathode electron temperature are fixed with
$T_{e}(0)=T_{e}(L)=2\leavevmode\nobreak\ $\mathrm{e}\mathrm{V}$$. All other
boundary conditions are free (spatial second derivative is zero).

Following the LANDMARK benchmarklan (2018) Test Case 3, we use the parameters
that result in bulk low-frequency oscillations, with the electron energy
anomalous loss coefficient inside the channel
$\nu_{\varepsilon,\text{in}}=10^{7}\leavevmode\nobreak\ $\mathrm{s}^{-1}$$ and
outside $\nu_{\varepsilon,\text{in}}=0.4\cdot 10^{7}\leavevmode\nobreak\
$\mathrm{s}^{-1}$$. The resulting plasma currents evolution shown in Fig. 8b,
assuming the typical inner and outer radius of a Hall thruster (“SPT
100”Morozov and Savelyev (2000))
$R_{1}=$3.5\text{\,}\mathrm{c}\mathrm{m}$,R_{2}=$5\text{\,}\mathrm{c}\mathrm{m}$$,
respectively. The input atom mass flow used in the simulation is
$\dot{m}=$5\text{\,}\mathrm{m}\mathrm{g}\mathrm{s}^{-1}$$ corresponds to
$I=\dot{m}q_{e}/m_{i}=$3.67\text{\,}\mathrm{A}$$. The average ion current (to
the exit plane) is $3.68\text{\,}\mathrm{A}$ which is consistent with the mass
flow rate. The electron current (to the exit plane) is of the same order (or
slightly larger). Other simulation parameters for this test case are:
$\delta_{B,in}=$1.1\text{\,}\mathrm{c}\mathrm{m}$$,
$\delta_{B,out}=$1.8\text{\,}\mathrm{c}\mathrm{m}$$,
$A=\pi(R_{2}^{2}-R_{1}^{2})$, $b_{v}=1$.

(a)

(b)

Figure 8: The magnetic field profile (a) used in full fluid and hybrid
simulations, with the channel exit located $2.5\text{\,}\mathrm{c}\mathrm{m}$
from anode (dashed line). Low-frequency oscillations in time (b) of total
current $I$, ion current $I_{i}$ (at $x=$5\text{\,}\mathrm{c}\mathrm{m}$$),
and electron current $I_{e}$ (at $x=$5\text{\,}\mathrm{c}\mathrm{m}$$) in the
full fluid model using specified parameters.

The self-consistent fluid model used here is essentially based on the
formulations and parameters suggested in Refs. Boeuf and Garrigues (1998);
Hagelaar  _et al._ (2002, 2004); Hagelaar (2007), also see the LANDMARK
benchmarklan (2018). The fluid model has been compared against the hybrid
model Hagelaar  _et al._ (2004); hal ; both models demonstrate similar results
Smolyakov  _et al._ .

###  V.2 Comparison of the reduced model with predictions of the full self-
consistent model

We compare now results of the full fluid model (in the regime described
above), with the reduced model. In the full model the ion velocity profile
$v_{i}$ and the ionization rate $\beta$ profile are self-consistent variables,
and for the reduced model we take averaged in time profiles of these
variables, Figs. 9a, 9b. We also include the ion back-flow recombination to
the atom boundary condition at the anode in the reduced model, described by
Eq. (30), in agreement with the full model.

(a)

(b)

Figure 9: Ion velocity (a) and ionization rate (b) coefficient profiles
obtained from the full fluid model as time-averaged over a few oscillation
periods. They are used in the reduced model as fixed profiles in time.

Time evolution of $n_{a}$, $n_{i}$, and the source $S=\beta n_{a}n_{i}$, are
shown in Figs. 10a, 10d, 10g, respectively. The same quantities from the
reduced model are shown in Figs. 10b, 10e, 10h. It is seen that the reduced
model results in oscillatory behaviour, qualitatively reproducing the
breathing mode oscillations, but with much higher amplitudes (about one order)
of evolving variables. The observed natural frequency in this model is also
higher, 14.9 kHz, compared to 10.2 kHz in the full model.

The observed discrepancy can be explained by the absence of the electron
temperature evolution and lack of self-consistent ionization rate dynamics in
the reduced model which is known to affect the breathing mode characteristics,
e.g. Ref. Wei  _et al._ , 2012b. One observation is that the average position
of the source term $S$ in our reduced model is closer to the anode (Fig. 10h)
(compared with the full model) that results in higher frequency. Higher
amplitude may also be explained by the shift of the ionization source toward
the region with higher ion back-flow velocity so that more ions return to the
anode and recombine enhancing the positive feedback loop. We found that simply
lowering the values of the ionization rate, precisely by taking $0.82\beta$ of
the original averaged $\beta$ profile (shown in Fig. 9b) agrees with the full
model much closer, see Figs. 10c, 10f, 10i, with the oscillation frequency of
11.8 kHz. While our intention is not really to achieve full quantitative
agreement, this exercise shows the sensitivity to the temperature effects. It
is also possible that the strong nonlinear dependence of the ionization rate
on temperature may suggest some weighting average for the ionization rate. We
leave the investigation of the proper averaging technique for future work.

(a) $n_{a}$ (Full)

(b) $n_{a}$ (Reduced)

(c) $n_{a}$ (Reduced, lower $\beta$)

(d) $n_{i}$ (Full)

(e) $n_{i}$ (Reduced)

(f) $n_{i}$ (Reduced, lower $\beta$)

(g) Source (Full)

(h) Source (Reduced)

(i) Source (Reduced, lower $\beta$)

Figure 10: Neutral density (a,b,c), ion density (d,e,f), and ionization source
$\beta n_{i}n_{a}$ (g,h,i) for the full fluid model, the reduced model, and
the reduced model with lower values of $\beta$ (see captions). Density values
are normalized to ${10}^{19}\text{\,}\mathrm{m}^{-3}$, and source values to
${10}^{24}\text{\,}\mathrm{m}^{-3}\mathrm{s}^{-1}$. The spatial domain is
limited to channel region only, $2.5\text{\,}\mathrm{c}\mathrm{m}$ from the
anode.

To show that the atom velocity $v_{a}$ plays a defining role in the observed
oscillation frequency, the full and reduced model (lower $\beta$) were
compared against each other, Fig. 11. Both full and reduced (with lower
$\beta$) model show approximately linear behaviour for the frequency, again
suggesting the feedback loop mechanism supported by ion back-flow and
advection of the neutral atoms. Interestingly, the ion current amplitudes are
comparable between the full and reduced (lower $\beta$) models, Fig. 12.

(a)  Figure 11: Oscillation frequency as a function of neutral velocity in
full and reduced (with lower $\beta$) fluid models.

(a)  (b)

Figure 12: Temporal and spatial ion current dynamics in full model (a) and
reduced model with lower $\beta$ (b), in Amperes (assuming cross section area
$A=$12.75\pi\text{\,}\mathrm{c}\mathrm{m}^{2}$$).

##  VI Summary

In this paper, we have proposed a reduced continuum 1-D model with the ion
velocity profile that has an ion back-flow (toward the anode) region. The
model consists of two one-dimensional equations for coupled dynamics of the
ion and neutral densities. While the 0-D simplification of this model,
commonly referred to as the predator-prey model, shows the oscillations, the
continuum 1-D model with a uniform ion velocity and boundary conditions
typical of Hall thruster show damped modes converging to the analytical
stationary state. Under externally imposed perturbations at the left boundary,
the 1-D continuum model shows limited convective amplification of the
perturbation (in the region where the local value of the neutral density is
smaller than plasma density) and subsequent convection (with a slow amplitude
decay) where the ion and neutral density perturbations are mostly decoupled.
Therefore the models with the uniform ion velocity, in the absence of
additional effects e.g. such as temperature variations Hara  _et al._ (2014),
do not predict the excitation of the breathing modes.

The key element of the proposed continuum model is the ion velocity profile
that has an ion back-flow (toward the anode) region. As it is shown here, such
a system becomes unstable and shows self-consistent nonlinear oscillations.
The comparison between this reduced model and the full self-consistent fluid
model shows very similar characteristics for oscillations of the ion and
neutral densities, and ion current. We have investigated the scaling of the
frequency and amplitude of the oscillations as a function of the width of the
back-flow region and neutral flow velocity which show generic $v_{a}/L_{b}$
dependence both in the reduced and full models. These results suggest that the
presence of the ion back-flow region may be a critical condition for the
breathing mode oscillations and the overlap of this region with the ionization
region creates the closed feedback loop necessary for the instability. These
results also indicate that the electron dynamics, in particular diffusion,
which leads to the appearance of the back-flow region, is important for the
breathing mode. A number of full models, reporting the breathing mode
oscillations, include the electron diffusion either explicitly via the
electron fluid equations, such as full fluid Barral and Ahedo (2009) or
hybridHagelaar  _et al._ (2004) formulations, or via full kinetic treatment as
in the full kinetic (particle-in-cell) treatment Coche and Garrigues (2014).
The ion back-flow near the anode occurs as a result of the electron diffusion
ensuring the current conservation to maintain quasineutrality. The ion
stagnation in the back-flow region triggering enhanced ionization is suggested
as the mechanism for the instability. It is also worth noting that the
electron diffusion is also a reason for the appearance of the singular sonic
point where the ion velocity is equal to the local sound velocity,
$v_{i}=c_{s}$. The smooth (regular) solution is obtained when the sonic point
is made regular by imposing some constraints on the current, mass flow rate,
and plasma parameters. These constraints made singular point regular and play
an important role for the existence and characteristics of stationary
solutions Dorf, Semenov, and Raitses (2003); Cohen-Zur  _et al._ (2002);
Ahedo, Martinez-Cerezo, and Martinez-Sanchez (2001); Makowski  _et al._
(2001). As it was shown in Smolyakov  _et al._ ; Romadanov  _et al._ (2020)
the constraints indeed define the stability of the resulting profile and
characteristic frequency of the breathing mode.

One has to note, however, that low frequency ionization oscillations were also
observed in the models without the electron diffusion, and hence, without
back-flow region Boeuf and Garrigues (1998); Morozov and Savel’ev (2000,
2000); Hara and Hanquist (2018); Hara, Boyd, and Kolobov (2012). It has been
suggested Chable and Rogier (2005) that the high frequency oscillations due to
the resistive axial modes Chable and Rogier (2005); Fernandez  _et al._
(2008); Koshkarov  _et al._ (2018) is a driving mechanism for a breathing
mode. Alternative mechanisms may be based on the temperature dependence of the
ionization coefficient Hara  _et al._ (2014), other electron temperature
effects Staack, Raitses, and Fisch (2004), and/or more complex interactions
between the low frequency modes involving ionization and higher frequency
modes in the ion transit time frequency range as it was suggested earlier
Kurzyna  _et al._ (2008); Smolyakov  _et al._ . The simulations with the full
model also demonstrate the existence of two different regimes of the breathing
mode Smolyakov  _et al._ , the so-called solo regime, where only the low
frequency mode present, and the regime with coexisting low frequency and the
high frequency ion transit time oscillations. The breathing mode oscillations
reproduced by our reduced model are similar to the solo regime, but are
different from the modes coexisting with high-frequency ion transit time
oscillations. The existence and nature of the breathing mode together with the
well pronounced high-frequency mode is an interesting question that will be
addressed in future studies with the full model.

In summary, we would like to conclude that the importance of the reduced
models is not only in their ability to predict reasonably well some essential
features of the ionization modes in Hall thrusters but also, and perhaps, even
more importantly, in pointing to some important physics involved in the modes
excitation and characteristics. The reduced models typically involve a smaller
number of adjustable parameters compared to the full models, e.g. as the
values of the anomalous mobility are poorly known. Therefore, the tests of the
properties and consequences of the reduced models in the experiments
Romadanov, Raitses, and Smolyakov (2018); Dale and Jorns (2019b, c, a) could
be an effective approach to test the crucial physics responsible for the
breathing modes oscillations.

###### Acknowledgements.

This work is supported in part by US Air Force Office of Scientific Research
FA9550-15-1-0226, NSERC Canada, and Compute Canada computational resources.
The authors acknowledge with gratitude the important discussions with I. V.
Romadanov, J. B. Simmonds, I.D. Kaganovich, and K. Hara.

## Data availability

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

## References

  * Zhurin  _et al._ (1993) V. Zhurin, J. Kahn, H. Kaufman, K. Kozubsky, and M. Day, “Dynamic characteristics of closed drift thrusters,” in  _23rd International Electric Propulsion Conference_ (1993).
  * Raitses and Fisch (2001) Y. Raitses and N. J. Fisch, “Parametric investigations of a nonconventional Hall thruster,” Physics of Plasmas 8, 2579–2586 (2001).
  * Wei  _et al._ (2012a) L. Wei, K. Han, C. Wang, H. Li, C. Zhang, and D. Yu, “Study on breathing mode oscillation suppression of self-excited Hall thrusters,” Journal of Vacuum Science & Technology A 30, 061304 (2012a).
  * Morozov and Savelyev (2000) A. Morozov and V. Savelyev, “Fundamentals of stationary plasma thruster theory,” in  _Reviews of Plasma Physics_ (Springer, 2000) pp. 203–391.
  * Hagelaar  _et al._ (2002) G. Hagelaar, J. Bareilles, L. Garrigues, and J.-P. Boeuf, “Two-dimensional model of a stationary plasma thruster,” Journal of Applied Physics 91, 5592–5598 (2002).
  * Barral and Ahedo (2009) S. Barral and E. Ahedo, “Low-frequency model of breathing oscillations in Hall discharges,” Physical Review E 79, 046401 (2009).
  * Boeuf and Garrigues (1998) J. Boeuf and L. Garrigues, “Low frequency oscillations in a stationary plasma thruster,” Journal of Applied Physics 84, 3541–3554 (1998).
  * Fife (1998) J. M. Fife,  _Hybrid-PIC modeling and electrostatic probe survey of Hall thrusters_ , Ph.D. thesis, Massachusetts Institute of Technology (1998).
  * Chable and Rogier (2005) S. Chable and F. Rogier, “Numerical investigation and modeling of stationary plasma thruster low frequency oscillations,” Physics of Plasmas 12, 033504 (2005).
  * Hara, Boyd, and Kolobov (2012) K. Hara, I. Boyd, and V. Kolobov, “One-dimensional hybrid-direct kinetic simulation of the discharge plasma in a Hall thruster,” Physics of Plasmas 19, 113508 (2012).
  * Fife  _et al._ (2012) 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,” in  _33rd Joint Propulsion Conference and Exhibit_ (2012) p. 3052.
  * Barral and Ahedo (2008) S. Barral and E. Ahedo, “On the origin of low frequency oscillations in Hall thrusters,” AIP Conference Proceedings 993, 439–442 (2008).
  * Hara  _et al._ (2014) K. Hara, M. J. Sekerak, I. Boyd, and A. Gallimore, “Perturbation analysis of ionization oscillations in Hall effect thrusters,” Physics of Plasmas 21, 122103 (2014).
  * Wang, Wei, and Yu (2011) C. Wang, L. Wei, and D. Yu, “A Basic Predator-Prey Type Model for Low Frequency Discharge Oscillations in Hall Thrusters,” Contributions to Plasma Physics 51, 981–988 (2011).
  * Dale and Jorns (2019a) E. Dale and B. Jorns, “Frequency scaling of the hall thruster breathing mode,” in  _AIAA Propulsion and Energy 2019 Forum_ (2019).
  * Dale and Jorns (2019b) E. Dale and B. Jorns, “Two-zone Hall thruster breathing mode mechanism, Part I: Theory,” in  _36th International Electric Propulsion Conference, University of Vienna, Austria_ (2019).
  * Dale and Jorns (2019c) E. Dale and B. Jorns, “Two-zone Hall thruster breathing mode mechanism, Part II: Experiment,” in  _36th International Electric Propulsion Conference, University of Vienna, Austria_ (2019).
  * (18) A. Smolyakov, O. Chapurin, I. Romadanov, Y. Raitses, and I. Kaganovich, “Theory and modelling of axial mode oscillations in Hall thruster,” in  _AIAA Propulsion and Energy 2019 Forum_.
  * Yamamoto  _et al._ (2002) N. Yamamoto, T. Nakagawa, K. Komurasaki, and Y. Arakawa, “Discharge plasma fluctuations in Hall thrusters,” Vacuum 65, 375–381 (2002).
  * Yamamoto, Komurasaki, and Arakawa (2005) N. Yamamoto, K. Komurasaki, and Y. Arakawa, “Discharge current oscillation in Hall thrusters,” Journal of Propulsion and Power 21, 870–876 (2005).
  * Yamamoto  _et al._ (2005) N. Yamamoto, S. Yokota, K. Watanabe, A. Sasoh, K. Komurasaki, and Y. Arakawa, “Suppression of discharge current oscillations in a hall thruster,” Transactions of the Japan Society for Aeronautical and Space Sciences 48, 169–174 (2005).
  * Hara and Kubota (2017) K. Hara and K. Kubota, “Direct kinetic simulation of ion acoustic turbulence in cathode plume,” in  _35th International Electric Propulsion Conference, Atlanta, Georgia, IEPC-2017-496_ (2017).
  * Wei  _et al._ (2015) L. Wei, L. Han, D. Yu, and N. Guo, “Low-Frequency Oscillations in Hall Thrusters,” Chinese Physics B 24, 055201 (2015).
  * Raitses  _et al._ (2006) Y. Raitses, D. Staack, A. Dunaevsky, and N. J. Fisch, “Operation of a segmented Hall thruster with low-sputtering carbon-velvet electrodes,” Journal of Applied Physics 99, 036103 (2006).
  * Romadanov, Raitses, and Smolyakov (2018) I. Romadanov, Y. Raitses, and A. Smolyakov, “Hall thruster operation with externally driven breathing mode oscillations,” Plasma Sources Science and Technology 27, 094006 (2018).
  * Fisch  _et al._ (2001) N. J. Fisch, Y. Raitses, L. A. Dorf, and A. A. Litvak, “Variable operation of Hall thruster with multiple segmented electrodes,” Journal of Applied Physics 89, 2040–2046 (2001).
  * Cohen-Zur  _et al._ (2002) 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, 4363–4374 (2002).
  * Ahedo and Rus (2005) E. Ahedo and J. Rus, “Vanishing of the negative anode sheath in a Hall thruster,” Journal of Applied Physics 98, 043306 (2005).
  * Dorf, Semenov, and Raitses (2003) L. Dorf, V. Semenov, and Y. Raitses, “Anode sheath in Hall thrusters,” Applied Physics Letters 83, 2551–2553 (2003).
  * Dorf, Raitses, and Fisch (2005) L. Dorf, Y. Raitses, and N. J. Fisch, “Experimental studies of anode sheath phenomena in a Hall thruster discharge,” Journal of Applied Physics 97, 103309 (2005).
  * Hagelaar and Pitchford (2005) G. Hagelaar and L. Pitchford, “Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models,” Plasma Sources Science and Technology 14, 722 (2005).
  * (32) “SIGLO database,” www.lxcat.net/SIGLO, [Retrived: June-2013].
  * Hagelaar  _et al._ (2004) G. Hagelaar, J. Bareilles, L. Garrigues, and J.-P. Boeuf, “Modelling of stationary plasma thrusters,” Contributions to Plasma Physics 44, 529–535 (2004).
  * Hagelaar  _et al._ (2003) G. 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, 67–75 (2003).
  * Hagelaar (2007) G. Hagelaar, “Modelling electron transport in magnetized low-temperature discharge plasmas,” Plasma Sources Science and Technology 16, S57 (2007).
  * lan (2018) “Low temperAture magNetizeD plasMA benchmaRKs,” https://www.landmark-plasma.com/ (2018), [Accessed: 26-March-2021].
  * (37) “HALL Ion Sources simulation software,” https://www.hallis-model.com/, [Accessed: 26-March-2021].
  * Wei  _et al._ (2012b) L. Wei, C. Wang, K. Han, and D. Yu, “Effect of ionization distribution on the low frequency oscillations mode in Hall thrusters,” Physics of Plasmas 19 (2012b), 10.1063/1.3676160.
  * Coche and Garrigues (2014) P. Coche and L. Garrigues, “A two-dimensional (azimuthal-axial) particle-in-cell model of a Hall thruster,” Physics of Plasmas 21, 023503 (2014).
  * Ahedo, Martinez-Cerezo, and Martinez-Sanchez (2001) E. Ahedo, P. Martinez-Cerezo, and M. Martinez-Sanchez, “One-dimensional model of the plasma flow in a Hall thruster,” Physics of Plasmas 8, 3058–3068 (2001).
  * Makowski  _et al._ (2001) 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 (2001).
  * Romadanov  _et al._ (2020) I. V. Romadanov, A. I. Smolyakov, E. A. Sorokina, V. V. Andreev, and N. A. Marusov, “Stability of ion flow and role of boundary conditions in a simplified model of the ExB plasma accelerator with a uniform electron mobility,” Plasma Physics Reports 46, 363–373 (2020).
  * Morozov and Savel’ev (2000) 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, 219–224 (2000).
  * Morozov and Savel’ev (2000) A. I. Morozov and V. V. Savel’ev, “One-dimensional hybrid model of a stationary plasma thruster,” Plasma Physics Reports 26, 875–880 (2000).
  * Hara and Hanquist (2018) K. Hara and K. Hanquist, “Test cases for grid-based direct kinetic modeling of plasma flows,” Plasma Sources Science & Technology 27 (2018), 10.1088/1361-6595/aac6b9.
  * Fernandez  _et al._ (2008) E. Fernandez, M. K. Scharfe, C. A. Thomas, N. Gascon, and M. A. Cappelli, “Growth of resistive instabilities in exb plasma discharge simulations,” Physics of Plasmas 15, 012102 (2008).
  * Koshkarov  _et al._ (2018) 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, 011604 (2018).
  * Staack, Raitses, and Fisch (2004) D. Staack, Y. Raitses, and N. J. Fisch, “Temperature gradient in Hall thrusters,” Applied Physics Letters 84, 3028–3030 (2004).
  * Kurzyna  _et al._ (2008) J. Kurzyna, K. Makowski, Z. Peradzyński, A. Lazurenko, S. Mazouffre, and M. Dudeck, “Where is the breathing mode? High voltage Hall effect thruster studies with EMD method,” AIP Conference Proceedings 993, 443–446 (2008).