Electrodeless plasma thrusters, such as Helicon (HPT) [1–3] or Electron Cyclotron Resonance (ECR) [4, 5] thrusters, represent very promising solutions in the plasma propulsion community to overcome the inherent lifetime limitations of more mature technologies, such as the electrostatic gridded ion thrusters (GITs) or the Hall thrusters (HTs). Nevertheless, these relatively new engines currently feature very limited thrust efficiencies, which rarely surpass 10–15% [6] and need to be increased significantly before they can be considered a serious competitor of GITs or HTs. In this maturity raising process, the correct experimental characterization of their emitted current-free, quasineutral and magnetized plasma plumes represents a fundamental step. In order to better interpret the observed experimental curves, numerical simulations represent a relatively cheap but powerful comparison tool that can help explore quickly the effects of different physical processes occurring either inside the thruster or in the near plume.

Numerical models can also be employed to simulate the measurements of a real diagnostics tool used in the experimental campaigns, through the so-called synthetic diagnostics simulations. This generally requires the modeling of the plasma - probe pair together, with their mutual influences. As a result, these synthetic simulations provide an even better understanding of the involved physics, as well as the characterization of some measurement errors induced inherently by the probe definition (geometry, bias voltages, etc … ). The usual diagnostic tools employed in the experimental characterization of plasma thrusters plumes include, among others, retarding potential analyzers (RPA) [7], Faraday probes (FP) [8], Langmuir probes (LP) [9, 10], emissive plasma probes [11], laser induced fluorescence techniques (LIF) [12]. Of particular relevance is the RPA, which has been used for decades to measure relevant plasmas properties in both space plasmas [13] and plasma thruster plumes [14, 15], such as the ion velocity distribution function (IVDF). An RPA measures the electric current, mainly carried by positively charged ions, which traverse a set of grids with holes and reach a collector plate. This collected current is measured as a function of the voltage applied to one of the grids, hereinafter referred to as retarding grid. This bias imposes an adverse potential field for the probed particles, acting as an energy filter. Consequently, by varying the grid voltage in a given sweep range (depending on the expected ion energy distribution), a characteristic I-V curve can be measured, and through a simple differentiation technique, the IVDF is finally obtained. Given its simplicity and usefulness, improvements of the RPA design that reduce measurement errors have been constantly proposed for the last decades [16, 17]. In this context, the RPA design optimization can greatly benefit from numerical synthetic simulations: in fact, the effects of varying the geometry of the grids and even advanced RPA concepts [7] can be quickly evaluated without having to build a costly experimental prototype.

Numerical models for the magnetized plasma plume expansion through a magnetic nozzle are of different types. A first approach is that of the kinetic approaches, which can be split into methods solving a simplified Boltzmann’s equation for the ion and electron distribution functions in a low dimensional space (typically 1D, [18]), and full particle-in-cell (PIC) models [19, 20], featuring particle ions/electrons, collisions through Montecarlo approaches [21], and generally limited to 2D to reduce the computational cost. An alternative approach, which is computationally cheaper, is that of hybrid codes, in which electrons are modeled as a fluid, while neutrals and ions are followed as macro-particles of a PIC sub-model featuring Montecarlo collisions. This enables the use of a much larger time step compared to full-PIC codes, since the Courant-Friedrichs-Lewy (CFL) [22] condition now applies to ions, which are much heavier and hence slower than electrons.

The 3D code EP2PLUS [23] makes use of this hybrid approach, and has already been used in different plasma thruster scenarios, from S/C-debris-plasma interaction in an ion beam shepherd scenario [24] to a plasma plume expansion under the geomagnetic field [25], including also simulations of 3D Hall thruster near plumes [26] and ion grid optics [27]. The magnetized electron model of Refs. [25, 26] is particularly appropriate to simulate the plasma expansion across the divergent part of the magnetic nozzle of an electrodeless plasma thruster. Such a model features a quasineutrality assumption and solves the coupled electric current continuity and electron momentum balance equations, retaining both the electron collisions term and the very relevant magnetic force term, j _e × B , where j _e is the electron current density and B is the magnetic induction field. The unmagnetized electron model of Ref. [27], on the other hand, is particularly suited for the synthetic simulation of an RPA: Boltzmann electrons are assumed and a non-linear Poisson’s equation is solved to properly assess the effects of space charge inside the RPA and of the plasma sheath forming in front of the electron repelling grid.

In this paper, we aim to better understand the physics of both a magnetic nozzle expansion and a retarding potential analyzer, and, at the same time, to partially validate their numerical models by direct comparison with experiments. In particular, this work attempts to reproduce numerically a plasma plume emitted by a Helicon thruster prototype [2], which has been characterized experimentally through the use of an RPA, a Faraday probe and a radiofrequency compensated Langmuir probe, although the main focus is on the RPA results. With the first magnetized plume model and some experimental inputs, a large-scale simulation is first run to estimate the plume properties downstream of the magnetic nozzle throat and, in particular, at the RPA location. Then, the unmagnetized model is applied to a smaller spatial scale simulation (corresponding to a single RPA orifice) to obtain a synthetic I-V curve and hence the IVDF as a function of the axial ion energy. Numerical results are then compared to experiments and the differences are discussed.

The paper is structured as follows: Section 2 introduces the experimental setup used to characterize the HPT plume, while Section 3 introduces the numerical models for both the magnetized and unmagnetized simulations. Then Sections 4, 5 report the numerical simulation results obtained with the two models, and Section 6 compares the obtained results with the experiments and discusses the observed differences. Finally, Section 7 outlines the main conclusions of the study.

The HPT prototype used in the frame of this work is the so-called HPT05M prototype, developed together by UC3M and SENER Aeroespacial [2, 28, 29]. This was built according to the classical architecture and is composed of the following elements: a single solenoid (S1 in Figure 1A), which generates an axisymmetric convergent-divergent magnetic field with its peak, B _max, at the geometric center of the solenoid; and a dielectric tube made of quartz, 30 mm inner diameter, 1.5 mm thick, and 150 mm long. The tube outlet coincides with B _max and the tube is axially aligned with the magnetic field. At the other end of the tube, there is a ceramic injector, multi-hole “shower-type”, made of Macor (machinable glass ceramic), which allows to inject the working noble gas (Xenon, Krypton or Argon).

Around the tube, there is a half-turn right-hand polarized helical antenna, which is centered between the magnetic throat and the mentioned injector plane. Solenoid S1 is fed with DC current, generating magnetic fields up to B _max = 1500 G, which is the nominal field considered here, while the antenna is fed with radio-frequency (RF) power at 450 W and 13.56 MHz. An L-type matching network, customized by SENER, is placed in between the RF power amplifier (SEREN IPS HR2000 series) and the transmission line (50 Ohm - coaxial cable) that connects to the antenna. This setup allows to work always at optimum conditions, or with null reflected power, so that 450 W can be considered as the transmitted power to the antenna-plasma system.

The thruster breadboard has been tested within a vacuum chamber designed specifically to characterize Electric Propulsion Plasma Thrusters up to 1.5 kW, and located at the UC3M facilities. The vacuum chamber consists of a stainless-steel 304 (non magnetized) vessel of 1.5 m inner diameter and 3.5 m long. This is equipped with three different vacuum technologies: a dry mechanical pump Leyvac LV80 with pumping speed of about 80 m³/h, a pair of turbo-molecular pumps, Leybold MAGW2.200iP with 2,000 L/s of pumping speed each, and three cryo-panels, Leyvac 140 T-V, all systems from Leybold GmbH. The total pumping speed is about 37,000 L/s Xe, reaching an ultimate pressure of 10^–7 mbar in dry conditions. The operational pressure is roughly 2⋅10^–5 mbar, at 20 sccm of Xe. Xenon gas flow rate is set/measured by a Bronkhorst EL-FLOW Select Xenon calibrated mass flow controller, with a resolution of 0.1 sccm and 100 sccm full range.

The RPA, a Semion unit by Impedans, is mounted on a multi-probe stand, which is held on a two-axes polar translation stage, designed by UC3M. Figure 1B shows a picture of the full setup under test, the HPT05M firing with Xenon and the arm system scanning the expelled plasma plume. Probes can be moved along the radial distance R and angle θ, on the horizontal plane, coplanar to x − z (see Figure 1A for further reference). The probes vertical position is adjusted before the tests in such a way that the horizonal plane x − z contains the thruster setup axis line. The radial and angular resolutions of the arm system are 1 mm and 1° respectively, with uncertainties ± 0.3 mm and ± 0.3°. The stages range is 0–400 mm for the radial direction R and ± 90° for the angular, sweeping fully a semicircle in front of the tested thruster unit. The RPA from Impedans [30] then features a collector plate and 4 grids: an orifice grid; an electron repelling grid; a retarding potential grid to filter the different ions, according to their axial energy; and a secondary electron suppression grid, which is used to reduce the spurious signal due to secondary electrons emitted from the collector plate due to ion impacts. Additional details on the RPA geometry are given in Section 3.2 where the RPA simulation setup is defined.

Apart from the RPA, the probes setup also includes a Faraday probe and a radio-frequency compensated Langmuir probe. The first one measures the ion current density profile j _i(θ, R = const) through the plasma plume. This Faraday probe (FP) is made of stainless steel with an electrode collector of 10 mm diameter surrounded by a guard ring, 20 mm outer diameter, and 10.5 mm inner diameter. Both electrodes are biased to −50 V with respect to ground. Concerning the RF compensated Langmuir probe [10] (LP), this is used to measure several plasma properties, including the ion density, the electron temperature and the plasma potential. Its main electrode is a tungsten rod, with an exposed length of 2 and 0.27 mm in diameter. The experimental results, in terms of the IVDF (RPA), the ion current density (FP) and the plasma density (LP), are plotted together with the numerical results for an easier comparison in Figures 11 and 12.

The operating conditions of the HPT05M have been maintained constant for all the experimental results presented next, if not mentioned otherwise. The Xenon mass flow rate has been set to 10 sccm and the RF power to 450 W (with null reflected power). The magnetic topology consists of a convergent-divergent magnetic field with a peak strength of about 1500 G at the throat. The field is measured within the plume region using a 3-axis Gaussmeter and the same arm system holding the diagnostics described above. The measured field intensity and direction are shown in Figure 6A for the considered simulation domain of the magnetized expansion.

The models for both the magnetic nozzle expansion and the synthetic RPA simulations follow a hybrid approach, so that ions and neutrals are considered as macro-particles of a PIC sub-model, while electrons are treated as a fluid, subject to conservation equations. Therefore, the simulation time step is given by the Courant-Friedrichs-Lewy condition [22] applied to ions, meaning that the fastest ion in the simulation must cross less than one cell per time step. A full description of the PIC algorithms can be found in Ref. [23] and is not repeated here. In the following, we focus primarily on the electron models considered in both scenarios and we only summarize the main aspects of the PIC sub-models, such as the relevant boundary conditions for ion/neutral macro-particles.

In this section, the reported results are averaged over 500 PIC time steps, in order to reduce the noise. The magnetic field intensity and direction and the corresponding effective Hall parameter are reported, at the meridional plane y = 0, in Figures 6A,B.

The other most relevant simulation results at y = 0 and at the throat cross section (z = 0) are then shown respectively in Figures 7, 8.

The electric potential and the corresponding electric field are shown in subplots (a) and (b) of both figures. Most of the electric potential decay occurs in the first centimeters of the expansion, downstream of the nozzle throat, while the radial profile is quite flatter than the one found for a simple ambipolar and unmagnetized expansion. The magnetic field has the effect of reducing the radial decay of the electric potential thus focusing better the plume expansion, as clearly shown by the electric field plots in both the meridional (y = 0) and nozzle throat (z = 0) sections.

The electron density (which is equal to the ion density in this quasineutral scenario with only singly-charged ions) is shown in Figure 7C and decays from an initial value close to 10¹⁹ m⁻³ at the throat to around 7⋅10¹⁵ m⁻³ downstream at the centerline.

In order to evaluate the effect of collisions in the expansion, the ratio between the number density of source ions n _{i, src} (i.e., ions injected from z = 0) and the total plasma density n _i = n _e is shown in subplot Figure 7D. While at the throat the near totality of ions come from the helicon thruster source, this percentage reduces progressively to around 60% downstream at the centerline, and much more quickly radially outward at the throat section. The responsible for this deviation are the ions produced by both ionization of neutrals and charge-exchange between source ions and neutrals.

Finally, Figures 7E,F, 8C,D show respectively the electron current density (with the direction of the electron flow) and the ion current density at y = 0 and z = 0. While ions are only dimly magnetized, electrons feature a dominant azimuthal current density (in the x-y plane), which is what ultimately produces the magnetic thrust ( j × B ≃ B _r j _θ 1 _z, with polar components referring to the x-y plane). This is evident at the throat cross section, but also in Figure 7E, where the size of the arrows becomes negligible far from the centerline because of the dominant y (azimuthal) component. It is underlined that the magnitude of the total electron and ion current densities is not exactly equal at the lateral domain boundaries (see Figures 7E,F). This is not in contrast with the imposed boundary condition, which consists in a zero electric current density in the direction normal to the boundaries, i.e. j e + j i ⋅ 1 n = 0 , thus allowing a circulation of electric currents parallel to them.

The electron number density along the centerline is reported by the solid line of Figure 11A, while the ion current density versus the angle θ is shown by the solid line of Figure 11B at a distance R = 40 cm from the nozzle throat. Finally, the simulated total thrust force is 5.2 mN, which has the following contributions: 3.9 mN of injected momentum flow from the upstream boundary (of all species), 1.2 mN of magnetic force on electrons (∫ j _e × B dV, with dV the differential volume), and 0.1 mN of magnetic force on ions (∫ j _i × B dV). As expected, the ion magnetic force contribution is much smaller than that of the electrons.

Figure 9 shows the RPA simulation results in terms of electric potential, ion density and ion trajectories, at the y′ = 0 plane and for three different retarding grid potentials ϕ _r relative to the local plasma plume: 0 (left column), 25.5 (middle column) and 49.5 V (right column).

At the densities and electron temperature of the local plasma plume (respectively less than 10¹⁶ m⁻³ and 1.5 eV), the Debye length is in the order of 0.1 mm, which is much larger than the side of the hole. This yields a nearly flat potential potential profile along x′, y′ in all cases as observed in Figures 9A–C. As the electric potential is gradually increased, the ion density map clearly changes (Figures 9D–F). An increasing peak of ion density forms at the retarding grid hole center, due to the reflection of an increasing fraction of ions. At ϕ _r = 49.5 V, the near totality of ions is reflected backwards, and the highest peak can be observed. This can be appreciated also in Figures 9G–I, showing the individual ion trajectories. As the potential is increased, ions start to be reflected back towards the plasma plume and do not reach the collector.

By computing the ion current reaching the collector as a function of ϕ _r, it is straightforward to obtain the I-V characteristic (i.e. the collector current evolution versus the retarding grid potential) of a single RPA hole, which is displayed in Figure 10A. This characteristic gradually decreases and shows the maximum negative slope at the location of the measured distribution peak. By simple differentiation, and applying Eq. 9, it is finally possible to obtain the simulated IVDF.

Some relevant information can also be found in Figure 10B, which shows the corresponding evolution of the currents collected by the 3 RPA grids. The current collected by the electron repelling grid (green line) increases monotonically, as an increasing fraction of ions gets reflected backwards by the retarding grid and ends up hitting this grid from the inside. Therefore, the current collected by the retarding grid (red line) decreases rapidly close the IVDF peak, as expected. Finally, the current collected by the SEE suppression grid has a non-monotonic behavior reaching a maximum value at around 40 eV and then decreasing. This anticipates the collector current decay for ϕ < 40 V, while delaying it for ϕ > 40 V. The consequence is that the simulated main IVDF peak appears shifted towards lower energies and with an overestimated width as shown in Figure 13, as further discussed in the next section.

The total thrust from the magnetic nozzle expansion simulation (5.2 mN) compares very well with the measured thrust value (5.1 mN), which was obtained directly through a specific thrust balance for HPTs, described in Ref. [35].

Figure 11 then shows a comparison between the experimental results and the magnetic nozzle simulation. Subplot A is for the plasma density along the plume centerline, while B is for the normalized ion current density along the azimuth angle θ, at a fixed radial position, R = 40 cm from the magnetic nozzle throat (or thruster outlet section).

Regarding the number density, the experimental I-V characteristic curve of the radio-frequency compensated Langmuir probe is post-processed according to the most convenient model of the Langmuir probe theory [9], taking into account the different plasma-probe regimes depending on the ratio between the local Debye length and the probe radius. Simulation results appear to underestimate the plasma density by a factor of 2, although the relative error seems to be reduced downstream. Such deviations, however, are justifiable by the expected Langmuir probe measurement errors (both bias and stochastic errors), which can be even larger than the observed differences. For what concerns the ion current density profiles (Figure 11B), both experimental and simulated profiles are normalized with respect to their maximum value occurring at θ = 0°. Since the Faraday probe points to the thruster outlet section for all angular positions, the measured ion current density is assumed to be always aligned with the radial direction R of the polar reference frame. Simulations agree very well with the measurements especially for θ < 10°, and tend to slightly overestimate the distribution at larger angles.

Figure 12 then shows a comparison between the experiments and the magnetic nozzle simulation of the IVDF expressed as a function of the ion kinetic energy in the radial direction R, at R = 40 cm from the thruster outlet and at 3 different θ angles (refer to Figures 1, 2).

The IVDFs are normalized in such a way that the results from experiments and magnetic nozzle simulations present the same peak value (1.0) at θ = 0. Very interestingly, both experiments and magnetic nozzle simulations agree very well in the energy of the observed IVDF peaks. The first one, at ∼ 45 eV, is due to source ions accelerated through the full electrostatic potential drop of the divergent magnetic nozzle of Figure 7A, while the second one, at ∼ 15 − 20 eV, is due to slow ions produced in the plume, due to either CEX or ionization, with simulations underestimating its importance compared to experiments. As we increase the θ angle, the relative weight of the low energy population with respect to the high energy one increases slightly, and this is confirmed by both experiments and simulations. A significant discrepancy can be observed at intermediate energies between the two mentioned peaks, with experiments giving a larger IVDF there. While this also might be caused by the assumption of Maxwellian and polytropic electrons, a very likely reason is the assumed Maxwellian distribution of injected ions at the upstream magnetic throat section of the simulation. This might indeed differ from the ‘unaccessible’ real ion distribution, which could present a more populated low energy tail. Another remarkable difference is that the magnetic nozzle simulations also predict a non-negligible IVDF at very low thermal energies, which is due to the existence of a background plasma created by CEX collisions with local neutrals. This is totally missed by the measurements, so that RPA experimental curves may not be fully reliable at very low energies, as also explained below, by comparison with the RPA simulation results. For the θ = 0 case, Figure 13 shows a comparison of the IVDF predicted by the RPA simulations with that of both the experiments and the magnetic nozzle simulations.

The RPA simulations take, as upstream boundary conditions for particles at z′ = 0, the velocity distribution function extracted from the nozzle simulation results, so that the targeted IVDF curve is the black solid line in both Figures 12A, 13. Two simulated IVDF curves are shown: the first one (grey dashed line with circle markers) is obtained directly from the I-V characteristic of Figure 10A, with f ∝ − dI _c /dϕ _r , while the second one (blue dotted line with triangle markers) is a corrected IVDF computed by applying Eq. 9 not to I _c alone but to the sum of the collector, SEE suppression and retarding grids currents: I _c′ = I _c + I _r + I _see, so that f ∝ − dI _c′/dϕ _r . Clearly, this correction permits to filter out most of the parasitic effects due to the evolution of the grids currents with ϕ _r as depicted in Figure 10B, so that the corrected IVDF reproduces almost exactly the input IVDF from magnetic nozzle simulations, thus verifying the correct implementation of the RPA model. The knowledge of the grids currents (if available) would thus permit to correct the RPA measured curve. On the other hand, as already anticipated in Section 5, the uncorrected IVDF (grey dashed line with circle markers) captures correctly both the intermediate ion energy peak and the very low energy tail (unlike the experimental curve). Nevertheless, the simulated curve features a slightly lower energy and an overestimated width of the main IVDF peak, as previously commented. This has the effect of artifically increasing the IVDF values at intermediate energy levels, which might partially explain the discrepancies between simulations and experiments. This observed parasitic effect, however, may not be the same in reality, since it is strongly dependent on whether or not the holes are perfectly aligned, which is not necessarily the case of the real RPA. Non-aligned holes would feature a lower global transmission, but would also suffer from lower parasitic effects caused by the ion current collected by the SEE suppression grid.

This paper has presented two different numerical models for (i) a plasma plume expansion through a divergent magnetic nozzle and (ii) the simulation of a retarding potential analyzer with an arbitrary ion velocity profile. Such models have been applied to a plasma plume emitted by a Helicon plasma thruster prototype, tested at UC3M facilities, and their results have been compared with the corresponding experimental measurements.

Both models are hybrid models, meaning that the heavy particles are followed with a particle-in-cell method, while electrons are treated as a fluid. In the model for the magnetic nozzle expansion, magnetization and collisional effects are retained for electrons, while in the RPA model, a simple electron Boltzmann distribution is assumed.

The numerical simulations of the magnetic nozzle expansion appear to be good at predicting the measured angular ion current density profile at 40 cm distance from the source (and hence the divergence angle) and also the energies of the two main peaks of the IVDF, observed at the RPA location (15–20 and 45 eV). This permits to conclude that the lower energy peak is essentially due to CEX ions or slow ions from neutrals ionization, both produced in the plasma plume. Nevertheless, the IVDF at intermediate energies between these two peaks is underestimated with respect to experiments, a fact that might be caused by either the model approximations (e.g. a polytropic electrons thermodynamics) or an ion distribution at the nozzle throat differing from a Maxwellian distribution, as assumed in the simulation. Therefore, the real distribution of ions at the throat probably includes an important low energy population tail.

The RPA synthetic simulation, on the other hand, has permitted to evaluate precisely the parasitic effects on the measured IVDF of the ion currents reaching the RPA grids. Since such currents vary with the applied retarding grid potential, the resulting simulated IVDF curve shows a main energy peak at a lower energy and with a larger width, compared to the expected IVDF.

Regarding future work, the magnetized electron fluid model for the nozzle expansion needs to be improved by properly including an electron energy conservation equation, thus overcoming the current limitations due to the polytropic assumption. Also, future simulations should cover a larger angle extension downstream, which would require, however, a much larger computational cost. Finally, regarding the RPA studies, it would be interesting to apply the model to an impinging plume with a relatively large angle relative to the RPA axis to evaluate any induced parasitic effect in that scenario.