Numerical and Experimental Investigation of Longitudinal Oscillations in Hall Thrusters

: One of the main oscillatory modes found ubiquitously in Hall thrusters is the so-called breathing mode. This is recognized as a relatively low-frequency (10–30 kHz), longitudinal oscillation of the discharge current and plasma parameters. In this paper, we present a synergic experimental and numerical investigation of the breathing mode in a 5 kW-class Hall thruster. To this aim, we propose the use of an informed 1D fully-ﬂuid model to provide augmented data with respect to available experimental measurements. The experimental data consists of two datasets, i.e., the discharge current signal and the local near-plume plasma properties measured at high-frequency with a fast-diving triple Langmuir probe. The model is calibrated on the discharge current signal and its accuracy is assessed by comparing predictions against the available measurements of the near-plume plasma properties. It is shown that the model can be calibrated using the discharge current signal, which is easy to measure, and that, once calibrated, it can predict with reasonable accuracy the spatio-temporal distributions of the plasma properties, which would be difﬁcult to measure or estimate otherwise. Finally, we describe how the augmented data obtained through the combination of experiments and calibrated model can provide insight into the breathing mode oscillations and the evolution of plasma properties.


Introduction
Hall thrusters are the most widely adopted electric propulsion technology for space applications. Their success is largely due to their high thrust efficiency and thrust-topower ratio, coupled with an overall robust design. During Hall thruster operation, a stream of electrons is generated by an external hollow cathode while neutral gas, typically xenon, is injected in the thruster annular ceramic channel through an anode. A magnetic circuit generates a predominantly radial magnetic field in proximity to the exit section of the discharge channel. The intensity of the magnetic field is tailored to ensure the magnetization of the electron population while leaving the ions unaffected. Part of the electrons generated at the cathode move towards the anode due to the potential difference imposed externally between the two. In the region where the magnetic field is high, the crossed electric and magnetic fields force the electrons to perform an E × B motion in the azimuthal direction. When the propellant neutral atoms reach the region of high azimuthal electron current, they are ionized by electron impact. The electrons generated by the ionization events join those coming from the cathode and drift towards the anode through successive collisions. The ions are instead accelerated by the applied electric field and exit the thruster at high speeds, generating thrust. Finally, a second population of electrons is emitted by the cathode to neutralize the outbound ion beam and preserve the charge balance of the system. More details on the physical processes involved in the plasma dynamics in Hall thrusters can be found, for instance, in [1,2].
Several oscillatory modes exist inside the plasma in the channel and in the near-plume of Hall thrusters. The characteristics and nature of these oscillations strongly depend on the thruster geometry, magnetic field topology and operating condition [3][4][5]. One of the most prominent unsteady modes of Hall thrusters is recognized as a low-frequency  oscillation of the plasma properties in the longitudinal direction, which produces periodic rises in the thruster's discharge current. This so-called breathing mode was first described in the 1970s [6] and, since then, has been reported almost ubiquitously in Hall thrusters.
An intuitive understanding of the breathing mode can be traced back to the periodic replenishment of the channel by the injected propellant, followed by a rapid ionization and subsequent acceleration of the generated plasma. The ejected plasma generates a surge in the discharge current and simultaneously depletes the channel of the propellant particles before the cycle starts again with the injection of new neutral atoms through the anode. This intuitive description falls short of capturing the nature of breathing mode in its entirety, and the core physical origin of the instability is still a subject of debate and further investigations. This is partially due to the complex interaction of the longitudinal modes with other oscillations that generate anomalous diffusion of the electrons through the magnetic field [1] and partially because of the difficulties in gathering high-frequency experimental data on the plasma properties inside Hall thrusters.
From a numerical perspective, longitudinal oscillations in Hall thrusters have been investigated since the advent of plasma simulations in the field [7,8]. Numerical simulations are carried out mainly using two different modeling approaches, i.e., the fluid and kinetic descriptions. Moreover, the kinetic approach can be in a direct kinetic (DK) or particle-in-cell (PIC) formulation. Several advantages can be obtained with hybrid models that combine different approaches for different species within the same computational framework. Although complex methodologies are available and currently used to carry out 2D (see, for instance, [9][10][11][12][13][14][15][16][17][18]) or even 3D simulations (see, for instance, [19,20]) of Hall thruster discharges (see also [21,22] for recent reviews of the literature), 1D models are still very appealing for studying the fundamental mechanisms in the dynamics of the plasma discharge. Indeed, despite the unavoidable lower accuracy and prediction capabilities in comparison with more complex models, they are easier to handle and definitely cheaper from a computational viewpoint. Moreover, due to their simplicity, it is possible to apply ad-hoc tuning against reference data, for instance, when they are available from experiments or from more complex simulations, so as to reproduce at least a set of operating points in a fairly accurate way, thus partially compensating the inaccuracy related to the simplification of the physical model employed. In the literature, several 1D models, comprising both hybrid kinetic-fluid models [23][24][25] and fully-fluid models [26][27][28][29][30][31][32], have been successfully employed to simulate the global characteristics of plasma discharges in Hall thrusters. Among the two modeling approaches, fully-fluid models are easier to tune against reference data, as they are entirely formulated in terms of deterministic spatiotemporal plasma characteristics. In hybrid approaches, in any case, tuning can be carried out on the fluid part of the model, as proposed, for instance, in [33]. Moreover, fully-fluid models are particularly suitable to be used for stability analyses, as they can be linearized in a rather straightforward way. As a result, most of the stability analyses documented in the literature, performed to further investigate the nature of the breathing mode, are carried out starting from simplified dynamical systems, as in [34], or with fully-fluid models (see, for instance, [35,36]).
From the experimental perspective, gathering data on the plasma parameters in the channel and near-plume of Hall thrusters poses significant technological challenges. This is mainly due to the harsh environment that invasive probes need to withstand, limiting the maximum residence time in the high plasma density regions of the domain. This is especially true when high-frequency data need to be collected in order to reconstruct oscillatory modes, since long time series are necessary. Additionally, invasive probes can significantly disturb the plasma flow and invalidate the measurements. Therefore, significant care must be taken in the design of the probe and of the test setup. For what concerns non-invasive probes, such as optical sensors, the annular geometry of Hall thrusters and the intrinsically complex test setup hinder their applicability. Sensors installed on rapidly moving arms to probe the plasma in Hall thrusters have been used since the late 1990s [37]. Since then, similar concepts have been proposed, including single and triple Langmuir probes, as well as emissive probes, in order to gather a time-averaged picture of the plasma properties in the thruster near-plume [38][39][40][41][42][43]. The plasma perturbations occurring upon probe insertion complicate the interpretation of experimental data, particularly within the channel and in the vicinity of the magnetic peak, where the electron azimuthal current is the highest. The factors influencing the onset and the impact of these perturbations were studied by several authors, seeking to improve the reliability of the measured data and the quality of the information they provide. Haas et al. [37] identified material ablation of the probe as an important source of perturbation and suggested to minimize the probe residence time in the hotter plasma regions to reduce probe perturbations. Nevertheless, perturbations remain up to a certain degree and, as discussed by Jorns et al. in Reference [42], a sharp transition of global plasma parameters and a downstream shift in plasma properties are observed upon probe insertion in the channel. Optical methods, such as Laser Induced Fluorescence (LIF) [44][45][46][47][48][49] and Thompson scattering [50], provide non-intrusive alternatives to fast moving ionic probes for the exploration of the steady and unsteady behavior of the near-plume plasma of Hall thrusters. In spite of being capable of providing useful information, these diagnostic systems can present signal-to-noise interpretation issues and, in most cases, require complex and expensive setups, limiting their extensive use. Lobbia et al. made noteworthy contributions to the use of invasive diagnostics to study the evolution of low-frequency oscillations in Hall thrusters [51][52][53][54], positioning a single Langmuir probe (with a noise compensating null probe) in multiple fixed locations of the thruster plume and rapidly biasing it to obtain distributions of plasma properties on a 2D grid. Since single probes require time to undergo several bias cycles, measurements were limited to the far plume, where the probe could withstand the plasma conditions indefinitely.
The present paper proposes an investigation of the breathing mode in a 5 kW-class Hall thruster, SITAEL's HT5k DM2, based on the synergic combination of a 1D fully-fluid model and experimental data. The latter comprises two different datasets: (a) the signal of the thruster discharge current acquired with an oscilloscope; (b) the measurement of the local plasma properties in the channel and near-plume carried out with a fast moving triple Langmuir probe and processed with a dedicated diagnostic method [55]. The numerical model uses standard modeling methodologies and contains physically meaningful free parameters associated with those modeling aspects that are mostly affected by uncertainty. Within this context, the main contribution of the paper is to assess the potential of the considered 1D model when the free parameters are calibrated against easily measurable experimental quantities. To this purpose, the model has been calibrated using only the first of the two datasets, the current measurement, which does not require a complex experimental setup and is easy to acquire. Results show that calibration is indeed possible, leading to an accurate recovery of the time evolution of the discharge current. Successively, thanks to the second dataset, i.e., the measurement carried out with the triple Langmuir probe, it has been possible to assess the level of accuracy reached by the model in the prediction of the values and oscillations of the plasma properties in the thruster nearplume when calibration is carried out solely on the current signal. Results show that the model provides a good estimation of the plasma properties, which, conversely, are difficult to measure. As a result, the proposed calibrated model can also be interpreted as an extrapolator of partial experimental data in Hall thrusters. Thus, more than a predictive tool, it can be used as a complementary tool, together with experiments, to estimate quantities that cannot be measured and to orient more detailed measurements in the experiments. Finally, in the manuscript we provide an example of how the augmented data resulting from the joint use of experiments and calibrated model can be used to explore the variability of the plasma properties within typical cycles associated with the breathing mode.
Concerning the outline of the paper, Section 2 describes the numerical formulation adopted to model the plasma flow in Hall thrusters, discussing the main assumptions and presenting the core equations of the model. Section 3 details the characteristics of the test item and summarizes the experimental setup and data processing technique, described at length in [55]. The comparison between the results of the calibrated plasma simulations and the experimental results is then presented in Section 4. Finally, Section 5 summarizes the conclusions of the present research.

Numerical Description of the Plasma Dynamics
A fully-fluid unsteady 1D model is proposed here to describe the dynamics of the plasma discharge in the thruster. The plasma is assumed to be quasi-neutral and mainly composed of singly-charged ions. Three species are thus included in the model, namely electrons, singly-charged ions and neutrals. A drift-diffusion approximation is used to treat the electron flux. Ions are supposed isothermal and unmagnetized, and neutrals are assumed to have constant axial velocity (u n ) in the domain. Furthermore, the magnetic field is assumed to be purely radial. With reference to Figure 1, the model is aimed at providing time-varying, averaged plasma properties over z-sections inside the channel and in the near-plume, from the anode surface (z = 0) up to a virtual line (z = z f ) where cathode boundary conditions are imposed. Part of the plume is included in the spatial computational domain of the model in order to provide some information also outside of the channel.  In the following, we concisely introduce the model equations dividing the discussion for each species included in the model. Since we are dealing with a 1D model, it is implicit that all the equations and quantities involved are averaged over z-sections.

Neutrals
Neutral atoms are assumed to have constant axial velocity u n = u n e z and constant temperature. Thus, only the continuity equation is included in the model as a classical linear transport/reaction equation with u n > 0 constant and assigned among the problem input data: ∂n n ∂t + ∂ ∂z (u n n n ) = −n e n n k I +ṅ w , where n n and n e are the number density of neutrals and electrons, respectively. As will be detailed later, plasma is assumed to be quasi-neutral, so that n e = n i where n i is the ion number density. The term k I in Equation (1) is the ionization-rate coefficient of the neutral propellant atoms due to electron impact, and the product n n k I is the corresponding ionization frequency. The ionization-rate coefficient k I is a function of the electron internal energy ε. More specifically, we use a tabulated version of k I using collision cross section data from the lxcat database [56] and processed with the Bolsig+ solver [57] in order to obtain rates for a Maxwellian electron energy distribution function as a function of ε.
In Equation (1), the termṅ w represents the neutral particles flowing in the domain from the lateral walls due to ion recombination. The model used to describe the plasma interaction with the channel walls and the expression forṅ w is presented in detail in Appendix A. Note that this term, as well as the other wall sink terms introduced in the electron momentum and energy equations (see Section 2.3), is modulated through a wall interaction coefficient, α, that accounts for the uncertainties in the description of the plasmawall interaction and is one of the calibration parameters that are used to fit the experimental reference data, as detailed later in Section 2.5.

Ions
Ions are assumed to be isothermal with constant temperature T i , which is assigned among the problem data. For this reason, only two scalar equations are included for the ion dynamics in the present 1D model, i.e., the continuity equation and momentum balance in the z-direction. Indicating with u i the ion velocity in the axial direction, the continuity equation writes as: ∂n i ∂t + ∂ ∂z (u i n i ) = n e n n k I −ṅ w .
Note that the decrease in ion density outside the thruster channel due to the plume expansion (see, for instance, [29]) is not taken into account here as it is not considered important in the dynamics of the low-frequency breathing mode, which is the main objective of the present analysis.
For the momentum balance in the z-direction, we assume unmagnetized, collisionless ions, and the resulting equation writes as follows: where p i is the ion pressure, m i is the ion mass and Φ is the electric potential (so that E z = − ∂Φ ∂z is the component of the electric field in the z-direction). We further specify that the ion pressure is related to the ion temperature, assuming an isotropic Maxwellian distribution for the random thermal ion velocity, thus leading to the ideal gas law: where k B is the Boltzmann constant. Note that in the momentum balance (3) we have neglected momentum transfer due to charge exchange collisions, which are sometimes included in similar models in the literature (see, for instance, [29]). While the formulation presented above allows for a finite, albeit constant, ion temperature, the energy of the ion fluid is largely dominated by the kinetic component, and thermal effects are not considered to have a significant impact on the dynamics of breathing mode oscillations. Therefore, in the simulation results presented in Section 4, the ions are assumed to be cold by setting T i = 0.

Electrons
The continuity equation for the electrons is identical to that of ions as a consequence of the quasi-neutral assumption, which implies the point-wise and instantaneous equality n i = n e : ∂n e ∂t + ∂ ∂z (u e n e ) = n e n n k I −ṅ w , where u e is the electron velocity component in the axial direction. Since the characteristic length and time scales of the phenomena under investigation are much larger than the Debye length and of the inverse of the plasma frequency, and considering that the wall sheaths are excluded from the domain, the quasi-neutral assumption is well justified. By subtracting the two conservation equations for ions and electrons, namely Equations (2) and (5), and multiplying by the electron charge e, the current conservation equation is obtained: ∂ ∂z (e n i u i − e n e u e ) = 0.
Equation (6) essentially states that the current density, defined as is a funcion of the sole time, i.e., J = J(t).
As concerns the momentum balance of the electrons, we use the drift-diffusion approximation, which implies that electrons are assumed to be at steady state within the characteristic evolutive time of ions, i.e., inertia terms are negligible in the momentum balance for electrons. This last assumption is consequent to the small ratio between electron and ion masses, i.e., m e /m i 1, which implies a low-Mach number approximation for electrons, i.e., electrons move at a much lower velocity with respect to their thermal speed, and this is shown in detail in the asymptotic analysis proposed in [58]. Finally, electrons are supposed to be magnetized, and the magnetic field is assumed to be purely radial, i.e., B = B r e r .
As a result, the momentum balance in the axial (z) and azimuthal (θ) directions becomes: m e n e ν e u e = − ∂p e ∂z + e n e ∂Φ ∂z + e n e u eθ B r , m e n e ν e u eθ = −e n e u e B r , where u eθ is the electron velocity component in the azimuthal direction and the coefficient ν e is the momentum transfer collision frequency. In the proposed model, ν e is given as the sum of several concomitant contributions. In particular: where ν ew is the electron-wall collision frequency, and ν c is a generic collision frequency taking into account all the collisions that electrons undergo with neutral atoms, including elastic, ionization and excitation collisions, whose values result from the product of the neutral density and the relevant reaction rate k x , which were obtained using Bolsig+ [57] and cross section data extracted from the lxcat database [56]. The wall collision frequency ν ew , which accounts for the electron momentum loss to the channel lateral walls, follows from the plasma-wall interaction model detailed in Appendix A, and it is modulated with the wall interaction coefficient α, as previously discussed. The anomalous collision frequency ν a is finally modeled to be proportional to the local cyclotron frequency of the electrons ω e = (e B r )/m e using the classical expression: where β is a free non-dimensional parameter among those included in the calibration (see Section 2.5).
The momentum balance in the azimuthal direction, Equation (9), can be solved in terms of u eθ and the result substituted in the axial balance (8) so as to obtain a single equation in u e : n e u e = −µ n e 1 e n e ∂p e ∂z The cross-field mobility µ in the previous equation is defined as: where µ 0 = e/(ν e m e ) is the unmagnetized mobility and Ω = ω e /ν e is the Hall parameter.
If we further assume that electron velocities follow an isotropic Maxwellian distribution of temperature T e , the electron pressure p e can be related to T e as follows: Thus, Equation (12) becomes: Unlike ions, the electron temperature is not constant and is among the unknowns of the problem. Thus, it is necessary to include the energy equation for electrons in the model. In particular, following Reference [29], we consider here the internal energy equation for the electrons, obtained by subtracting the kinetic energy equation from the equation for the total electron energy. As a result, the contribution of the electron potential Φ disappears in the internal energy balance. Using the relation between internal energy and electron temperature for a Maxwellian distribution (ε = 3 2 k B T e ), the conservation of internal energy writes as follows: In this equation, the term D is a diffusion term, the term H represents the conductive heat flow, S coll represents ionization and excitation losses, and S w (which is also modulated with the wall interaction coefficient α) models the energy losses at the walls and is discussed in detail in Appendix A. More specifically, the collisional energy loss S coll = −n n n e K is expressed in terms of the collisional energy loss coefficient K, which is a function of the electron energy. In the same way as the ionization rate, K(ε) values were extracted from the Biagi database [56] and were calculated as a function of the electron temperature for a Maxwellian energy distribution using Bolsig+ [57].

Model Formulation and Numerical Discretization
The 1D model is based on the equations introduced in the previous sections, which are further manipulated as detailed here. The first manipulation carried out concerns the ion axial momentum Equation (3). In particular, following the literature [28], the electric potential in the equation is eliminated by means of Equation (15), which gives: This simplification enhances the numerical stability in the solution of the hyperbolic (inviscid) part of the ion equations, as amply demonstrated in the literature. The main reason is that the correct acoustic speed to consider in the hyperbolic part of the ion equations is the one of Equation (17), i.e., the one related to the following equivalent sound speed c i for ions: A second manipulation of the equations consists in combining the expression for the current density, Equation (7), and the electron momentum Equation (12) and integrating over the entire domain to compute the discharge current density: where we assumed n = n i = n e and where ∆V is the potential difference between the first point of the investigated quasi-neutral domain (z = 0) and the cathode (z = z f ). Specifically, if V d is the total discharge voltage applied between anode and cathode, ∆V takes into account the potential drop (φ a ≥ 0) that occurs in front of the anode sheath: where j e,a = J − j i,a is the electron current density at the anode, j i,a is the ion current density at the anode, and j e,th = en k B T e 2πm e is the current density associated with the electron thermal flux, assuming a Maxwellian distribution function at the anode sheath edge. The expression for the sheath potential drop of Equation (21) follows a classical description of a charged electrode immersed in a plasma and is limited to zero since the sheath found in front of the anodes of Hall thrusters is considered to be electron repelling [29].
Once J is known (and thus the total discharge current I by multiplication with the channel cross section area), the electron velocity can be deduced by charge conservation: Finally, the electric potential can be computed, if needed, directly by using the electron momentum Equation (15), as usual when working with the quasi-neutral plasma approximation.
The final system of equations is strongly coupled and non-linear. However, when semidiscretized in time following a line-method approach, i.e., thinking about a discrete-in-time advancement of the model, the following fractional step method is proposed to advance the equations from a generic time level t k to t k+1 = t k + ∆t. Provided that all plasma quantities are known at time t k , sub-problems are identified and solved in consecutive order: 1.
The neutral continuity equation is solved providing n k+1 n ; 2.
The ions continuity (2) and momentum (17)  The electron temperature T k+1 e is computed solving Equation (16) and considering the electron velocity u e at time level k; 4.
The current density J k+1 is computed using Equation (19);

5.
The electron velocity u k+1 e is computed, solving the charge conservation Equation (22).
The first three subproblems above are governed by PDEs, and appropriate boundary conditions must be thus supplied, while boundary conditions on the potential are already included in the integral equation at point 4, as discussed above.
Subproblem 1 is a scalar transport equation with constant convection velocity u n > 0. Consequently, only one boundary condition on n n is required, which is the neutral density specified at the anode, i.e., at z = 0. This boundary condition is easily recovered by knowing the thruster channel cross section area (A) and the total injected mass flow rate (ṁ) at the operating point under investigation; considering also the ion recombination at the anode, the condition reads as: Subproblem 2 is a hyperbolic system of two equations, thus a variable number of boundary conditions is needed at each of the two boundaries of the domain depending on the local characteristic velocities. In the simulations carried out here, the solutions are always supersonic at the outlet boundary (i.e., z = z f ), so that no boundary conditions are required on that point. Conversely, at the anode, the local characteristic velocities are such that there can be characteristic lines entering the domain. In particular, in order for the electron-repelling sheath to form at the anode, it is imposed that the ions must enter at the anode sheath edge with a velocity higher than or equal to the sound speed c i . Subproblem 3 is an elliptic 1D problem, which thus requires two boundary conditions on T e , one per domain boundary. At the anode, the electron energy flux is specified through a classical analysis of an electron-repelling sheath in front of a biased electrode, and is set equal to j e,a e (2 k B T e + e φ a ), where the values of the plasma parameters at the sheath edge are approximated with those at the center of the first cell adjacent to the anode. At the cathode boundary, instead, a constant temperature of 2 eV is imposed.
Concerning the numerical solution of the model's PDEs listed above, all of them are discretized in space using a finite-volume formulation. In particular, the neutral dynamics (point 1) are solved by a first-order upwind discretization. Ion dynamics (point 2), which is governed by a hyperbolic system of two equations/unknowns with a flux function satisfying the homogeneity property (see [59]), is solved by the Steger-Warming [60] flux vector splitting method. The same splitting method is used to impose characteristics-based boundary conditions on the system. The electron energy equation (point 3) is solved by a second-order finite-volume formulation with a central evaluation of gradients at the cell interfaces. Finally, the integral equation for the evaluation of current density (point 4) is carried out consistently with the finite-volume formulations adopted for the PDEs. Regarding time discretization, since the splitting of the whole coupled problem into the five consecutive subproblems described above is a first-order splitting in time, all subproblems are advanced in time by a first order Euler implicit/explicit scheme. The implicit scheme is used for the electron energy equation, while the remaining equations are advanced in time explicitly.
A six-core/12-thread i7-8750H CPU PC with 16 Gb of RAM and a solid state drive was used for the computation. Adopting the numerical scheme described above, the simulation time needed to compute 1 ms of the plasma dynamic was about 1310 s. Note that the code is not optimized, and there are ample margins of improvement in this respect.

Model Calibration
The model described in previous sections includes a number of calibration coefficients. These account for the unknowns in the physical description and allow to accommodate the simplifications introduced in the model. In detail, the three elements that are considered as calibration coefficients are: (i) the neutral velocity, u n , (ii) the wall interaction coefficient, α, and (iii) the anomalous diffusion coefficient, β, that can, in general, be a function of axial location, i.e., β = β(z).
The uncertainties on the value of the neutral velocity u n come from the absence of a detailed simulation of the flow injection and of the anode thermal behavior coupled with the assumption of a constant velocity across the domain. For what concerns the wall interaction coefficient α, a physical closure is avoided and the coefficient is retained in the model to account for the unresolved radial gradients in the plasma profiles and to counteract the uncertainties in the plasma-wall semiempirical descriptions (see Appendix A). Moreover, the coefficient allows to decrease the plasma-wall interaction, simulating partial shielding of the channel walls from the plasma due to the local inclination of the magnetic field, e.g., eroded and end-of-life conditions, such as the case under investigation in the present study. Finally, β represents the impact on the electron mobility of plasma turbulence and azimuthal oscillatory modes [1]. Following experimental and numerical investigations found in the literature [28,29,32,61], β was assumed to have two different values, inside and outside of the channel, with the outside anomalous coefficient significantly higher (by a factor of 100) than the inside one and with a smooth transition of variable steepness between the two.
The three parameters were progressively varied in order to allow the simulation results to match as closely as possible the experimental reference data. The overall discharge current signal was taken as the reference calibration quantity for the model. More in detail, while varying the calibration coefficients, when an oscillatory limit-cycle solution was found, the main parameters that were monitored and compared with experiments were: (i) the average, minimum and maximum value of the discharge current during the breathing mode cycle, (ii) the root mean square value of the oscillatory component of the discharge current signal and (iii) the dominant frequency component in the power spectral density of the discharge current signal.

Test Item, Facility and Global Diagnostics
The device under investigation is the second development model (DM2) of SITAEL's HT5k, a 5 kW-class Hall thruster. This thruster has been designed with a configurable soft alloy magnetic circuit, so as to generate fields with topologies spanning from conventional to magnetically shielded ones [38,62]. The present work deals only with a conventional magnetic topology, called M1 (see [38] for more details), with almost radial magnetic field lines at the channel exit. Nevertheless, it is worth noting that this configuration has a chamfered ceramic channel that reproduces an end-of-life condition. The thruster is coupled to an externally mounted HC20 high-current hollow cathode [63]. The thruster has a nominal discharge power of 4.5 kW. However, the current analysis focuses on a case in which the thruster operates at nearly 2.6 kW, with a discharge voltage of 300 V and a mass flow rate of 8.4 mg/s, as this condition is characterized by stronger breathing mode oscillations. Table 1 summarizes the operation parameters of the condition under study. The test campaign was carried out in SITAEL's IV4, a 2 m diameter and 4.2 m long vacuum facility [64]. This chamber is a nonmagnetic stainless-steel vessel equipped with an oil-free pumping system capable of a pumping speed of 70,000 L/s for xenon, thanks to a combination of cold heads, turbomolecular pumps, and cryogenic pumps. The system typically achieves a base pressure below 10 −7 mbar (<7.5 × 10 −8 Torr), and is capable of maintaining the pressure below 10 −4 mbar (<7.5 × 10 −5 Torr) during thruster operation. Additionally, the discharge current measurements were made using a LEM LA25-NP current probe and a Tektronix DPO 4104 digital oscilloscope with a sampling rate of 5 MHz.

Fast Diving Triple Langmuir Probe
To have a better insight into the behavior of the breathing mode and to assess the local effects on the plasma properties, the experimental setup included a fast diving triple Langmuir probe. The use of a triple probe, despite reducing the space resolution in the perpendicular plane, allows achieving a good temporal resolution without the need to perform rapid voltage sweeping, which introduces complications in the results interpretation due to the presence of important capacitive effects, as is the case with single Langmuir probes used in previous studies (see, for instance, [51,54]). This diagnostic is described in detail in previous publications from the authors [39,40,55], therefore, only a summarized discussion is included in the following.
Triple Langmuir probes are plasma diagnostics that consist of three conductive electrodes mounted on the tip of an isolating ceramic bar. Using a proper electric arrangement of the electrodes, this type of probe allows to obtain instantaneous measurements of the plasma density n, plasma-to-ground potential V gp and electron temperature T e , [55,65]. The three electrodes are typically denominated Bias (B), Common (C) and Float (F). A potential difference is applied between C and B while the F electrode is left floating. For the case under analysis, a bias of 30 V is applied between electrodes B and C. This value was selected to provide good accuracy of the measurements (in particular, for the electron temperature range of interest) and, at the same time, to comply with safety requirements and to avoid test setup complications. Triple Langmuir probe measurements consist of the simultaneous acquisition of the current flowing between B and C (I BC ), the voltage between F and C (V FC ) and the voltage between ground and C (V gC ). The three acquired values, coupled with a model of the plasma interaction with the electrodes, are sufficient to reconstruct the local values of plasma density, potential and electron temperature.
The three probe electrodes are made of a 75% tungsten 25% rhenium alloy, with a diameter of 0.178 mm and a length of 2.6 ± 0.1 mm. The electrodes are installed inside a 1/8 alumina insulating tube. Although it is true that the high secondary electron emission of alumina would result in a higher perturbative effect in the high-temperature plasma regions, its high mechanical strength and availability in the required format led to the selection of this ceramic for the probe body. A minimum distance of 2 mm between the electrodes was established to avoid undesired interaction effects. The Debye length of the probed plasma is always smaller than 0.1 mm, so the separation distance of 2 mm allows to always have at least 20 Debye lengths of separation between the electrode tips, which can be considered as the best trade-off between minimizing electrodes interaction and maximizing the probe spatial resolution. Additionally, the selected electrode length of 2.6 mm allows neglecting tip effects.
The probe was mounted on a mechanical arm, shown in Figure 2, that allowed for its rapid insertion and extraction from the plasma region. During its motion, the probe performed a circular trajectory with a radius of 350 mm. Data acquisition was only carried out in the final 0.27 rad arc of the probe motion in the near-plume and channel region of the thruster. Taking into account that the arm is 350 mm long, this implies that the probe tip had a maximum deviation from the channel centerline of 0.2 mm inside the channel and 9.5 mm in the plume. In order to minimize the probe exposure to the harmful high plasma density and temperature environment, the arm was kept in a parking position when non-operational. To perform the measurement, the arm was moved in and out of the acquisition region in 200 ms by a high-speed magnetic actuator. While moving, an encoder recorded the position of the probe with a resolution of 0.3 mm.
The adoption of a triple probe allows gathering instantaneous measurements of plasma properties without the need of voltage sweeps, at the cost of assuming that the electrons present in the plasma have a non-drifting Maxwellian velocity distribution function (EVDF). As explained by the authors in References [39,55], the assumptions on the EVDF are acceptable to a first order, particularly in the plume, where the electron drift velocity is negligible with respect to their thermal velocity. Moreover, a set of correction parameters, introduced by Saravia et al. in Reference [39], allows mitigating the impact of the electron drift on the measurements. It is important to note that the characteristic time of the probe motion (O(1) Hz) is orders of magnitude longer than that of the oscillatory phenomena under study (O(10 4 ) Hz), making it possible to assume that the probe is steady during several plasma oscillations, allowing to make considerations on the local unsteady behavior of the plasma even if the probe is moving. The probe signals were conditioned using an analog electronic box (see [55] for more details) and acquired with the same oscilloscope used to acquire the discharge current signal, i.e., a Tektronix DPO 4104 digital oscilloscope, set at a sampling rate of 5 MHz. The data storage was triggered using a photocell activated by the passage of the arm. The signal conditioning box was based on the AD215 120 kHz bandwidth, low distortion, isolation amplifiers, which limited the effective maximum frequency of the system. This bandwidth is much higher than the frequency expected from the breathing mode and its first harmonics, so that the information present in the gathered data is sufficient to properly measure the phenomena under study.

Time-Resolved Bayesian Reconstruction of Plasma Parameters
The reconstruction of the plasma properties is performed using a Bayesian inference methodology, described in detail in [55], which yields probability distributions for the plasma parameters and, as a consequence, permits calculating the most likely value and the corresponding uncertainty. Bayesian data analysis is a branch of Bayesian probability theory that permits a robust fusion of data coming from different sources while consistently keeping track of uncertainties. Within this framework, unknown parameters are treated as probability distributions and experimental data sets are combined with physical models and prior knowledge to infer the values of the investigated physical parameters. The method makes use of the Bayes theorem to update the probability distributions of the parameters by considering the likelihood of experimental measurements in the light of a physical model, and relies on the sum rule to marginalize the distributions and analyze each variable separately [66,67].
The process applied in the present investigation uses a parametric solution of the Lafambroise model of the particle collection in cylindrical probes [68], a Gaussian likelihood to uncertainties (model, measurements, etc.) and an implementation of the nested sampling algorithm to explore the parameter domain and calculate the parameter probability distributions [66,69].
Given that the electrons have a significant drift velocity in the regions close to the channel exit, the assumption of a non-drifting Maxwellian population is not completely fulfilled. Hence, a correction is introduced considering changes in the floating potential observed by the probe electrodes caused by the presence of the electron drift velocity. The value of the correction is calculated as a function of the position, using data obtained with a null bias applied between the B and C electrodes, and acquired with different probe geometric arrangements, so as to take into account the mutual screening effects between the different electrodes, as detailed in References [39,55]. The correction factors are then introduced as prior information in the analysis of unsteady data.
To perform the analysis, the raw data time series (V gC , V FC and I BC ) sampled at 5 MHz are downsampled, calculating the average and the covariance matrix of series of 20 samples. This process reduces the frequency of the new data set to 250 kHz, as shown in Figure 3. The resulting time series are then introduced in the Bayesian inference process together with the prior distributions of the sought parameters and the correction factors. This process leads to the estimation of the joint probability distributions and then marginalizes each parameter, yielding the time series of n, T e and V gp , as illustrated in Figure 4.

Results and Discussion
The detailed analysis of the experimental measurements is beyond the scope of this discussion; the reader is invited to read the work of Giannetti et al. [55] for additional details on the probe data. In the following, we focus on the simulation results obtained after model calibration, and on the comparison with the reference experimental data, highlighting the main physical phenomena that the model is capable of reproducing, as well as the main shortcomings.
In order to simulate the HT5k DM2 (M1) operation, the thruster salient features were imported into the model. These included the core geometrical properties of the channel, as well as the radial magnetic field profile along the channel centerline. Moreover, the xenon operating condition detailed in Table 1 was imported as input data for what concerns the injected mass flow rate and the applied discharge voltage. Finally, the simulation spatial domain was set to extend from the anode surface up to a full channel's length downstream of the exit section, i.e., z f = 2L ch , where cathode boundary conditions were applied. With the domain and all required inputs fixed, the model calibration procedure detailed in Section 2.5 was followed to find the best combination of the three parameters α, β and u n .
The best fit of the discharge current experimental data was found for u n = 395 m/s, α = 0.115, and β equal to 0.075 inside the channel and equal to 7.5 in the plume (keeping a ratio of 100 between the two). Figure 5 depicts a direct comparison between the discharge current measured by the oscilloscope during thruster testing and the discharge current predicted by the simulation of the operating point under investigation. Moreover, Table 2 reports the quantitative comparison of the main parameters of the current signal that were monitored during calibration in order to obtain the best match between the simulations and the experiment.  The superposition of all main characteristics of the current signal is remarkable, with all control parameters exhibiting a relatively small difference between simulations and the reference experiment. This is strong evidence that the proposed numerical model incorporates all the necessary physics to correctly reproduce the core mechanism of breathing mode oscillations.
The neutral speed is typically linked to the anode temperature since the flow is considered sonic at injection. Following the accommodation coefficients discussed in [43], typical anode temperatures found in Hall thrusters translate to typical neutral injection speeds in the 100 to 300 m/s range. Therefore, the value of u n required to match the experimental results is relatively high. During the calibration process, it was found that the value of the neutral velocity has a direct impact on the dominant frequency of the discharge current oscillations when present. This is coherent with an intuitive understanding of breathing mode as an ionization instability: if the neutral speed increases, the time needed by the neutral flow to reach the high temperature region decreases and, thus, the plasma surges become more frequent. Since the experimental data presented a relatively high dominant frequency (see Table 2), a correspondingly high neutral speed was needed. Relaxing the assumption of constant neutral velocity over the domain, and thus solving the neutrals momentum equation, could result in calibrated injection speeds more in line with expectations.
For what concerns the wall interaction correction factor α, a significant reduction of the particle flow to the walls, when compared to that predicted by classical models, was needed to correctly reproduce the experimental results, and especially to match the total oscillation amplitude. This trend was expected, in part, since the radial gradients of the plasma properties are not resolved in the current model, and the plasma density is bound to decrease towards the walls due to the ion acceleration in the pre-sheath region. Moreover, a significant uncertainty exists in the representativeness of fluid sheath models, such as the one described in Appendix A, which, according to multiple studies, tend to overestimate the wall losses [70,71]. A further justification lies in the specific configuration of the thruster model under investigation. As described in Section 3.1, the M1 version of the HT5k DM2 thruster was designed to mimic end-of-life conditions of the Hall thruster, with a chamfer at the end of the ceramic channel that diminishes the local incidence angle of the magnetic field lines on the channel walls. As demonstrated in the literature [38,72,73], this can effectively produce a shielding effect of the channel walls from the plasma. It is not possible to reproduce this feature in a 1D model in which magnetic field lines are assumed to be radial and orthogonal to the channel walls, but this is reflected in the reduction of the wall interaction coefficient α.
Finally, it was found that the anomalous collisionality parameter β has a strong impact on the dynamic behavior of the solution. Particularly, slight variation in the anomalous diffusion parameter, especially inside the channel, implied significant differences in both the average current value and the amplitude of the oscillations. A similar trend was also observed by Hara et al. in [29]. Interestingly, also the steepness of the transition in the anomalous diffusion coefficient profile from the inside value to the outside one was found to have a significant influence on the oscillation dynamics. The resulting profile of the anomalous collision frequency is presented in Figure 6. It is worth noting how this profile recovers values of the anomalous collisionality close to the ones obtained by other authors and calibrated on experimental data [29,46,61].
With the model now calibrated on the experimental profile of the discharge current, it is possible to investigate the dynamics of all intensive plasma properties predicted by the simulations in the channel and near-plume. Figures 7 and 8 depict the evolution of all unknowns of the system over the entire z − t plane for a simulated time of 0.2 ms. Figure 7 details the dynamics of the heavy species (ions and neutrals), while Figure 8 reports the behavior of the properties related to the electron equations.
The results give a striking picture of the effects of breathing mode oscillations on the plasma. Neutral atoms are generated at the anode boundary through the mass flow injection of propellant. They move through the channel (see diagonal lines in Figure 7 (middle)) until they reach a region with a sufficiently high electron temperature and density to achieve effective ionization. At this point, a surge of plasma density occurs in the channel, indicated by the periodic density peaks of Figure 7 (left). This dynamic exchange between neutrals and plasma is very clear in the near anode region, where the plasma density and neutral density oscillations are almost in phase opposition. The newly generated highdensity plasma is then moved by the electric field both downstream, exiting the thruster at high speed, and upstream towards the anode, where a weak reversed electric field is established to counteract the local electron pressure and ensure the anode sheath boundary conditions are respected. This ion backflow recombines at the anode, increasing the local neutral density and further fueling the breathing mode feedback mechanism, before the cycle starts again. Observing that the oscillations are found for all properties over the whole plasma domain with the same characteristic frequency validates the interpretation of breathing mode as a global self-sustained instability. The same picture was deduced from the experimental results of [55] where a wavelet analysis of the measurements was performed over the probed region.  The measurements performed and the data processing technique presented in [55] allowed the reconstruction of local oscillations of plasma density, potential and electron temperature in a region spanning from 15 mm upstream of the thruster exit plane up to 70 mm downstream. Although the domain investigated with the simulations was limited to just one channel length downstream of the channel exit section, it is interesting to compare the local measurements of the plasma properties with the results of the simulation in order to understand the degree of accuracy achieved by the adopted formulation in the reconstruction of the local plasma behavior after calibration on the sole discharge current signal. This comparison is presented in Figure 9, where the experimental oscillations are reported only for the region where the probe was not perturbing the plasma flow (z > L ch + 2 mm). It should be noted that the plasma potential measurements refer to the voltage difference between plasma and ground, while the simulations are relative to the cathode plasma potential, therefore, the comparison is only quantitatively valid neglecting the cathode reference potential. The order of magnitude of all relevant properties and the general trend of the profiles and oscillations over the probed domain is successfully reconstructed by the simulations: the plasma density average value and oscillation amplitude is higher closer to the anode and decreases moving away from the thruster exit plane, the potential drop is mostly concentrated in the high magnetic field region of the thruster, close to the channel exit. Consequently, the highest value of the electron temperature is found in the same region. The numerical model also manages to simulate the oscillation of the acceleration region during the breathing cycle, with the fraction of potential drop happening outside of the channel varying during the oscillation period.
Two major differences are found between the measurements and the simulation results. First, the plasma density average value and oscillation amplitude seem to be underestimated by the code in the near-plume region. The numerical results predict the peak of plasma density inside the channel, further upstream with respect to the measurements. Second, the predicted temperature profile, and especially the peak value, is higher than the measurements, even though the peak is in a region where the measurements may be in part affected by the perturbations induced by the probe. The adopted formulation appears to overestimate the electron temperature and to move the acceleration region more inside the channel with respect to the empirical evidence. For what concerns the density profile, the neutral back ingestion from the chamber and from the cathode during testing could supply additional slow particles in the near-plume region. In general, the discrepancies observed between numerical results, and experimental data are also tied to the core assumptions of the simplified model. First, the code presents values averaged over z-sections of the channel, while the measurements refer to the local values on the channel centerline. Moreover, the plume expansion is not modeled in detail in the simulations, leading to potential local inaccuracies in the representation of the plasma dynamics downstream of the channel exit. Considering the strong sensitivity of the local plasma parameters to the calibration coefficients, both discrepancies are also linked with the specific profile of the anomalous collisionality that we have selected and with the value of u n and α needed to match the current profile, but are also a consequence (at least partially) of the application of a 1D axial formulation to an end-of-life configuration, as previously discussed.
It is worth mentioning that, in general, the plasma turbulence effect on electron transport is a function of the local plasma properties (see, for example, Reference [74]) that vary during the breathing mode cycle. Therefore, β could also be a function of time, i.e., β = β(z, t). While we have decided to fix the axial profile of the anomalous collisionality to limit the degrees of freedom in the model calibration, experimental investigations in the literature [46] have highlighted how β varies significantly during the breathing mode cycle.
The adoption of a time-varying profile for the anomalous diffusion coefficient could allow for a closer representation of the experimental results within the boundaries of the presented formulation and will be the subject of future investigations.

Conclusions
In this paper, we have presented a synergic experimental/numerical investigation of the breathing mode in a 5 kW-class Hall thruster. In particular, we have shown the potential of using an informed 1D fully-fluid model of the plasma flow as a complementary tool with respect to the experiments. On the basis of available measurements, the current signal in particular, the model is aimed at estimating evolving plasma properties that are difficult to measure. Thanks to the availability of a detailed experimental database, which has been collected using a diagnostic technique based on the use of a fast-diving triple Langmuir probe, we have provided here an assessment of the prediction capabilities of the calibrated model.
More specifically, the model was set up to simulate one operating point of the HT5k and calibrated to reproduce the measured discharge current profile. Three calibration parameters implemented in the numerical description allowed to properly tune the model to the experimental observations: (i) the neutral injection velocity, (ii) the wall interaction coefficient and (iii) the anomalous diffusion profile. Following calibration, the discharge current profile predicted by the simulation closely matched the experimental results, exhibiting relatively small differences on all control metrics. This promising result confirms the capabilities of the adopted numerical description to reproduce the core mechanism of breathing mode oscillations.
The time evolution of the main plasma properties predicted by the calibrated simulations over the entire domain was presented and discussed, highlighting the dynamics of the longitudinal oscillation and its underlying global and self-sustained nature.
Finally, the measurements obtained for the plasma density, potential and electron temperature oscillations in the near-plume were compared to the numerical results. The model successfully managed to recreate the order of magnitude and overall trend of the properties' profiles. Nevertheless, the simulations underestimated the values of density, overestimated the electron temperature, and predicted the acceleration region to be further upstream with respect to the experimental results.
Overall, the adopted approach of informing a reduced-order numerical description of the plasma dynamics with high-frequency measurements of the discharge current has demonstrated to be a promising path to gather additional insight on the dynamics of oscillatory modes in Hall thrusters. The possibility of comparing the simulation results with measurements of the local plasma properties represents a unique framework for the identification of the main physical processes that are needed to model the discharge of Hall thrusters. Future investigations will explore the application of the described approach to other thruster configurations and operating conditions, focusing on the adoption of a time varying profile of the anomalous diffusion and on the inclusion of the effects of the magnetic field topology in the model. Since experimental measurements of the discharge current are often easily available from the functional characterization of a Hall thruster, the proposed combined numerical-experimental approach may represent an important step toward the understanding of the characteristic physical processes of Hall thruster discharges. Funding: The work described in this paper has been funded by the European Space Agency in the framework of the contract 4000113279/14/NL/KML "Low-Erosion Long-Life Hall-Effect Thruster" and by the European Union under the H2020 Programme ASPIRE-GA 101004366. The views expressed herein can in no way be taken to reflect the official opinion of the European Space Agency.

Data Availability Statement:
The data presented in this study are available in Reference [55].

Acknowledgments:
The authors wish to express their gratitude to Ugo Cesari, Nicola Giusti, Luca Pieri, Stefano Caneschi and Carlo Tellini for their valuable assistance in preparing and performing the experimental campaign. Fruitful discussion with Fabrizio Paganucci and Francesco Califano are gratefully acknowledged.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Plasma-Wall Interactions
The expressions for the ion recombination at the walls, the electron-wall collision frequency, and the electron energy losses at the walls are determined based on the formulation of Hobbs and Wesson for a 1D sheath with cold ions and in the presence of secondary electron emission [75]. Those quantities are then modulated using the wall interaction coefficient α, which belongs to the set of parameters used for the calibration (see Section 2.5) and are included in the model for z < L ch .
The particle balance at the dielectric walls of the discharge channel can be expressed as: where Γ iw and Γ ew are, respectively, the ions and primary electrons fluxes to the dielectric walls, and σ is the effective secondary electron emission yield for the wall material considered (BN-SiO2 in this case). Different empirical formulas for the secondary emission yield at low plasma temperatures (up to 50 eV) can be found in the literature [2,21,30,31]. In this work, we adopted the expression recommended in Reference [2]. Also accounting for the possibility of sheath space charge saturation (that for xenon occurs when σ = 0.983), this expression writes as: σ = min 1, 36 · 0.123 T 0.528 eV , 0.983 .
Neglecting the minor modification to the Bohm velocity in the presence of secondary electron emission, discussed in Reference [75], the particle flux to the walls in a control volume delimited by two channel sections at a distance dz is: where R 1 and R 2 are the inner and outer channel radii and where we have introduced the wall interaction calibration coefficient α. The frequency of ions impacting the walls (ν iw ) can be defined as the flux of particles to the wall per unit density and per unit volume and, therefore, writes as: Consequently, the source/sink term representing the particle-wall interaction in the continuity equations can be expressed as: Note that the same expression applies to the electron population due to the zero net current condition to the walls and the quasineutrality assumption. Using Equation (A1), the equivalent electron collision frequency accounting for the momentum losses of the electron fluid to the walls can be defined as: Concerning the electron energy losses, assuming a Maxwellian electron velocity distribution function at the sheath edge, the mean energy of primary electrons lost at the walls is 2k B T e + eφ w , where the sheath potential drop is given by: Therefore, assuming that the secondary electrons thermalize with the main electron population, and neglecting the energy with which they are emitted from the wall, the term accounting for the electron power loss to the walls in the electron energy equation writes as: