domain: spi phys

text:
O Chapurin
email: alex.chapurin@usask.ca

A I Smolyakov

G Hagelaar

J.-B Boeuf

Y Raitses

Fluid and hybrid simulations of the ionization instabilities in Hall thruster

Low-frequency axial oscillations in the range 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 quasi-neutral 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 behavior and the other one, where the low-frequency oscillations coexist with high-frequency oscillations in the range 100-200 KHz, with the characteristic time scale of the ion channel flyby time, 100-200 kHz. The other test case demonstrates the effect of a finite temperature of injected neutral atoms, which has a substantial suppression 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 1972 1 ), 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.

The plasma discharge in Hall thrusters is supported by the electrons drifting in the closed azimuthal E × B direction. At the same time, ions create the thrust accelerated by the electric field in the axial direction; ions are effectively unmagnetized due to large gyroradius. 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 of great practical importance and 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 E×B electron drift are typically prone to various drift instabilities, both due to the fluid 2 and kinetic mechanisms [START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. i. insights from particle-in-cell simulations[END_REF][START_REF] Janhunen | Nonlinear structures and anomalous transport in partially magnetized E×B plasmas[END_REF][START_REF] Janhunen | Evolution of the electron cyclotron drift instability in two-dimensions[END_REF][START_REF] Villafana | 2D radial-azimuthal particle-in-cell benchmark for E×B discharges[END_REF][START_REF] Jimenez | 2D3V Particle-in-cell simulations of Electron Cyclotron Drift Instability and anomalous electron transport in E×B plasmas[END_REF] , which drive high cross-field electron currents; more details can be found in Ref. 8.

Among the plethora of wave phenomena in a Hall thruster device, low-frequency oscillations propagating in the axial direction stand as one of the most common and observed in most types of Hall thrusters 1 . They appear as the axial discharge current oscillations with frequencies of 5-50 kHz [START_REF] Esipchuk | Plasma oscillations in closed-drift accelerators with an extended acceleration zone[END_REF] . 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 its extent 1 cm and atom velocity 150 m/s gives 15 kHz. 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 → decrease of ionization → increase of the neutral density in the exhaust region → increase of the electron conductivity in that region → increase of the current and ionization → neutral depletion → 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 [10][11][12] 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 [12][13][14] . A simple model was proposed that the ion backflow region, which occurs near the anode as a result of the large contribution of the electron diffusion current (due to the density gradient) and quasineutrality constraint, provides a critical excitation mechanism for the breathing mode 13 . Linear resistively unstable modes 15 and fluctuations of electron temperature and power absorption 12,14 were also investigated as possible triggers of the breathing modes. Ref. 16 provides a general overview of physical mechanisms and stabilizing methods for 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 interrelated, 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 11,[17][18][19][20][21] . However, some calibration and adjustment of the parameters are required to satisfactorily reproduce the breathing modes characteristics observed experimentally 22 . 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. Additionally, the external circuit effects can play an important role in low-frequency plasma oscillations 14,23 .

One can propose an active control algorithms for the discharge current stabilization by varying the external circuit components 24 . Self-stabilization of the low-frequency discharge current was observed when the electromagnet current is coupled to the discharge current 25 .

In this paper we do not consider these effects, focusing on the plasma dynamics not including the external circuit.

While many time-dependent numerical models for breathing modes are based on fully fluid formulations, the hybrid modeling was also undertaken using the kinetic description for ions and neutrals 15,19,20,[26][27][28] . The extent to which the ion and neutrals kinetic effects influence the breathing mode excitation and characteristics remain a mute point of many studies. One of the goals of this paper is the analysis of 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 ions and atoms with the two fluid moments, conservation of mass and momentum, and electrons are considered in a drift-diffusion approximation with a full electron energy balance. In the hybrid model, ions and atoms are kinetic via particlein-cell method, and electrons are fluid, modeled in the same way in both approaches. Both models include plasma recombination at the anode and neutral dynamics with ionization due to electron-neutral impact. The 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 framework 29 is used. The hybrid code was developed in the LAPLACE laboratory, France 18,30,31 . This work is based on one of the test cases in the LANDMARK (Low temperAture magNetizeD plasMA benchmaRKs) project 32 and follows the setup initially proposed in LANDMARK.

Previously, a comparison between fluid and hybrid models for the axial direction of Hall thruster configuration was presented in Ref. 20. 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 ion transition through the ion-sound barrier.

Nor did it include the full electron energy balance.

One of the important findings 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 of ∼14 kHz that coexists with the high frequency ion "transient-time" oscillations of ∼150 kHz 9 . In the second regime, with low electron energy losses, pure 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 high frequency oscillations as the excitation of the resistive modes (convective instability with the characteristic ion fly-by frequency) 15,[33][34][35] . 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) 34 . Similar features are shown in this work, while 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 that of the breathing modes at the lower end of its spectrum. Some axial thruster models with ionization but without electron diffusion [START_REF] Morozov | Fundamentals of stationary plasma thruster theory[END_REF]36 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 according 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 13 . It was demonstrated that such configuration is prone to low-frequency oscillations, where the ion backflow region is necessary. Recently, similar conclusion has been reached in Ref. 37 where it was formally proved that the sign-alternating ion velocity profile with a positive slope, i.e. negative ion velocity near the anode and positive towards the exit, indeed is a necessary condition for the excitation the oscillations.

Another important result is the demonstrating the effect of atom temperature in the solo regime. We find that a finite energy spread of injected atoms strongly suppresses the amplitude of the oscillations compared to the injection cases with the same velocity, further called 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 the three test cases with a detailed comparison between the models.

II. FLUID AND HYBRID MODELS OF LOW-FREQUENCY DYNAMICS

This section presents a detailed description of the full fluid model and the hybrid model.

Two models share the same electron fluid equations (drift-diffusion approximation and energy evolution). The models are considered in the electrostatic and quasi-neutral 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, and radial atom losses were not included.

The models also include 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 cm is assumed in the axial direction of a Hall thruster (xdirection), with the channel exit in the middle where the radial magnetic field has its maximum, Fig. 1. The profile of the magnetic field magnitude given by B = B 0 exp -(xx 0 ) 2 /2δ 2 B , where x 0 = 2.5 cm is the channel exit location, and δ B is the characteristic width coefficient for the magnetic field profile, which are set δ B,in = 1.1 cm, δ B,out = 1.8 cm, respectively for the inner and outer regions, according to the LANDMARK setup 32 . The chosen magnetic field profile corresponds to a typical radial magnetic field of the SPT 100 thruster [START_REF] Morozov | Fundamentals of stationary plasma thruster theory[END_REF]19,21 , with the maximum near the channel exit.

6

A. 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, a constant flow velocity V a along the channel is considered, and the continuity equation with the source term is used to describe their dynamics:

∂n a ∂t + V a ∂n a ∂x = -βn a n e , (1)
where n a is the atom number density, n e is the electron number density, β(ε) is the ionization rate coefficient that depends on the electron energy ε = (3/2)T e , where T e is the electron temperature (β is shown in Appendix B). The ionization due to electron-atom impact produces a pair of ion and electron with a loss of neutral atom, hence the sink term -βn a n e (same and opposite sign source terms are included in the ion and electron continuity equations).

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:

∂n i ∂t + ∂ ∂x (n i V i ) = βn a n e , (2)
∂V i ∂t + V i ∂V i ∂x = e m i E + βn a (V a -V i ) , (3)
where n i is the ion density, V i is the ion flow velocity, E is the axial electric field, e is the scalar elementary charge, m i is the ion mass (Xenon, 131.293 amu). 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). Note that the radial sheath losses are often included in low-frequency models of a Hall thruster 21 . Such losses are of the same order of magnitude as the ionization losses. We neglect this effect in the present study, focusing on the main coupling mechanism between ions and neutrals, such as the ionization and the recombination at the anode.

The magnetized electron species are described with the first three fluid moment equations:

(electron inertia is neglected)

∂n e ∂t + ∂ ∂x (n e V ex ) = βn a n e , (4) 0
= - e m e E - eB m e (V e⊥ × ẑ) - 1 n e m e ∂ (n e T e ) ∂x -ν m V e⊥ , (5) 3 2
∂ ∂t (nT e ) + 5 2

∂ ∂x (n e V ex T e ) + ∂q e ∂x = -n e V ex ∂φ ∂x -n e n a K -nW, (6)
where n e is the electron density, V e⊥ = (V ex , V eθ ) is the electron flow velocity perpendicular to the magnetic field (x and θ are the axial and the azimuthal coordinates), m e is electron mass, B is the external radial magnetic field, ν m is the total electron momentum exchange frequency, W is an anomalous energy loss coefficient, K is the collisional energy loss coefficient (it is shown in Appendix B) and q e is the electron heat flux. A phenomenological anomalous electron energy loss coefficient (e.g., due to radial sheath losses) W is introduced

19 as W = ν ε ε exp (-U/ε), (7)
where ε = 3T e /2, U = 20 eV, ν ε is the anomalous energy loss frequency. The range of ν ε = (0.3 -1.4) × 10 7 s -1 is found to adequately reproduce experimentally observed electron temperature profiles, with the average peak electron temperature in the range (24-38) eV.

The heat flux across the magnetic field is

q e = - 5 2 µ e nT e ∂T e ∂x . (8)
The electron momentum conservation equation ( 5) is simplified assuming no pressure gradients nor equilibrium electric fields other than in the axial direction, hence the axial electron velocity (denoted further as V e ) can be expressed as

V e = -µ e E - µ e n e ∂(nT e ) ∂x , (9)
where the electron mobility µ e is the well-known classical electron mobility across the magnetic field:

µ e = e m e ν m 1 1 + ω 2 ce /ν 2 m , ( 10
)
where ω ce = eB/m e is the electron cyclotron frequency. Eq. ( 9) commonly called the driftdiffusion equation. The model of electron transport as based on the assumption of the following total electron momentum exchange collision frequency:

ν m = ν en + ν walls + ν B , (11)
where the electron-neutral collision frequency ν en , electron-wall collision frequency ν walls , and anomalous Bohm frequency ν B are given with

ν en = k m n a , (12)
ν walls = α10 7 [s -1 ], (13)
ν B = (β a /16) eB/m e . ( 14
)
where k m = 2.5 × 10 -13 m -3 s -1 , α and β 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 α in = 1, α out = 0, the anomalous contribution is set to β a,in = 0.1, β a,out = 1. The anomalous loss (electron) energy frequency coefficient in Eq. ( 7) is a constant, set to ν ε,out = 10 7 s -1 outside the channel. Inside the channel ν ε,in is variable and serves as the input parameter of the model, see Section III.

The step-like behaviour of some model parameters (inside and outside the channel) can lead to numerical difficulties, but the obtained results were tested for convergence and the problems presented in this paper. A possible advantage of implicit approach is resolution of the Debye sheath (although approximately), where the sheath potential drops over one computational cell, assuming the cell size is larger than the Debye length 45 . We believe that bulk oscillations in Hall thruster, which is the main object of this study, are well reproduced with the full quasi-neutral approach given above.

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 equation ( 9). Full quasineutrality n = n i = n e is enforced and the self-consistent electric field is found via the electron driftdiffusion equation ( 9), given by

E = J T enµ e - V i µ e - 1 n ∂nT e ∂x , (15)
where the total current density J T = en (V i -V e ). Here J T is constant in space (divergenceless current), which can be seen by combining continuity equations for ions and electrons with quasi-neutral assumption. We will use the integral approach, consisting of the evaluation the total current density J T via the constraint L 0 Edx = U 0 (then it is substituted to Eq. ( 15) to evaluate the electric field), which yields to

J T = U 0 + L 0 V i µ e + 1 n ∂p e ∂x dx L 0 dx enµ e , ( 16
)
where L is the system length, U 0 = 300 V is the applied voltage (without sheath voltage).

In this model, the effects of external circuit and fluctuations of the applied voltage are not not included.

The fluid model is solved with the following boundary conditions. A constant mass-flow rate ṁ = 5 mgs -1 and the full recombination of plasma that flows to the anode determines the value of n a at the anode boundary,

n a (t, x = 0) = n a0 + n a,rec = ṁ m i AV a - n i V i (t, x = 0) V a , ( 17
)
where

A = π(R 2 2 -R 2 1
) is the anode surface area, with the inner and outer radii R 1 = 3.5 cm and R 2 = 5 cm. The value of the ion velocity is imposed at the anode as V i (t, x = 0) = -b v T e /m i with the parameter b v = 0-1, the Bohm velocity factor that can be varied. Further, we show that Case 2 is highly sensitive to the parameter b v , and other cases are insensitive. In general, the electron energy balance (6) with the heat flux included requires requires two boundary conditions. Electron temperature near the anode observed in experiments 43,44 is of few electron volts, thus we have fixed both anode and cathode electron temperatures to T e (t, x = 0) = T e (t, x = L) = 2 eV to stay close to the experimental conditions 43 .

As noted in Ref. 17, 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 = 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 current 46 . Note that the total current J T , given by Eq. ( 16), is a function of U 0 with the integral dependence on all main plasma parameters, thus the problem is inherently nonlocal, which has no analogy with the standard de 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).

B. 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 18,30,31 . 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 sampling of ionization events (via the null collision method 45,47 ) using the macroscopic ionization rate β(ε), obtained in the same way as in the fluid model.

Neutral atoms are injected with a constant flow rate ṁ either with the constant velocity

V a (monokinetic), thus f (v x ) = δ(v x -V a )
, or with a half-Maxwellian velocity distribution function at the left wall (anode),

f (v x ) = 2v x v 2 T a exp(- v 2 x v 2 T a ), v x > 0, ( 18
)
where v 2 T a = 2k B T a /m a is the atom thermal speed (m a = m i ), T a is the atom temperature (in K). For the half-Maxwellian injection (18) the average flow velocity is v T a / √ π. Ions are assumed singly charged and unmagnetized, thus they are only accelerated by an electric field. Ions are produced according to the ionization rate coefficient β(ε), i.e. self-consistently with the electron temperature evolution and the local atom density. Both atoms and ions are lost at the boundaries. In this model the ion velocity is not forced to the Bohm velocity, like in the fluid ion model. With the quasi-neutral 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:

∂f i ∂t + v ix ∂f i ∂x + e m i E ∂f i ∂v = S(x, v x ), ( 19
)
∂f a ∂t + v ax ∂f a ∂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

III. SIMULATION RESULTS AND COMPARISON

The fluid and the hybrid models described above were studied for three test cases (denoted as Cases 1, 2, and 3). Cases 1 and 2 will demonstrate two distinct regimes of low-frequency oscillations, and Case 3 shows the effect of atom temperature. Note that Cases 1 and 2 use monokinetic atoms with the velocity V a = 150 m/s. Cases 1 and 2 are chosen with the following observation: larger values of anomalous (electron) energy loss frequency ν ε,in in Eq. ( 7) 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 values of ν ε,in reveal only the large amplitude low-frequency oscillations, we call it the solo regime. In Case 1 ν ε,in = 0.95 × 10 7 s -1 , and in Case 2 ν ε,in = 0.4 × 10 7 s -1 . This is the only parameter distinguishing Cases 1 and 2 (for all cases reported in this paper we keep ν ε,out = 10 7 s -1 ). The difference in the electron temperature distributions might be due to the difference in the electron flow velocity near the anode. Fig. 3 shows the time-averaged electron velocity components, evaluated according to Eq. ( 9). Note substantially higher electron velocity due to the pressure gradient in Case 1 for x < 0.4 cm, which also results in larger total electron flow velocity near the anode. Thus, in Case 1, convection reduces the electron temperature in the ionization region (increasing it near the anode). Indeed, the average gradient parameter

L -1 n = ∂ x n i /n i
is about five times larger for Case 1, which drives this electron current. Recall that the presheath region is formed due to diffusive electron current by inducing ion current (total current is conserved), generating the negative electric field in this region. All features presented above clearly distinguish Case 1 and Case 2. Finally, Case 3 demonstrates the effect of finite atom temperature, where atoms are injected into the system with the half-Maxwellian distribution, Eq. ( 18), and all the other parameters are the same as in Case 2. The main effect of atom temperature is found to be 14 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. 9). Small averaging window of length 2.3 µs is applied to the ion and electron currents (to filter out high-frequency noise).

PLEASE CITE THIS ARTICLE AS

model contains more noise (statistical noise due to the use of macroparticles), but the highfrequency component is clearly seen at around 125 kHz. In the fluid model the high-frequency mode is shifted towards higher frequencies, centered at about 175 kHz. Besides the currents, a more rigorous comparison between the two models is shown in 16

Figs. 7a-7d, with the time-averaged axial profiles of the main 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 main discrepancy lies in the peak plasma density in the ionization (source) region at about 1.3 cm from the anode; the hybrid model results in the higher value of the peak plasma density. Also, the ion velocity in the plume (x > 2.5 cm) is slightly higher in the hybrid model. 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 nullified pressure term gives the closure in the ion fluid model). 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 when accelerating by the axial electric field. To check the validity of this assumption, 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:

∂V i ∂t + V i ∂V i ∂x = e m i E - 1 n i ∂p i ∂x + βn a (V a -V i ) , (21)
where the fluid moments, such as the ion density n i , the ion flow velocity V i , and the ion pressure p i were evaluated from the ion distribution function in the following way:

n i = f i dv ix , (22)
V i = 1 n i v i f i dv ix , (23)
p i = m i v i v i f i dv ix , (24)
where v i = v i -V i is the random component of particle velocity. The time-averaged profile of the ion temperature T i = p i /n i is shown in Fig. 8b, revealing values of the ion temperature up to 3.5 eV in the ionization region, with the average over the whole domain of 1.7 eV.

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 between the left-and the right-hand sides of equation ( 21).

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 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. Ion momentum terms, norm. It is interesting to further inspect each term in Eq. ( 21). The unsteady term ∂V i /∂t is negligible due to time averaging, while the other terms have comparable values inside the channel, see Fig. 9b. In the near-anode region, the ion pressure (T i ≈ 0.1 eV) and the ionization friction are negligible, so the ions accelerate towards the anode in the weak negative electric field ballistically. Then in the ionization (source) region, we see that all terms are comparable. The ion pressure and the collisional drag compensate the ballistic acceleration, so the ion inertia remains low. Finally, in the acceleration zone (x > 1.5 cm) and beyond the inertial and the ballistic terms start to dominate. It is seen that the ion pressure term changes the sign at x ≈ 1.6 cm (due to ion density profile, Fig. 7b) and contributes to the ion acceleration.

v i ∂v i ∂x e m i E βn a (v a -v i ) 1 n i m i ∂p i ∂x (b)
Based on the results above, the ion pressure force term -T i ∂ x n i was added to the fluid model with the temperature kept constant for simplicity. For the self-consistent treatment, an ion energy balance equation shall be included (or equation of state). The average in space and time value of the ion temperature evaluated in the hybrid model for this case is 1.7 eV. However, it is found to damp the oscillations and produce a fully stationary solution with the ion pressure force mentioned above. With a lower value of T i = 1.2 eV, the oscillatory solution is obtained. This results in a better agreement between the two models, the total current amplitude decrease and the main low-frequency mode increase, Figs. 10a and10b. The low-frequency mode increased to 13.9 kHz, and the high-frequency peak is shifted to a lower value of 153 kHz. The ratio of the ion current to the total current shows an improved agreement, 45%, the same as in the hybrid model. However, the time-averaged total current value of 8.7 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. As in the previous case, the total current amplitude is higher in the fluid model, Figs. 11a and11b, but the main oscillation frequency in two models is similar, Figs. 12a and12b.

21

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269

The time-averaged profiles agree well between the two models, with the only notable discrepancy in the plasma density profile, which is ∼30% lower inside the channel region for the fluid model, see Figs. 13a-13d. 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 ≈ 0).

Indeed, the average (space and time) ion temperature, calculated similarly as in Case 1 from the hybrid model, is 0.3 eV (about 5 times smaller than in Case 1) and not exceeding 1 eV in the domain.

The main stages of the breathing mode dynamics are illustrated with the plasma and atom density distributions at various moments in time, It is found that the oscillation dynamics in Case 2 highly depends 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, which is seen in Fig. 17. In contrast, 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. We modeled this behaviour with the free boundary condition for the ion velocity at the anode in the fluid model, enforcing ∂ 2 x V i (t, x = 0) = 0. In the presence of ion backflow (like we see in Case 1,2) in the quasi-neutral approximation, the ion velocity at the anode is defined self-consistently by the characteristic at the location where V i = 0; thus, a fixed boundary is not formally required. Moreover, depending on operating conditions, various sheath regimes can form in a Hall thruster 43 , thus the ion velocity can differ from the Bohm velocity at the anode.

The obtained result in the fluid model with the free boundary condition 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. 18, closer to the hybrid mode result. The time-averaged atom and plasma density profiles are in a better agreement, Figs. 19a and19b.

C. Case 3: Effects of finite temperature of neutral atoms

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 K. Unlike Case 2, the oscillation amplitude in this case is much smaller in the hybrid model, Fig. 21a, with the amplitude of ∼1 A (in comparison to Case 2 with ∼13 A). In the previous cases, injected atoms in the hybrid model were kept monokinetic (zero thermal energy). The average injection velocity for the half-Maxwellian is V 0 = v T a / √ π = 142 m/s, which is close to 150 m/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 (the ionization frequency does not depend on the atom particle velocity, but their mean free path does), observed both in experiments 49,50 and simulations 51 . 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:

∂ (n a V a ) ∂t + ∂ ∂x n a V 2 a = -βn a n i V a - T a m i ∂ (n a ) ∂x , (25)
where the closure is given with the constant temperature T a = 500 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. 20. Thus, Eq. ( 25) along with the continuity equation, these oscillations are the so-called resistive modes 15,35 . We believe these modes are observed in the models without electron pressure, many of which were done in the primary studies and modeling of the axial dynamics of Hall thrusters [START_REF] Morozov | Fundamentals of stationary plasma thruster theory[END_REF]15 . We conjecture that the resistive modes of high frequencies play an important role in the excitation of the low-frequency modes.

∂n a ∂t + ∂ (n a V a ) ∂x = -βn a n i , (26
For Case 2, which shows only the low-frequency oscillations, we already performed extensive studies, Ref. 13, 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 the recombination from the anode, Eq. ( 17), in comparison to the hybrid model where this velocity was unconstrained and found to be smaller on average (thus leading to smaller atom yield). Replacing the Bohm condition with a free boundary condition in the fluid model 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. The fluid model extended with the atom momentum balance equation ( 25) with a constant temperature resulted in a better agreement between the two models, reproducing the selective ionization effect and predicting a stationary solution. The hybrid model preserved the breathing oscillations (with similar frequency to Case 2) but with a much smaller amplitude. We were not able to find higher amplitude low-frequency oscillations for this case in the hybrid model (scanning parameters ν ε or β a ). However, a small variation in the position of the maximum magnetic of the field profile (shifting it inside the channel), higher amplitude breathing modes were observed with thermal atoms. This important behaviour needs further attention and is left for future studies.

The magnitude and profile of the anomalous mobility have strong effect on the characteristics of the breathing mode. For the purposes of the comparison of the fluid and hybrid models we used the empirical values of the anomalous collisions frequency which were often used in previous works 14,18,22,32 , so that the anomalous mobility follows the magnetic field profile according the equations (10-14). Our simulations reveal strong sensitivity of the breathing oscillations on the location of the magnetic field maximum. Anomalous electron mobility and diffusion remains the critical parameters that are poorly known in theory and experiments. Further work on calibrating the models against experimental data is required 22,52 .

Appendix A: Resistive modes

Here we present the cases with more prominent resistive modes and show their effect on the ion heating, which is the extension to Case 1. It is found that the larger value of the anomalous (electron) energy loss frequency, ν ε,in = 1.2 × 10 7 s -1 , leads to the solution where the resistive mode dominate, see Fig. 24a. It reveals larger amplitude and clearly distinct high-frequency oscillations (168 kHz), with the small low-frequency modulation (14.4 kHz), Fig. 24b shows 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. 23, where the resistive modes appear for values ν ε,in = 0.9 s -1 and higher. Now, we show that the frequency of the resistive mode scales with the average electron collision frequency νm , Fig. 25a, and thus the higher electron mobility ∼1/ν m leads to the higher frequencies. At the same time, the ion heating effect is stronger for the higher frequency of resistive modes, Fig. to solve the algebraic-differential system resulting from the method of lines utilized in the BOUT++. The upwind terms are discretized with the WENO 3 rd order scheme and other terms with the 2 nd order central-difference scheme.

Figure 1 .

 1 Figure 1. The magnetic field profile used in simulations, with the channel exit located at the middle, 2.5 cm from anode (dashed line).

  main behaviour is reproduced in both models (fluid and hybrid). Moreover, the fluid model implemented in the BOUT++ was tested against the fluid code MAGNIS[38][39][40] for the same model resulting in a perfect match between two. Numerical details and convergence tests for the fluid model are presented in Appendix C.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 check the validity of the quasi-neutral approximation 17 by using the Poisson equation ε 0 ∂E/∂x = e(n in e ) (ε 0 is the permittivity of free space). With the typical values of the electric field E ≈ 10 4 V/m, the size of acceleration zone ∼1 cm, average plasma density n 0 = 10 17 m -3 , one obtains the difference (n in e )/n 0 ≈ 5 × 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 in the Debye sheath near the anode, it still allows the presheath region to form if the electron pressure is included[41][42][43][44] , 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. An alternative to the fully quasi-neutral approach used in this paper is solving the Poisson equation with implicit (or semi-implicit) methods, which allows to ignore the Debye length scales, since full resolution of Debye scales would be excessive for the large-scale

V

  a (constant) and the source term can be expressed 19 as S(x, v x ) = βn e n a δ(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 techniques 48 . Solutions to Eqs.(19,20) effectively obtained by solving the motion equations for the corresponding particle type (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.In the hybrid model, in all test cases presented below, the average total number of macroparticles in the system was ∼5 million for the ions and ∼0.5 million for the atoms. The fluid equations (for electrons) in the hybrid model are resolved with the finite volume methods 45 , with the same spatial discretization and time step as the PIC method, set to 300 cells and ∆t = ∆x/v max , where v max = 5 2eU 0 /m i , resulting in ∆t = 1.62 ns for U 0 = 300 V. The same time step is used as the collisional time step in the Monte Carlo process, and for a given ∆t the condition ν ion ∆t 1 is well satisfied for a typical value of the ionization frequency ν ion ≈ 5 × 10 5 s -1 (T e = 20 eV, n a = 10 19 m -3 ).

Figure 2 .

 2 Figure 2. Averaged in time profiles of ion velocity (a) and electron temperature (b) for Case 1 (ν ε,in = 0.95 × 10 7 s -1 ) and Case 2 (ν ε,in = 0.4 × 10 7 s -1 ). The result is obtained with the hybrid model.

Figure 3 .

 3 Figure 3. Electron flow velocity components given by Eq. (9) for Case 1 and 2, averaged in time (denoted with angle brackets) over few periods of the low-frequency component (result obtained with the hybrid model).

Figure 4 .Figure 5 .

 45 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 loss coefficient ν ε (result obtained with the hybrid model). The range of ν ε shown here is found to correspond adequately to the electron temperature profiles observed in experiments. For example, for values ν ε < 0.3 × 10 7 s -1 the peak electron temperature reaches values of higher than 50 eV.

Figure 6 .

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

Figure 7 .

 7 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), for Case 1.

Figure 8 .

 8 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).

Figure 9 .

 9 Figure 9. Difference between left and right hands sides of the ion momentum balance equation (22) 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).

Figure 10 .

 10 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 with the ion pressure term included.

Figure 11 .

 11 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 cm.

Figure 12 .

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

Figure 13 .

 13 Figure 13. 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), for Case 2.

Fig. 15 .

 15 After atoms reach the ionization zone and undergo ionization, plasma density increases and is quickly depleted (∼1 km/s) to the left (due to the backflow region with negative ion velocity) and the right of the ionization zone. Then all of the ion flux that reached to the anode recombines and forms the peak in neutral density at the anode, increasing the number of atoms advecting to the ionization zone; this process repeats. In our previous work, Ref. 13, this setup was studied in detail with the reduced fluid model, consisting of the continuity equations for neutrals and ions with a prescribed ion velocity profile. It was shown that the ion backflow region near the anode is sufficient to excite the low-frequency oscillations, and the recombination is not necessary. The role of the recombination in this case consists in increasing the oscillations amplitude. , [km/s] f i (x, v x )/10 10 , [s/m 2 ]

Figure 14 .

 14 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).

Figure 15 .

 15 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 .Figure 17 .

 1617 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 . The result is obtained with fluid model.

Figure 18 .

 18 Figure 18. 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.

Figure 19 .

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

  Fig. 21b, but at the same time, the fluid model resulted in a completely stationary solution, Fig. 21a. To better identify the role of the atom distribution in the hybrid model, we configured the shifted Maxwellian distributions 48 for atoms with various temperatures denoted as T a,s , and

Figure 21 .

 21 Figure 21. The total current in fluid and hybrid simulations (a); comparison of time-averaged spatial profile of atom flow velocity for both models (b).

Figure 22 .

 22 Figure22. 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).

Figure 23 .Figure 24 .

 2324 Figure 23. 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.

Figure 25 .

 25 Figure 25. Frequency of breathing and resistive modes as a function of averaged (over time and space inside the channel) total electron momentum exchange frequency νm (a). Average ion temperature (in time and whole domain) as function of νm (b). The parameter varied in this study was the anomalous Bohm coefficient β a,in inside the channel, directly affecting the total electron frequency, Eq. (11).

Figure 26 .

 26 Figure 26. Ionization rate β and electron ionization energy loss coefficient K.

  scheme:

  Figure 27. Ion density spatial profile, averaged over 10 time periods of its main low-frequency mode, where the parameters corresponds to Case 2. Convergence with grid size (a) with the fixed ODE solver tolerances: relative is 10 -5 and absolute is 10 -9 . Convergence with the CVODE tolerances (b) with the fixed number of grid points: n x = 256.

  This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

  This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.∂ t (n a v a ) + ∂ x n a v 2 a -∂ x (nT ax ) / (m i n a )This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269

ACKNOWLEDGMENTS

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

kHz

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269

DATA AVAILABILITY

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

Appendix B: Ionization rate and energy loss functions

The ionization rate β(ε) = σ(v)v as a function of electron energy ε is obtained by averaging over Maxwellian EEDF with the cross-sections from the SIGLO database 53 . The collisional energy loss coefficient K(ε) is generated by BOLSIG+ 54 as a table-valued function.

Both functions are shown in Fig. 26.

Appendix C: Numerical methods for fluid model

The time-dependent fluid model, given by equations (1,2,3,6), is solved via the BOUT++ 29 computational framework. The CVODE solver from the SUNDIALS package 55 was used

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. -11 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0094269