Collisionless electron cooling in a plasma thruster plume: experimental validation of a kinetic model

A central challenge in the modeling of the near-collisionless expansion of a plasma thruster plume into vacuum is the inadequacy of traditional fluid closure relations for the electron species, such as isothermal or adiabatic laws, because the electron response in the plume is essentially kinetic and global. This work presents the validation of the kinetic plasma plume model presented in (Merino et al 2018 Plasma Sources Sci. Technol. 27 035013) against the experimental plume measurements of a SPT-100-ML Hall thruster running on xenon presented in (Giono et al 2018 Plasma Sources Sci. Technol. 27 015006). The model predictions are compared against the experimentally-determined axial profiles of electric potential, electron density, and electron temperature, and the radial electric potential profile, for 6 different test cases, in the far expansion region between 0.5 and 1.5 m away from the thruster exit. The model shows good agreement with the experimental data and the error is within the experimental uncertainty. The extrapolation of the model to the thruster exit plane and far downstream is consistent with the expected trends with varying discharge voltage and mass flow rate. The lumped-model value of the polytropic cooling exponent γ is similar for all test cases and varies in the range 1.26–1.31.


Introduction
The interaction of the plasma thruster plume with the surrounding elements of the spacecraft poses a major concern for the integration of electric propulsion systems like gridded ion engines (GIT) and Hall thrusters (HTs) in space missions [1][2][3]. The main threat are the charge-exchange ions generated in the plume near-region, which can impinge on neighboring surfaces and interact mechanically, electrically, or chemically with them. These ions originate from collisions between primary ions and low-velocity neutrals outside of the thruster, and their motion is dominated by the local electric field. For this reason, determining accurately the electric potential in the plume near-region is a key step toward the predictive simulation of charge-exchange ion effects.
The electric potential map ( ) f x arises from the expansion of the primary ions and the neutralizing electrons, but depends especially on the electron thermodynamics, i.e. on the evolution of the electron temperature T e in the plume, with hotter electrons leading to larger electric fields ( f D~T e e ). Due to computational constraints, the majority of existing simulation codes follow a fluid description for the electrons, using a simple electron cooling law in the form µ g-T n e e 1 as the closure relation, where n e is the electron density and γ is a polytropic exponent. The limit γ=1 means an isothermal electron expansion and results in Boltzmann's relation. Unfortunately, in the near-collisionless regime of plasma thruster plumes, the fluid approach is unjustified, and a fully kinetic model must be used to obtain physically-correct solutions. The inadequacy of the fluid approach is evidenced, in particular, in the case of an isothermal electron species, which leads to an infinite potential fall as the electron density drops to zero. Moreover, even accepting the fluid approximation and a polytropic closure with γ>1, the parameter γ must be regarded as an extra degree of freedom of the fluid model, which has to be determined for each thruster and each operating condition empirically.
In order to address these shortcomings, a kinetic electron model for plasma plumes expanding into vacuum [4] (and the open-source AKILES code [5]) was recently developed to compute self-consistently the near-collisionless electron expansion. The model is valid in the far expansion region, once local thruster electric and magnetic fields, neutralizer effects, and collisions have become negligible. The results of this model show the gradual electron cooling and development of anisotropy in the plume, and can be used to inform existing plume/spacecraft interaction codes like EP2PLUS [6] or SPIS [7], providing more accurate electron fluid closure relations than the currently-used isothermal and polytropic laws, or at the bare least, choosing a value for γ that is consistent with the kinetic electron response.
As with any physical model, contrasting the predictions of the kinetic model against actual experimental measurements is a necessity to validate it and gain confidence on its accuracy. To this end, a set of dedicated experiments on the far plasma plume of a HT have been carried out at ESA-ESTEC, which were reported earlier [8,9]. The present paper compares the kinetic model with those experimental results, analyzing the differences found for six different thruster operating conditions. The validity of the model, and its different predictions, for these operating points is assessed by comparing the axial plasma profiles of the model and the experiments. While a kinetic solution is clearly superior to any polytropic approximation, due to their widespread use the paper concludes with an estimation of the polytropic cooling exponent γ in the different experimental cases, and discusses the limitations of such approximations to model the electron response in plasma thruster plumes.
The rest of the paper is structured as follows. For the sake of self-containment, section 2 summarizes the hypotheses and the main aspects of the electron kinetic model [4]. Likewise, section 3 reviews the experiments [9]. Section 4 explains the comparison procedure followed, and presents the comparison results to determine the error of the model. The discussion of the approximated polytropic exponent is then carried out in section 5. Finally, the main conclusions of this work are gathered in section 6.

Kinetic plasma plume model
The model of Merino et al [4] describes the steady state expansion of a plasma jet into vacuum, composed of singlycharged ions and electrons. The reader is directed to that reference for a full account of the model derivation; in this section, only the major aspects of the model, such as its assumptions and parameters, are summarized.
The plasma plume is considered to be collisionless, quasineutral, unmagnetized, and slowly-diverging (i.e. paraxial). The electric potential f is assumed to decrease monotonically downstream and radially, and is modeled as x are slowly-varying functions of x, which represent the potential along the plume axis (r=0) and the radial widening of the plume as it expands, to be determined as part of the solution. In essence, expression (1) prescribes the shape of the radial profile for f to be parabolic. The resulting electric field f - causes a small acceleration on the already-hypersonic ions that stream out of the thruster, but effectively confines most of the electrons, except for the highest-energy tail of their distribution. These free, escaping electrons are responsible for balancing out the ion current in the plume, making it globally currentfree. A typical confined electron trajectory involves many radial reflections between two axial reflections, or before returning to the plasma source.
Under these conditions, the mechanical energy E, the angular momentum about the plume axis p θ , and the adiabatic invariant J r given by the radial action integral [10]  . Ions are treated as cold species, and their bulk density and velocity are integrated directly from their paraxial continuity and momentum equations e e r , is averaged over the radial electron motion as explained in [4] e e r and solved for using the paraxial Vlasov equation This equation states thatf e is constant along x for each ( ) q E p J , , r , as long as , 0 x r defines the turning manifold at which electrons are axially reflected, and is given by where U eff is an effective potential for the axial motion of the electrons Observe that the turning manifold depends on J r and p θ only through , which allows to reduce the dimensionality of the model.
The turning manifold in equation (6) splits the electron phase space into four distinct regions, according to the connectivity of electron trajectories with the boundaries of the domain (x=0 and = ¥ x ): (i) reflected electrons, i.e. lowenergy electrons that are emitted from the plasma source and eventually reflected back to it by the potential fall along the plume; (ii) free electrons, i.e. high-energy electrons emitted from the source that can reach infinity downstream; (iii) doubly-trapped electrons, i.e. isolated electrons whose trajectories do not connect neither with x=0 nor with = ¥ x ; and (iv) empty regions, which connect only with = ¥ x and are therefore not occupied.
At the upstream boundary (x=0), conditions on the ion density n i0 and axial bulk velocity u i0 are provided. Likewise, the forward-going (v x >0) part of the radially-averaged electron distribution functionf e is required at that position. In the present work, a semi-Maxwellian electron population is prescribed at x=0. These boundary conditions suffice to obtain the solution for the ions and for the reflected and free electrons, but additional information is required to solve for the isolated doubly-trapped electron population. The model assumes that either the transient set-up process of the plume [11] and/or the occasional collisions help populate these regions with the same distribution as in neighbouring electron regions (fully-occupied doubly-trapped regions, i.e. α=100% in the nomenclature of [4]).
After assuming an initial guess for the functions f x (x) and h(x), the ion and electron models are used to obtain the ion and electron density n i and n e along the plume axis. The quasineutrality condition, -= n n 0 i e and the current-free condition are then used to compute the error of the solution and set up an iterative scheme that allows finding the self-consistent solution for f x (x) and h(x).
The resulting model has two remaining degrees of freedom, which can be parametrized in terms of the initial ion Mach number, M , i0 and the square root of the ion-to-electron mass ratio μ: In this work the value μ 2 =2.39×10 5 for xenon is used.
Observe that T e0 and n e0 depend on the backward-going part of the distribution functionf e at x=0, which is not known a priori and must be computed as part of the solution. Finally, note also that the initial values of f x (x) and h(x), i.e. f 0 and h 0 , as well as the constant C introduced in equation (1), can be fixed after solving the model, as long as the solution is rescaled appropriately.

Plasma plume measurements
This section overviews the main aspects of the experimental campaign carried out in April-May 2017 at the ESA/ESTEC Electric Propulsion Laboratory [9]. The measurements have been obtained for the plume of a 1.5 kW-class SPT-100-ML Hall thruster. This thruster has a channel outer diameter of 100 mm and an inner diameter of 69 mm [12]. The magnetic field at the exit is about 150G, and the thruster provides a thrust of approximately 90mN when used at its nominal operating point (300 V discharge voltage and 4.5 A discharge current). A schematic cutaway of the thruster can be found in Gascon et al [13]. The cylindrical Langmuir probes (0.2 mm diameter, 5 mm long tungsten wire, aligned parallel along the plume axis) were installed on two translation stages mounted inside the large CORONA vacuum chamber at ESA-ESTEC (2 m in diameter and 4 m in length). The first translation stage was mounted on the chamber rotating arm, allowing to scan angular position from −20°to +90°at distances between 550 and 750 mm from the thruster, while the second translation stage was placed further downstream, providing measurements along the plume axis 850 to 1550 mm from the thruster exit. A detailed account of the experiments and the setup used can be found in [9]. The current-voltage characteristic (i.e. I-V curve) of the Langmuir probe was measured using a Keithley 2440 sourcemeter, sweeping the probe voltage from −15 to +35V. A single method was used to obtained the plasma potential f, based on the maximum of the I−V curve first derivative. The plasma parameters, namely the electron density n e , and the electron temperature T e , were obtained from two different methods, (i) by fitting the appropriate expression of the current onto the I-V curve, and (ii) by using the Druyvesteyn method [14]. Three sources of uncertainty affect the derived plasma parameters: error on the measurement, i.e. repeatability and measurement uncertainty, error due to the direction of the voltage sweeps, i.e. hysteresis uncertainty, and error due to the method employed, i.e. method uncertainty. Measurement and hysteresis errors were quantified by taking multiple measurements at a fixed location and operating point (case A, ∼750 mm from the thruster along the plume axis), performing five measurements with both up and down sweeps, and obtaining the plasma parameters with each method. Independent measurements were consistent within 0.04 V for f,±1.5% for n e , and ±0.04 eV for T e ; differences between up and down sweeps were also consistent within 0.09 V, ±1.5%, and±0.08 eV, respectively. All errors were calculated as the±1σ deviation from the average value.
The method/model error arises from plasma phenomena not included in the standard probe model, such as the existence of non-zero ion flow and temperature, and partial electron magnetization. The two different methods were used to minimize this uncertainty. Both methods reproduce the same trend for n e and T e along the plume axis; however, a systematic offset in n e and T e was observed between the two, between 40%-50% for n e and 1-2 eV for T e . This offset indicates that the method error mostly consists of an absolute uncertainty on the true plasma parameters. The (relative) method uncertainty after removing the offset and taking the average of both methods was within 6%-9% for n e , and 0.10-0.16 eV for T e , depending on the operating case. Figure 1 shows an example of the measured electron density and temperature in the plume region obtained by fitting the I−V curve, for a discharge voltage of 300 V and a  The obtained values of electron temperature T e also present smooth variations, except for cases B, E and F. Noisier I−V curves were measured for cases B and E, most likely due to a less stable plasma for these discharge parameters. This significantly affected the plasma parameters obtained using the Druyvesteyn method, in particular the electron temperature, and explains the larger variations seen along the plume axis for these cases. For this reason, these two operating points also show the largest discrepancies between the two analysis methods. Faraday cup (FC) measurements were used in the test campaign to obtain the angular distribution of the ion current density in the plasma plume. Details about FC design and operation can be found in [15]. The graphite probe collimator aperture is 10 mm in diameter; it defines the collection area. The cup is 16 mm in diameter and 20 mm in length. Electrical insulation between the cup and the collimator is achieved with PEEK spacers. The FC is mounted onto a floating aluminum holder fixed to the chamber rotating arm. The FC faces the center of the thruster exit plane, i.e. it is parallel to the ion current streamlines. The HT to probe distance is x F =797 mm. The current was recorded with a Keithley 2400 calibrated sourcemeter. The latter is also used to polarize the cup at −100 V for all measurements. The FC was positioned in the horizontal plane that contains the SPT-100-ML Hall thruster centerline. Ion current density measurements were performed over a hemisphere. In this work the Hall thruster was not aligned with the rotating arm pivot. As a consequence, the FC angle with respect to the thruster exit plane does not correspond to the arm angle, but at 0 • when the arm is aligned with the thruster axis. Correction has been applied to all measured current angular distributions to retrieve the effective distributions. Figure 2 shows raw and corrected ion current density angular profiles for case F. Figure 3 shows the current-voltage curve at various angles. The results show the great advantage of using Faraday Cups to measure the ion current density, since the saturation current does not depend on the voltage. This means that the plasma sheath expansion effect is cancelled; thus current density measurements are more accurate with FCs compared to other probe types. Although the error in the ion current value is hard to estimate, it is known that several effects must be taken into account: secondary electron emission in the cup, plasma perturbations due to probe-plume interaction, point source assumption and alignment. More information is available in [15] and references therein.

Comparison results
The procedure used to compare the experimental results with the model can be summarized as follows. Firstly, the available experimental far plume data for each test case is prepared for the comparison. Since the spatial locations where the Faraday cup data and the Langmuir probe data were obtained are not identical, interpolation is carried out when needed to evaluate a plasma property. While the absolute values of the available plasma properties, f, n e and T e , are not essential for the comparison and validation process, the kinetic model must reproduce correctly the relative values of each of these variables. In order to limit the method error and offset between the I-V curve method and the Druyvesteyn method, the average of n e and T e obtained from Langmuir probe measurements processed with each of them is computed and used as the reference data.
Secondly, we determine the input parameters of the model for each case. The value of M i0 is obtained by using the ion current density from the Faraday probe at point F, x F =797 mm, i.e., j i (x F ). The value of ( ) n x e F at that point is then interpolated from the Langmuir probe data, and ion velocity is determined as   energy equation As it can be observed, the reference electron density n e , ion velocity u i , electric potential f, and electron temperature T e increase with the discharge voltage V d , with the exception of case B, which has the worst quality experimental results and is slightly off this trend. The increase of n e with V d can be attributed to a less divergent plume as V d increases. The estimated ion velocity at x F following the procedure above is seen to be consistent in all cases with the estimation based on . Lastly, it is seen that electron density n e increases and u i decreases with increasing  m, while there is no clear trend with  m for f and T e . The behavior of the ion Mach number M i0 , which results from the combination of the behavior of u i and T e , is also non-trivial; however, overall it can be stated that it increases with the discharge voltage, V d .
Thirdly, the beam-width function, ( ) ( ) h x h x P , is computed from the experimental data using the paraxial continuity equation: The beam-width function was inspected for all test cases and seen that it could be accurately approximated as a straight line in the measurement region with a relative RMS error of less than 0.1%. This is consistent with a nearly-conical expansion in the far plume region. Hence, a linear h(x) fit with a constant slope h x d d (dependent on the operating case) is used in the comparison in the range where experimental data are available. The value of h x d d , which represents the plume divergence rate, is seen to decrease with increasing discharge voltage V d . This is consistent with a more collimated plasma beam as the discharge voltage increases.
Fourthly, the remaining parameters of the model, i.e. ( ) f h x n , , P e 0 0 , and C, are determined for each test case by minimizing the mean square error in potential, density and temperature with the experimental data at the plume axis, for all data points between x=550 and x=1550 mm where Langmuir probe data are available. These parameters are shown in rows 10-14 of table 1. With the exception of the previously mentioned noisy case B, the general trend of the upstream potential, f 0 , and electron density, n e0 , is to increase with the discharge voltage as could be expected.
It should be noted that the kinetic model returns the axial, radial, and azimuthal electron temperature Once all the model parameters are fixed, simulations have been carried out for each experimental case. Below, the comparison and discussion of the results is presented.

Radial electric potential profile
One of the hypotheses of the kinetic model is that the electric potential profile f has a parabolic shape in the radial direction. The radial electric potential profile at x 1 =550 mm is shown in figure 4 for test case A. This figure is representative of the rest of the test cases, which are not shown in the sake of conciseness. As it can be observed, the experimental data are well approximated by the parabolic electric potential profile up to about r=150 mm; from there outward, the measured potential decreases at a lower rate than the one considered in the model.
As the majority of the electron population has low energy E, most of the electrons cannot travel to large values of r, and hence they do not perceive the difference between the radially-parabolic and the measured electric potential. The deviation from a radially-parabolic electric potential profile is expected to introduce an error in the model due to the change in the mathematical expression of the radial action integral, J r , of the higher-energy electrons that can travel far in the radial direction before being radially reflected. In essence, the higher value of the potential at large r means that those electrons with high J r will display a larger radial oscillation period than in the parabolic case. Future work must assess the importance of this assumption further by considering other radial profiles for the electric potential model.

Plasma properties along the plume axis
The central part of the plots in figures 5-7 compares the electric potential f, the electron density n e , and the electron temperature T e between the fitted model and the experimental results in the region of measurement for test cases A through F. We can conclude that the potential, electron density, and electron temperature on the plume axis are well described by the model in all cases. The comparison error in electric potential is below 0.6 V everywhere. This error is smaller than 10%-20% of the local value of T e in all cases, except in the noisy case B where it reaches 30%. The relative error of electron density n e is less than 5%-10% for all the experimental cases.
Regarding the electron temperature measurements, they are generally more noisy than the other variables. This is evident from the erratic experimental curves in these figures, especially in cases B, D, E. Notwithstanding this, the model clearly reproduces the average trend of the experimental data.
The data show that electron cooling is taking place in the near-collisionless expansion, which is one of the major predictions of the kinetic model. It is worth noticing that the experimental error in the electron temperature T e affects the determination of the constant C in the model, a fundamental sizing parameter for the radial electric potential profile, and thus an error in its determination affects the electric potential and electron density as well. Since the parameter fitting procedure followed here ensures that the error is as small as possible in f, the electron density n e , and the electron temperature T e simultaneously, these three fits are not independent from each other. The differences between model and experiment are in all cases comparable to the expected experimental error of a Langmuir probe in a plasma thruster plume, and indeed, most of the differences have a noisy oscillatory behavior, which can be attributed to the measurements themselves.  floated with respect to the vacuum chamber walls, to which the potential is referenced (f=0). Two major differences between the free plasma expansion in space and in a vacuum chamber are the higher background pressure and the limited plume length in the latter case. These effects make it more difficult to validate the kinetic model with laboratory experiments. Again, it is noted that the validity of the model relies on the expansion being non-magnetized and collisionless. Thus, as the model is extrapolated toward the exit plane of the thruster (x=0), the error committed is expected to increase, as these assumptions are not longer satisfied. In fact, h(x) is not expected to have the same slope h x d d in this near region as in the far region. Notwithstanding this, the model can be used to estimate f 0 , n e0 , and T e0 , the upstream values of the electric potential, electron density, and the electron temperature. Indeed, the former two were obtained as part of the parameter fitting process in section 4 and the three of them are shown in table 1. These values should roughly correspond to the potential, density and temperature close to the thruster, but they cannot be directly related to any specific location for the reasons above. As it can be observed, with the exception of test case B, f n , e 0 0 and T e0 increase with increasing discharge voltage V d . Increasing the xenon mass flow rate increases the electron density n e0 . However, it also increases the background pressure in the vacuum chamber, which translates into a larger plasma collisionality. This goes against the collisionless assumption in the kinetic model, so it can be expected that cases E and F will be worse represented by the model.
As in the laboratory the presence of the chamber wall interrupts the plasma expansion to infinity, the plasma forms a sheath at the chamber wall downstream to ensure a net current-free plume, such that the total potential fall from the thruster to the chamber wall, i.e. including the sheath, is essentially the same as the total potential fall in space from the thruster to infinity. In other words, f ¥ should roughly correspond to the potential of the vacuum chamber, f=0. As it can be observed in table 1, f ¥ is greater than 0 in all test cases, but within one electron temperature T e0 from it, i.e. f < . This weak consistency further supports the validity of the model, but also provides evidence that it cannot be extrapolated reliably to the sheath and wall region, as expected.

Approximate fluid closure law
As previously mentioned, most of the existing plasma plume simulation codes follow a fluid description for the electrons and many of them use as closure relation a simple law relating the electron density and temperature by means of a simple polytropic exponent, γ, i.e. 1 . The kinetic results of the model [4] can be used to compute an approximated, effective electron polytropic cooling exponentḡ for Figure 6. Experimental (red segment lines) and fitted model (blue lines) curves for the plasma potential, f, electron density n e , and electron temperature T e , along the plume axis (r=0), for test cases C and D. The difference between the two lines in the region of measurement is displayed in the inset plots.
increasing differences in the plume radial periphery. The deviation from the assumed profile introduces a computation error in the radial action integral of electrons, J r , which is expected to affect only high energy electrons, i.e. only a minority of the population. Kinetic models with other radial potential profiles, closer to the measured ones, can be developed to better approximate this aspect of the problem.
The agreement on the fundamental plasma parameters along the plume axis (potential f, electron density n e , and electron temperature T e ) has been evaluated and discussed using the available measurements, in the range x=550-1550 mm. The error in potential and electron density is small in all six test cases and within the experimental uncertainty. The evolution of T e clearly follows the trend recovered experimentally too, but the measurements are noisier.
The departure of the laboratory test conditions from inspace conditions presents a hurdle toward the validation of the model. In particular, the higher background pressure and presence of the vacuum chamber wall that interrupts the plume expansion to infinity can affect the results of the comparison. It is noted that the electron response in the collisionless plume is global, meaning that an effect downstream can alter the expansion upstream, and vice-versa. However the model shows a fitted f ¥ value that is only slightly higher than 0, consistent with the potential of the vacuum chamber wall.
Finally, an approximated, effective electron polytropic coefficient was computed that gives the same total potential fall than the model. The values computed for the different operating conditions tested are in the rangeḡ = 1.26 e -1.31, in agreement with other values reported in the literature.