Skip to main content


Front. Phys., 14 April 2022
Sec. Low-Temperature Plasma Physics
Volume 10 - 2022 |

Magnetic Nozzle and RPA Simulations vs. Experiments for a Helicon Plasma Thruster Plume

  • 1Institute for Plasma Science and Technology (ISTP, CNR), Bari, Italy
  • 2Equipo de Propulsión Espacial y Plasmas (EP2), Universidad Carlos III de Madrid, Leganés, Spain
  • 3Departamento de Bioingeniería e Ingeniería Aeroespacial, Universidad Carlos III de Madrid, Leganés, Spain

The experimental characterization of electrodeless plasma thrusters with a magnetic nozzle is fundamental in the process of increasing their maturity to reach the industrialization level. Moreover, it offers the unique opportunity of validating existing numerical models for the expansion of a magnetized plasma plume, and for the synthetic simulation of diagnostics measurements, like those of a retarding potential analyzer, which provides essential information regarding the ion beam energy distribution function. Simulations to experiments comparison ultimately enables a better understanding of the physical processes behind the observed experimental curves. In this work, input experimental data of a Helicon plasma plume is used to simulate both a magnetic nozzle expansion in the divergent field region, and the corresponding measurements of a retarding potential analyzer, through dedicated small-scale simulations of this diagnostics tool. Magnetic nozzle simulation and experimental results agree well in terms of the angular distribution of the ion current at 40 cm distance from the source, and also in the prediction of the energies of the two main peaks of the ion energy distribution function: a first one at 45 eV due to source ions, and a second one, at 15–20 eV, due to ions from charge-exchange and ionization collisions in the plume. Finally, the small-scale simulation of the retarding potential analyzer permits to assess the parasitic effects caused by the ion current collected by the different analyzer grids. The inclusion of the retarding and electron suppression grids currents in the overall I-V characteristic is shown to correct almost entirely these effects on the obtained ion velocity distribution.

1 Introduction

Electrodeless plasma thrusters, such as Helicon (HPT) [13] 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, je × B, where je 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.

2 Experimental Setup

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, Bmax, 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 Bmax 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).


FIGURE 1. (A) Experimental setup sketch (not in scale), and (B) HPT05M under test, firing with Xenon. The RPA is mounted on the multi-probe movable arm, which can be moved according to the polar frame Rθ, centered at {z, x, y} = 0, and coplanar with xz. The azimuth angle θ is measured with respect to z axis. Both reference frames are centered at the plasma discharge tube outlet, which coincides with the center of the solenoid S1. The magnetic nozzle throat where B = Bmax is on the z = 0 plane, and the magnetic field is convergent for z < 0 and divergent for z > 0. In subplot (B), the probes arm holder includes the RPA at the center and two different Faraday probes at its sides.

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 Bmax = 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 m3/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 xz (see Figure 1A for further reference). The probes vertical position is adjusted before the tests in such a way that the horizonal plane xz 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 ji(θ, 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.

3 Numerical Models and Setups

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.

3.1 Magnetic Nozzle Expansion Model and Simulation Setup

The domain considered for the simulation of the magnetic nozzle expansion is schematically shown in Figure 2 at the y = 0 plane, while the corresponding numerical mesh is given in Figure 3 at two different cross sections (y = 0 and z = 0). The mesh consists of 71 × 71 × 101 nodes and is a deformed structured mesh with a circular cross section and a growing spacing along z in order to reduce respectively numerical boundary effects affecting fluid electrons, and PIC noise downstream (reduced by the increasing cell size), for a fixed computational cost. The simulation parameters for both the PIC and electron sub-models are then reported in Table 1.


FIGURE 2. Simulation domain for the magnetic nozzle expansion simulation, shown at the y = 0 plane. The expansion is along the z axis, while the magnetic nozzle throat is at the left upstream boundary.


FIGURE 3. Structured mesh used by both PIC and electron fluid in the magnetic nozzle simulation, shown at (A) y = 0 and (B) z = 0. For the sake of clarity, fixed computational coordinates lines are shown every 5 cells along all directions. The expansion is along the horizontal z axis in subplot (A), while subplot (B) represents a mesh cross-section at the injection boundary.


TABLE 1. Simulation parameters for the magnetic nozzle simulation.

Regarding the particle-in-cell model, and referring to Figure 2 and Table 1, both Xenon singly charged ions and neutrals are simulated (no doubly charged ions are considered) and injected from the magnetic nozzle throat located at z = 0, with sonic conditions. So, ions are injected with an average axial velocity equal to Bohm’s velocity uB=Te/mi = 3.0 km/s and an injection temperature Ti = 0.5 eV, while neutrals with an average sonic velocity equal to 0.25 km/s and an injection temperature of 580 K. Ions and neutrals further collide between themselves through Montecarlo Collision algorithms (MCC), which include the effect of charge-exchange collisions [23]. These transform fast ions and slow neutrals into respectively slow ions and fast neutrals. Additionally, neutrals are ionized by the electron fluid, using the algorithms presented in Ref. [23] and based on ionization rates that are a function of the local electron temperature, thus representing an additional volumetric source of slow ions, since these are created with the velocity distribution of their predecessor neutrals. A background neutral density of 2.5⋅1017 m−3 is also considered for collisions, in order to properly represent the operating stationary conditions of the vacuum chamber during the experiments (see Section 2). No neutral recycling from the chamber walls is considered, since the relative reduction of the background neutrals density (due to CEX and ionization collisions) is negligible throughout the simulation duration. When ion/neutral macro-particles exit from the lateral or downstream boundaries, they are simply removed from the simulation. In order to respect the CFL condition, the time step is 0.25 µs, while the overall simulation time is 1.25 ms, sufficiently large to have slow neutral particles cross the whole domain. Finally, a total of around 45 million macro-particles are simulated in stationary conditions.

Coming now to the electron fluid model, this is the same as the one described in Ref. [25], but with a different numerical discretization strategy [31]. The model assumes a quasineutral plasma everywhere, i.e. ne = ni, a polytropic and isotropic electron thermodynamics, Teneγ1, with γ representing the electron polytropic coefficient, and a tensor electron conductivity K. For a quasineutral plasma, the total electric current density j = je + ji satisfies the continuity equation with no source terms:


Then, assuming both isotropic and polytropic electrons, the generalized Ohm’s law can be derived from the electron momentum balance equation as [25]:


where K is the normalized conductivity tensor (defined below), σe = e2ne/(meνe) is the electron scalar conductivity, jc is an effective current density grouping collisional effects with heavy species [25], and a “residual thermalized potential” Φ has been introduced. The latter is defined so that ∇Φ = ∇ϕ − ∇he/e, being ϕ the electric potential and ∇he = ∇pe/ne the barotropic function gradient (which is an exact differential when polytropic electrons are considered [25]). The effective current density is defined as jc=ene/νes=1Lνesus, where the summation extends over the L considered heavy species, and νe=s=1Lνes+νan is the total electron momentum transfer collision frequency, which also includes an anomalous collision frequency νan. The models for νes are described in Ref. [23] and include the effect of both ionization collisions, and elastic collisions with neutrals and ions, as shown in Table 1. The anomalous collision frequency νan follows a Bohm’s like transport [32] and is computed as νan = αanωce, with ωce = eB/me the electron gyro-frequency. Here, αan = 2%, which permits to limit the effective Hall parameter χ = ωce/νe to a maximum value of 50, thus improving the solver’s convergence, although it has been verified that the physical solution is only dimly affected by the chosen value. Finally, K is


where 1b=bx,by,bz is the unit vector along the applied magnetic field, so that B = B1b.

The gradient of the thermalized potential Φ in Eq. 2 measures the correction to be applied to the Boltzmann’s electric field due to the effects of both magnetic field and collisions on the electron fluid [25]. Assuming the same reference point for Φ, he and ϕ, we have


hence, the electric potential can be retrieved from the knowledge of he(ne) (known from the PIC) and of the unknown Φ. For an electron polytropic coefficient γ, the barotropic function depends only on the quasineutral plasma density ne [33]:


with ne0, Te0 the quasineutral plasma density and electron temperature at the reference point for potential.

The numerical approach for the discretization of Eqs 1, 2 is a 1st order hybrid Finite Volume (FV)–Finite Differences (FD) scheme on a staggered mesh [31]. In particular, Eq. 1 is discretized with FV, while Eq. 2 with FD schemes. Figure 4 shows a 2D representation of the unknowns locations, which are the electric current densities at cell face centers and the thermalized potentials at cell centers, including ghost cells at the boundaries.


FIGURE 4. Sketch of a 2D xy cross section of the deformed staggered mesh, showing the locations of the unknowns. The solid bold black line represents the external boundary, black dots the mesh nodes, blue circles the cell centers at which Φ is solved for, and red arrows the normal components of j, solved at cell faces centers.

This scheme is conservative due to the Finite Volume discretization of the continuity equation with a negligible global continuity error, so that there is no need for a very fine mesh as considered in Refs. [25,26] and the computational cost is significantly lower. In addition, FD schemes are centered everywhere (no need for forward/backward schemes at the boundaries), which leads to a lower discretization error. After the discretization on the staggered mesh, Ohm’s law, Eq. 2, is substituted into the continuity equation, Eq. 1, leading to a system of equations in the unknown thermalized potential:


with [M] a square matrix with maximum rank obtained by imposing the Dirichlet condition on Φ at a reference node (here x = y = z = 0), and {R} a right hand side vector.

Finally, boundary equations are explicitly imposed on the normal electric current densities at all external boundary cell faces. In particular, a local current free condition is applied: j1n = 0 (see Figure 2), where 1n represents the unit vector normal to the local boundary and pointing towards the plasma domain, as shown in Figure 4.

3.2 RPA Synthetic Simulations Model and Setup

The RPA from Impedans [30] has been briefly introduced in Section 2. Here a simplified domain is considered, which covers only the region downstream of the grounded orifice grid, including a single hole of the electron repelling, retarding potential and secondary electron emission (SEE) suppression grids, plus the collector plate at the downstream boundary, as shown in Figure 5. Notice that the coordinates origin O′ is now at the location of the RPA, and not at the magnetic nozzle throat, as in the magnetic nozzle simulation of Section 3.1.


FIGURE 5. (A) Schematic view of the simulation setup in the x′ − z′ plane (the y′ − z′ cross section is analogous) and (B) evolution of the electric potential (relative to the local plasma) along the RPA centerline. In subplot (A), electrons are present only in the light red region, before the electron repelling grid, while they are absent in the rest of the simulation domain.

The parameters considered for the RPA simulation are then reported in Table 2. The single simulated hole is a square with 20 µm side (an RPA is constituted by thousands of such holes), the thickness of the grids is 0.05 mm and the distance between them is 0.2 mm [30]. The considered mesh is a parallelepiped mesh with 15 × 15 × 93 nodes and a physical extension along x′, y′, z′ equal to 28 μm × 28 μm x 2.3 mm. This domain extension along x and y yields a grid open area equal to 50% (also known as transmission T), as specified by the RPA manufacturer, so that Ts=Tr=Tsee=0.5.


TABLE 2. Simulation parameters for the magnetic nozzle simulation. Grids potentials are relative to the local plasma plume potential ϕp, just outside of the RPA.

Regarding the PIC model, only singly charged ions are simulated and are injected from the upstream boundary of Figure 5 with a generic distribution function, obtained from the conditions at the RPA N.1 location of the magnetic nozzle simulation (see Figure 2). When ions hit any material surface (of either the grids or the collector) they are simply removed from the simulation, while if they cross the lateral boundaries, they are reflected periodically, meaning that they are moved to the opposite side of the simulation domain, with unaltered velocities. This allows us to simulate the interaction between different holes inside the RPA. For the chosen mesh spacing, the CFL condition now requires a PIC time step equal to 2 ns. At stationary conditions, around 3 million macro-particles are used.

Electrons are simulated as an isothermal fluid, with a temperature Te0 = 1.5 eV, i.e., the electron temperature obtained in the magnetic nozzle simulation at the RPA N.1 location. In order to properly simulate the charge density, electrons are assumed to be present only in the region upstream of the electron repelling grid (see Figure 5), with a density given by the Boltzmann’s relation:


with ne0 the electron density at x′ = y′ = z′ = 0 where the plasma is assumed quasineutral and the electric potential is equal to ϕp, the undisturbed plasma potential just outside of the RPA. Neither magnetic field nor collisions effects are thus retained in this electron fluid model. This is justified by the fact that at the RPA the plasma is weakly collisional and at the same time, the thruster magnetic field is vanishing. As fully described in Ref. [23], the electric potential is solved through a non-linear Poisson’s equation:


where ni is known from the PIC, and ne is a function of the electric potential through Eq. 7 (hence the non-linear character of the equation). The boundary conditions are either homogeneous Neumann conditions on the lateral and upstream boundaries (∂ϕ/1n = 0) or Dirichlet conditions at the grids and at the collector (the assumed potentials are reported in Table 2). In particular, the retarding grid potential is varied within the range [-0.1,62.5] V, through 43 consecutive simulations with a duration of 3 µs each (therefore during this time all potentials are maintained fixed).

The ion current hitting the collector can thus be obtained as a function of ϕr, and from this, the axial ion velocity distribution function (IVDF) as a function of the axial kinetic energy Kz=1/2mivz2=e(ϕrϕp) is obtained as [34]:


where Ac is the collector area and T is the total system transmission. The latter can be roughly computed as the product of the transmissions of the grids, T=TsTrTsee=, provided that the holes diameter (20 µm) is much smaller than the distance between the grids (200 µm). In fact, if the distance goes to zero and the holes are aligned, the global transmission should tend to the electron repelling grid transmission alone (0.5). Notice that Eq. 9 provides the axial IVDF as a function of the kinetic energy (and not of the axial velocity) and it is often referred to as the ion energy distribution function (IEDF), although, if integrated in energy, it provides the incoming ion flux (instead of the density): Kz=0f(Kz)dKz/mi=niuzi.

4 Magnetic Nozzle Simulation Results

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.


FIGURE 6. Simulation results at y = 0: (A) magnetic induction field intensity and direction, (B) Hall parameter.

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.


FIGURE 7. Simulation results at y = 0: (A) electric potential, (B) electric field, (C) electron density, (D) ratio between source ion density and plasma density, (E) electron current density magnitude (with mean flow direction), and (F) ion current density magnitude. Arrows in subplots (B), (E), and (F) indicate the direction of the vectors and their size is proportional to the ratio between their xz component and their total magnitude.


FIGURE 8. Simulation results at z = 0, showing (A) electric potential, (B) electric field with direction, (C) electron current density magnitude with mean flow direction, and (D) ion current density magnitude with direction. The size of the arrows in subplots (B), (C), and (D) is proportional to the ratio between the xy component of the vector and its total magnitude.

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 1019 m−3 at the throat to around 7⋅1015 m−3 downstream at the centerline.

In order to evaluate the effect of collisions in the expansion, the ratio between the number density of source ions ni, src (i.e., ions injected from z = 0) and the total plasma density ni = ne 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 ×BBrjθ1z, 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. je+ji1n=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 (je ×BdV, with dV the differential volume), and 0.1 mN of magnetic force on ions (ji ×BdV). As expected, the ion magnetic force contribution is much smaller than that of the electrons.

5 RPA Synthetic Simulation Results

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


FIGURE 9. Simulation results at y′ = 0 for (A,D,G) ϕr = ϕp, (B,E,H) ϕr = ϕp + 25.5 V, (C,F,I) ϕr = ϕp + 49.5 V. Subplots (A,B,C) show the electric potential, subplots (D,E,F) the ion density, and subplots (G,H,I) the 3D ion trajectories projected into the y′ = 0 plane. The end trajectory point is shown by an empty circle. Red trajectories correspond to collected particles, black solid lines to particles that hit the grids, and black dashed lines to particles that are reflected backwards out of the RPA. Dotted lines are used to highlight periodic reflections.

At the densities and electron temperature of the local plasma plume (respectively less than 1016 m−3 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.


FIGURE 10. (A) I-V characteristic obtained from the numerical RPA simulations, and (B) evolutions of the ion currents collected by the electron repelling grid (green line), retarding grid (red line), and SEE suppression grid (blue line). Each circle corresponds to a different simulated retarding grid potential.

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.

6 Comparison Between Simulations and Experiments

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


FIGURE 11. Comparison between magnetic nozzle simulations (black solid line) and experiments (green dots) in terms of (A) plasma density along the centerline, and (B) ion current density versus the angle θ at R = 40 cm distance from the nozzle throat.

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


FIGURE 12. Comparison between magnetic nozzle simulations (red solid line refers to slow ions, blue solid line to injected ions and black solid lines to the sum of the two populations) and experiments (green circles) of the IVDF as a function of the kinetic energy KR along the R direction, at R = 40 cm, for (A) θ = 0°, (B) θ = 5° and (C) θ = 10°.

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 45eV, 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 1520 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.


FIGURE 13. Comparison between magnetic nozzle simulations (black solid lines), experiments (green circles) and RPA simulations (blue dotted and grey dashed lines) of the IVDF along the z direction (θ = 0). Regarding RPA simulations, the grey dashed line with circle markers is obtained from differentiation of the simulated collector current Ic, while the dotted blue line with triangle markers from the differentiation of the total current Ic′ reaching at least the retarding grid: Ic′ = Ic + Ir + Isee.

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 ∝ − dIc/dϕr, while the second one (blue dotted line with triangle markers) is a corrected IVDF computed by applying Eq. 9 not to Ic alone but to the sum of the collector, SEE suppression and retarding grids currents: Ic = Ic + Ir + Isee, so that f ∝ − dIc/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.

7 Conclusions

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.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

FC prepared, ran and postprocessed the simulations described in this work. JN carried out the totality of the experimental characterization campaign of the Helicon plasma plume. GR and FC implemented the RPA simulation model. AM implemented the new solution algorithm for the magnetized plume model. All authors contributed in the writing of the paper text. Part of FC’s contribution was completed while he was affiliated with Equipo de Propulsión Espacial y Plasmas (EP2), Universidad Carlos III de Madrid.


This work has been funded by European Union’s Horizon 2020 research and innovation program under grant agreement No 870542 HIPATIA (HelIcon PlasmA Thruster for In-space Applications), and the Spanish National research plan under the project ESPEOS, PID2019-108034RB-I00/AEI/10.13039/501100011033.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.


1. Takahashi K, Charles C, Boswell RW, Takao Y, Fruchtman A, Navarro-Cavallé J, et al. Commentary: On Helicon Thrusters: Will They Ever Fly? Front Phys (2020) 8:277. doi:10.3389/fphy.2020.00277

CrossRef Full Text | Google Scholar

2. Navarro-Cavallé J, Wijnen M, Fajardo P, Ahedo E. Experimental Characterization of a 1 Kw Helicon Plasma Thruster. Vacuum (2018) 149:69–73.

Google Scholar

3. Trezzolani F, Manente M, Toson E, Magarotto M, Moretto D, Pavarin D, et al. Development and Testing of a Miniature Helicon Plasma Thruster. In: IEPC-2017-519: 35th International Electric Propulsion Conference. Atlanta, GA: Electric Rocket Propulsion Society (2017).

Google Scholar

4. Packan D, Elias P, Jarrige J, Merino M, Sánchez-Villar A, Ahedo E, et al. The MINOTOR H2020 Project for ECR Thruster Development. In: IEPC-2017-547: 35th International Electric Propulsion Conference. Atlanta, GA: Electric Rocket Propulsion Society (2017).

Google Scholar

5. Inchingolo MR, Merino M, Navarro-Cavallé J. Hybrid Pic-Fluid Simulation of a Waveguide Ecr Magnetic Nozzle Plasma Thruster. In: Space Propulsion Conference 2021 (March 17-19: Association Aéronautique et Astronautique de France) (2021). p. 00192.

Google Scholar

6. Takahashi K, Takao Y, Ando A. Performance Improvement of a Magnetic Nozzle Plasma Thruster. In: 36th International Electric Propulsion Conference (2019).

Google Scholar

7. Lai ST, Miller C. Retarding Potential Analyzer: Principles, Designs, and Space Applications. AIP Adv (2020) 10:095324. doi:10.1063/5.0014266

CrossRef Full Text | Google Scholar

8. Brown DL, Walker ML, Szabo J, Huang W, Foster JE. Recommended Practice for Use of Faraday Probes in Electric Propulsion Testing. J Propulsion Power (2016) 33:582–613.

Google Scholar

9. Lobbia RB, Beal BE. Recommended Practice for Use of Langmuir Probes in Electric Propulsion Testing. J Propulsion Power (2017) 33:566–81. doi:10.2514/1.B35531

CrossRef Full Text | Google Scholar

10. Sudit ID, Chen FF. Rf Compensated Probes for High-Density Discharges. Plasma Sourc Sci. Technol. (1994) 3:162–8. doi:10.1088/0963-0252/3/2/006

CrossRef Full Text | Google Scholar

11. Sheehan JP, Raitses Y, Hershkowitz N, McDonald M. Recommended Practice for Use of Emissive Probes in Electric Propulsion Testing. J Propulsion Power (2017) 33:614–37. doi:10.2514/1.b35697

CrossRef Full Text | Google Scholar

12. MacDonald NA, Cappelli MA, Gildea SR, Martínez-Sánchez M, Hargus WA. Laser-induced Fluorescence Velocity Measurements of a Diverging Cusped-Field Thruster. J Phys D: Appl Phys (2011) 44:295203. doi:10.1088/0022-3727/44/29/295203

CrossRef Full Text | Google Scholar

13. Knudsen WC, Spenner K, Bakke J, Novak V. Pioneer venus Orbiter Planar Retarding Potential Analyzer Plasma experiment. IEEE Trans Geosci Remote Sensing (1980) GE-18:54–9. doi:10.1109/tgrs.1980.350261

CrossRef Full Text | Google Scholar

14. Polansky J, Wang J, Ding N. Experimental Investigation on Plasma Plume Potential. IEEE Trans Plasma Sci (2013) 41:3438–47. doi:10.1109/tps.2013.2277724

CrossRef Full Text | Google Scholar

15. Lemmer KM, Gallimore AD, Smith TB, Austin DR. Review of Two Retarding Potential Analyzers for Use in High Density Helicon Plasma. In: 30th International Electric Propulsion Conference (2007).

Google Scholar

16. Davidson RL, Earle GD. A Design Approach for Improving the Performance of Single-Grid Planar Retarding Potential Analyzers. Phys Plasmas (2011) 18:012905. doi:10.1063/1.3533657

CrossRef Full Text | Google Scholar

17. Enloe CL, Shell JR. Optimizing the Energy Resolution of Planar Retarding Potential Analyzers. Rev scientific Instr (1992) 63:1788–91. doi:10.1063/1.1143339

CrossRef Full Text | Google Scholar

18. Martínez-Sánchez M, Navarro-Cavallé J, Ahedo E. Electron Cooling and Finite Potential Drop in a Magnetized Plasma Expansion. Phys Plasmas (2015) 22:053501.

Google Scholar

19. Kojima T, Morita T, Yamamoto N. Analysis of Plasma Detachment in the Magnetic Thrust Chamber Using Full Particle-In-Cell Simulation. High Energ Density Phys (2020) 36:100814. doi:10.1016/j.hedp.2020.100814

CrossRef Full Text | Google Scholar

20. Andrews S, Di Fede S, Magarotto M. Fully Kinetic Model of Plasma Expansion in a Magnetic Nozzle. Plasma Sourc Sci Technol (2022). doi:10.1088/1361-6595/ac56ec

CrossRef Full Text | Google Scholar

21. Bird G. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. The Oxford Engineering Science Series. Oxford, UK: Oxford University Press (1994).

Google Scholar

22. Hockney R, Eastwood J. Computer Simulation Using Particles. Boca Ratón, FL: CRC Press (1988).

Google Scholar

23. Cichocki F, Domínguez-Vázquez A, Merino M, Ahedo E. Hybrid 3D Model for the Interaction of Plasma Thruster Plumes with Nearby Objects. Plasma Sourc Sci. Technol. (2017) 26:125008. doi:10.1088/1361-6595/aa986e

CrossRef Full Text | Google Scholar

24. Cichocki F, Merino M, Ahedo E. Spacecraft-plasma-debris Interaction in an Ion Beam shepherd mission. Acta Astronautica (2018) 146:216–27. doi:10.1016/j.actaastro.2018.02.030

CrossRef Full Text | Google Scholar

25. Cichocki F, Merino M, Ahedo E. Three-dimensional Geomagnetic Field Effects on a Plasma Thruster Plume Expansion. Acta Astronautica (2020) 175:190–203. doi:10.1016/j.actaastro.2020.05.019

CrossRef Full Text | Google Scholar

26. Cichocki F, Domínguez-Vázquez A, Merino M, Fajardo P, Ahedo E. Three-dimensional Neutralizer Effects on a Hall-Effect Thruster Near Plume. Acta Astronautica (2021) 187:498–510. doi:10.1016/j.actaastro.2021.06.042

CrossRef Full Text | Google Scholar

27. Perales-Díaz J, Cichocki F, Merino M, Ahedo E. Formation and Neutralization of Electric Charge and Current of an Ion Thruster Plume. Plasma Sourc Sci. Technol. (2021) 30:105023. doi:10.1088/1361-6595/ac2a19

CrossRef Full Text | Google Scholar

28. Navarro-Cavallé J, Wijnen M, Fajardo P, Merino M, Ahedo E. Experimental Performances of a 1 Kw Hpt by Means of Plasma Diagnostics. In: IEPC-2017-447: 35th International Electric Propulsion Conference. Atlanta, GA: Electric Rocket Propulsion Society (2017).

Google Scholar

29. Navarro-Cavallé J, Wijnen M, Fajardo P, Ahedo E, Gómez V, Giménez A, et al. Development and Characterization of the Helicon Plasma Thruster Prototype Hpt05m. In: IEPC-2019-596: 36th International Electric Propulsion Conference. Vienna, Austria: Electric Rocket Propulsion Society (2019).

Google Scholar

30. Gahan D, Dolinaj B, Hopkins MB. Retarding Field Analyzer for Ion Energy Distribution Measurements at a Radio-Frequency Biased Electrode. Rev Scientific Instr (2008) 79:033502. doi:10.1063/1.2890100

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Leveque RJ. Finite Volume Methods for Hyperbolic Problems. In Cambridge Texts in Applied Mathematics. Cambridge: Cambridge University Press Vol. 31 (2002). doi:10.1017/CBO9780511791253

CrossRef Full Text | Google Scholar

32. Taccogna F, Garrigues L. Latest Progress in Hall Thrusters Plasma Modelling. Rev Mod Plasma Phys (2019) 3:12. doi:10.1007/s41614-019-0033-1

CrossRef Full Text | Google Scholar

33. Merino M, Cichocki F, Ahedo E. A Collisionless Plasma Thruster Plume Expansion Model. Plasma Sourc Sci. Technol. (2015) 24:035006. doi:10.1088/0963-0252/24/3/035006

CrossRef Full Text | Google Scholar

34. Böhm C, Perrin J. Retarding-field Analyzer for Measurements of Ion Energy Distributions and Secondary Electron Emission Coefficients in Low-Pressure Radio Frequency Discharges. Rev scientific Instr (1993) 64:31–44.

Google Scholar

35. Wijnen M, Navarro-Cavallé J, Fajardo P. Mechanically Amplified Milli-newton Thrust Balance for Direct Thrust Measurements of Electric Thrusters for Space Propulsion. IEEE Trans Instrumentation Meas (2021) 70:3505318. doi:10.1109/tim.2020.3037305

CrossRef Full Text | Google Scholar

Keywords: helicon, plasma plumes, retarding potential analyzer (RPA), particle in cell (PIC), fluid models, hybrid models, magnetic nozzle (MN)

Citation: Cichocki F, Navarro-Cavallé J, Modesti A and Ramírez Vázquez G (2022) Magnetic Nozzle and RPA Simulations vs. Experiments for a Helicon Plasma Thruster Plume. Front. Phys. 10:876684. doi: 10.3389/fphy.2022.876684

Received: 15 February 2022; Accepted: 11 March 2022;
Published: 14 April 2022.

Edited by:

Vladimir I. Kolobov, CFD Research Corporation (United States), United States

Reviewed by:

Kazunori Takahashi, Tohoku University, Japan
Ralf Schneider, University of Greifswald, Germany

Copyright © 2022 Cichocki, Navarro-Cavallé, Modesti and Ramírez Vázquez. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Filippo Cichocki,