Numerical Study of Electron Cyclotron Drift Instability: Application to Hall Thruster

We study the cross-field electron transport in Hall thruster induced by ExB cyclotron drift instability. The investigation tool, consisting of one-dimensional Particle-in-Cell model in the azimuthal drift direction, has been subjected to a convergence study to verify the effects of numerical parameters. The instability evolves keeping the discrete nature of its cyclotron harmonics only for low wavenumbers. A resonance broadening mechanism induced by the electron heating and velocity distribution deformation makes the high-wave numbers disappear in the long non-linear stage. A large wavelength modulation, comparable to the entire azimuthal domain considered is always superimposed. The saturation mechanism is conducted by ions that, due to friction with electrons and trapping in the potential well, heat up and rotate along the electron drift direction. With the best numerical parameters found, the scaling of anomalous mobility with the most important physical quantities (gas and plasma density, magnetic field, accelerating axial electric field and ion mass) has been obtained. Results show that the electron cross-field mobility has a strong dependence from the plasma density and ion mass: for large plasma density, the system undergoes an abrupt change entering in a mode dominated by fluctuation-induced transport, while lighter ions present larger mobility. The scaling of the dominant wavelength detected is compatible with the first harmonic and no transition towards ion acoustic like instability has been observed.


INTRODUCTION
Hall thrusters (HTs) [1][2][3][4] are ExB plasma-based devices used for Space application. In the most common coaxial version, electrons are emitted from an external cathode (working also as ion plume neutralizer), while the anode neutral propellant injector is fixed upstream at the beginning of the channel. A radial (the cylindrical reference system will be used in the rest of the paper) magnetic field lowers the electron mobility along the axial direction leading to a double effect: (1) generate an axial electric field in proximity of the exit plane that accelerates the un-magnetized propellant ions (usually noble gases) without the necessity of any extraction grid; therefore, HTs do not suffer from the limitation of the Child-Langmuir law and are characterized by a relatively high acceleration efficiency; (2) generate an azimuthal closed ExB drift v de and increase the electron residence time inside the channel in order to guarantee an high level of ionization.
Nevertheless, the measured electron cross-field mobility results hundreds times larger [4] than its classical collisionalinduced value Here, m e and q are the electron mass and charge, while ν m and ω ce are the electron collisional and cyclotron frequencies, respectively. Despite numerous studies on different physical mechanisms involved in HTs, the problem of electron anomalous transport across the magnetic field remains an open issue. Among the different mechanisms proposed [5][6][7][8][9] to explain the electron anomalous transport, the most plausible candidate is the electron ExB drift instability [10,11] generated by the unbalanced electron Hall current. It results in azimuthal modulations of electron/ion density and azimuthal electric field whose frequency and wavelength of the order of MHz and mm, respectively, have been observed in different experiments [12]. The correlation between electron density n e and azimuthal electric field E θ , equivalent to a friction force density [10,11] f ei = q n e E θ (2) between the rotating electrons and almost immobile ions along the azimuthal direction, gives an extra contribution to the electron cross-field mobility (derived by electron momentum conserving equation taking into account azimuthal instabilities [10]) where n e0 and E z0 represent the azimuthally averaged electron density and accelerating axial field, respectively. This instability has been analyzed in detail by different 1D-azimuthal Particle-in-Cell (PIC) [13,14] models [15][16][17][18][19][20] assuming fixed axial electric and radial magnetic fields in Cartesian geometry (often used due to the high value of the external diameter to channel with ratio) and using a periodicity length along the azimuthal coordinate always smaller than the real one ( Figure 1 represents the corresponding simulation domain).
In these models steady state solutions are reached only if particle losses along the virtual axial boundaries are included. In fact, particles are still tracked along the axial direction (even if their position is not used to self-consistently calculate the axial electric field) and each electron (resp. ion) leaving the simulation domain on the anode (resp. cathode) side is replaced by a new electron (resp. ion) at the cathode (resp. anode) side loaded with an half-Maxwellian velocity distribution with the initial temperature. This particle re-injection mechanism, necessary to avoid undesirable secular effect and to mimic the axial transport, works like an extra collision: reducing the axial-replacement length L z,acc is equivalent to increase the collisional frequency in a 1D-azimuthal model. Usually, the value chosen corresponds FIGURE 1 | Sketch of the simulation domain used for 1D-azimuthal PIC simulation of HTs. L y is the size of the azimuthal domain, while L z,acc represent the distance traveled by electrons and ions before to be re-injected from the opposite side.
to the acceleration region, where both the axial electric and radial magnetic fields can be considered uniform. Probably, no replacement can be still used only if one is interested to just characterize the linear stage of the instability occurring in <1 µs. However, in order to represent the following non-linear evolution and saturation, an axial re-placement is essential.
All these studies have confirmed a significant fluctuation amplitude of ExB drift instability characterized during the linear stage by high-frequency (of the order of MHz) and discrete number of short-wavelengths (the dominant one of the order of mm) corresponding to cyclotron harmonics where m is an integer. An important feature observed is the difference in the ion and electron density fluctuations. While the electron density perturbation has lower amplitude sinusoidal character, the ion density perturbation is characterized by much larger amplitudes and contains shorter wavelengths (nonquasineutral) with λ ≤ λ De (λ De being the electron Debye length) producing highly-localized cnoidal structures with sharper crests and flatter troughs more. The saturation mechanisms of the instability have been identified in the electron heating, deformation of the electron distribution function projected along the azimuthal direction and ion-wave trapping. The non-linear stage evolution of the instability instead has given controversial results: some models [17] show a transition of the mode into the regime of the ionsound instability typical of un-magnetized plasmas, while other models [18] suggest that the magnetic field remains important in the mechanism of the instability, electron heating, and transport. This discrepancy can be attributed to the intrinsic limitation of missing the axial coordinate. In fact, in [17] a re-injection criterion is used while in [18] the particles are kept throughout the simulation as with an infinite acceleration length.
In all the models, the computed cross-field mobility has very large values (one order of magnitude larger) compared with that measured in HTs. The last point can be related to the fact that the acceleration field E z0 is externally imposed and kept fixed, while one expects that an increase of the local mobility would yield a decrease of the local electric field, which would reduce the drift velocity and control the possible instabilities leading to the anomalous electron transport.
Recent models [20][21][22][23][24][25] have added the radial coordinate resolving the direction along the magnetic field and accounting for the sheath boundary conditions. They have shown that the full ion sound regime is not realized, but a more effective regime of discrete cyclotron resonance instabilities occur. The magnetic field remains to be a defining feature of ExB drift instability for HT operational regime. On the contrary, the scaling of wavelength, frequency and amplitude of the azimuthal field obtained within azimuthal-axial models [19,[26][27][28][29][30] is consistent with the ion acoustic like character of the instability (ω pi and T e are ion plasma frequency and electron temperature, respectively). All the azimuthal models show that results are very sensitive to the numerical parameters used. Instabilities require a very good resolution both in space/time and particle number, in order to have the modes well represented and running up to nonlinear and saturation regimes. In particular, numerical noise is responsible for the destruction of electron cyclotron resonance and facilitates the transition to an ion acoustic instability. For this reason, in this work, we have preferred to develop a 1Dazimuthal model by using a high-statistical level. In the next section, the model is explained in detail and the effect of the different re-injection methods is presented. Section Reference Case and Effect of Periodicity Length is devoted to the effect of the azimuthal length of the simulation domain on the evolution of the wave spectrum. The scaling of the mobility as a function of plasma density, neutral density, magnetic field, electric field, and ion mass has been conducted in section Parametrical Study and Scaling Laws. Finally, the conclusions have been drawn in the last section Conclusions.

DESCRIPTION OF THE MODEL AND NUMERICAL CONVERGENCE STUDY
The numerical method used is based on electrostatic 1Dazimuthal Particle in Cell / Monte Carlo Collision (PIC-MCC) methodology [13,14]. Uniform radial magnetic B 0 and axial electric field E z0 are kept fixed during the simulation that starts loading the particles (electrons and ions) uniformly along the azimuthal and axial directions of the simulation domain (represented in Figure 1) and with a prescribed velocity distribution. Electrons and ions are both initialized with a shifted-Maxwellian distribution: electrons at a temperature T e0 and with an azimuthal drift v de = E z0 /B 0 , ions with a temperature T i0 and an axial drift v id = 2qE z0 L z,acc /M (M is the ion mass) corresponding to the velocity acquired by the ions in the acceleration region. Differently from the other components, the axial component of ion velocity is not updated during the simulation and considered as frozen (approximation well-posed since ions are treated as collisionless). All the physical parameters used are summarized in Table 1.
The cartesian approximation is adopted neglecting the effect of the presence of the centrifugal force. From now on we will use y to represent the azimuthal coordinate θ . The time-step used is t = 2x10 −12 s. The self-consistent electric field is solved along the ExB direction with a Thomas tridiagonal algorithm modified to implement periodic boundary conditions [31]. The reference zero potential value is put at the boundaries.
Electron-neutral collisions (elastic scattering, electronic excitation, and ionization) are included. In order to keep the averaged plasma density equal to the initial condition selected, we assume that during the ionization event, only the energy loss and scattering is taken into account, while then new electron-ion pair is not produced.
In the present model, we consider the particle re-injection across the axial boundaries. Differently from [17] the re-injected particles are not randomly distributed along the azimuthal direction but they keep their azimuthal positions. They just change their velocity components, exactly like in an inelastic collisional event. This choice is dictated by the fact that, with a random y-location re-injection, a large charge imbalance is artificially created amplified by the fact that the model is one-dimensional. A certain space-charge density (difference between electrons reaching the anode side and ions crossing the cathode one) is created along the ExB direction that generates Hall parameter β 2.70 × 10 3 an artificially high azimuthal electric field. In the 2D-radialazimuthal extension of the model [25], the random re-injection has much less impact due to the fact that the charge injected is also distributed along the radial dimension. Furthermore, we have verified that results obtained with random azimuthal reinjection are sensitive to the couple of parameters (L y , L z,acc ) chosen, as already observed in [17]. For small L y (≤1 cm) or large L z,acc (>2 cm), the reference potential or the low rate of particles re-injected help to keep the azimuthal electric field low. Figure 2 shows the temporal evolution of the averaged (over the particle population) electron (blue curve) and ion (red curve) energy when re-injected particles are loaded with random yposition (Figure 2A) and when they keep their own azimuthal position had just before the re-injection ( Figure 2B) using L y = 2.67 cm and L z,acc = 1 cm. In the first case, both electron and ion energy increase reaching already after 3 µs unphysical values due to the large azimuthal electric field induced by the random re-injection.
Due to the high sensitivity shown by previous models on the numerical parameters used, we have studied the numerical convergence. The global quantities taken as monitoring observables are the electron cross-field mobility and electron density fluctuation level. In case of 1D-azimuthal PIC representation, the mobility can be simply extracted from the averaged axial velocity where v ip,ez is the axial component of the electron velocity and N p is the total number of particles used in the simulation. The physical condition used are that summarized in Table 1. Figure 3A shows the electron cross-field mobility and Figure 3B shows the amplitude of fluctuations in electron density computed as a function of number of particles per cell N ppc (from 20 to 1,000) for four different numbers of mesh nodes, N g = 250, 500, 1,000, and 2,000, corresponding to four different values of cell size to Debye length ratio y λ D = 1.44, 0.72, 0.36, and 0.18. Results correspond to a time larger than 10 µs and the mobility and fluctuation amplitude reported were also time averaged over the last 2 µs. In addition, the time series of electron averaged energy for different N ppc with N g = 1,000 and for different N g with N ppc = 500 have also been reported in Figures 3C,D, respectively. For N ppc >500 all the cases give results close to the averaged value of µ PIC = 3.4 m 2 V −1 s −1 within a deviation <5% and fluctuations in the electron density around 22% of the equilibrium density. The electron heating growth rate does not change if using N ppc larger than 500 and N g larger than 1,000. Therefore, for the rest of the work, we have fixed the numerical parameters at N ppc = 500 and y λ D = 0.36.

REFERENCE CASE AND EFFECT OF PERIODICITY LENGTH
In this section we present results obtained using physical data reported in Table 1. It is considered as reference case for the physical parametrical scan study reported in the next section Parametrical study and scaling laws. Figure 4 shows an overview of the simulation run. The history of the global averaged electron / ion energy and electrostatic potential energy density u E = 1 2 ε 0 E 2 y are displayed over a total time of 5 µs. In the graph one can identify four stages: (1) the linear growth, a very fast process terminating already after the first microsecond; (2) the first non-linear stage ending at t = 3.5 µs and during which the electron energy appears to saturate; this phase corresponds to a redistribution of the spectral energy toward lower harmonics (see Figure 6B) and to the ion heating; (3) the second non-linear phase when the electron and ion energy rise again; (4) the final quasi-steady state reached after t = 4.5 µs. Figure 5 show the azimuthal profiles of Figure 5A electron (red curve) and ion (blue curve) density and Figure 5B oscillating electric field at t = 12 µs, after the saturation, in the non-linear regime. Profiles confirm large amplitude coherent mode driven at the fundamental cyclotron harmonic λ 1 = 1.72 mm [Equation  =1.44 (number of mesh nodes N g = 250), 0.72 (N g = 500), 0.36 (N g = 1,000) and 0.18 (N g = 2,000). Electron averaged energy evolution for different N ppc using N g = 1,000 (C) and for different N g using N ppc = 500 (D). All the cases correspond to physical parameters reported in Table 1 using azimuthal length of L y = 2.67 cm.
FIGURE 4 | Electron (red curve), ion (blue curve) averaged energy and electrostatic potential energy density (black curve) in the simulation up to t = 5 µs. Results refer to the reference case using physical parameters reported in Table 1.
(4) with m = 1 gives λ ExB−EDI,m = 1 = 1.76 mm] with a marked cnoidal character for the ion density modulation, effect of the localized ion trapping (see end of this section), together with the presence of a large scale envelope. Figure 6A shows the time evolution of the azimuthal profile of E y and Figure 6B shows the time evolution of the corresponding power spectral density (PSD) into wavenumber k y = 2π/λ components. The instability starts during the linear stage (for t<2 µs) with the growth of the most instable cyclotron harmonics (seven different harmonics are easily distinguishable from the spectrum) whose dominant one has the following characteristics: wavelength k = 3,650 m −1 [corresponding to the first cyclotron harmonic m = 1 of Equation (4)], frequency ω = 4.8 MHz and phase velocity v φ = 8 km/s in the ExB drift direction. However, the values found are also in good agreement with Equation (5), the theoretical estimation of an ion acoustic character found in Lafleur et al. [11] using the wavenumber with the maximum growth. Using an electron temperature of T e = 60 eV (computed in the non-linear phase, see Figure 7C) Equation (5) gives k IAI = 3,890 m −1 , while Equation (6) gives ω = 3.2 MHz. In the next section Parametrical study and scaling laws, a dedicated study is conducted to distinguish the true character of instability.  Table 1. For t>2 µs, the highest cyclotronic harmonics are strongly attenuated. The resonance character of the ExB EDI based on the linear treatment and Maxwellian behavior of electron distribution function, disappears by a resonance broadening mechanism [32]. Only the first harmonic remains detectable making the electrons still keeping their magnetization even if much weaker: electrons trajectories are deformed by the azimuthal wave of electric field and slightly elongated in the direction parallel to E z0 (see Figure 7A). A simple criterion for resonance destruction [18] (corresponding to impose an electron displacement z e larger than the half-wavelength over one cyclotron period): highlights the dependence on the electron temperature for the different k-modes: the highest wavenumbers are faster damped requiring a smaller electron temperature threshold. Electrons can also move to the right, but in average they diffuse toward the left as evident from the asymmetry of the probability distribution function of the electron guiding centers jump size z e,gc reported in Figure 7B.

Figure 7C
shows the origin of the resonance broadening, i.e., the deformation of the electron velocity distribution function: its projection along the azimuthal direction at the end of the simulation (red curve) is compared with the initial Maxwellian distribution at T e0 = 10 eV (blue curve) and with its Maxwellian best fitting (green curve). The latter shows that electrons heat up to a level of T e = 60 eV due to their continuous trapping and detrapping induced by the bouncing into the potential well-created by wave field. The electron heating create a strong phase-space mixing of the electron distribution function leading to a marked low-energy range depletion and electron velocity distribution function (EVDF) flattening.
Furthermore, a very low k-number component (much lower than the fundamental cyclotron mode k 1 ) appears corresponding to wavelengths comparable to the entire simulation domain size (k y L y /2π = 1 in Figure 6B). Already previous 1D-azimuthal and 2D-radial-azimuthal models [18,23,25] have shown non-linear development of large-scale modes, comparable to the entire azimuthal size of the simulation domain. For this reason, we have investigated the effect of the periodicity length imposed to the model studying different cases with L y = 0.67, 1.34, 2.68, and 5.36 cm (the latter corresponding to one quarter the real length of the SPT100 HT), using the physical parameters of Table 1 and numerical parameters fixed at the end of section Description of the model and numerical convergence study. The wavenumber spectrum of the oscillating field E y in the non-linear stage (t = 12 µs) is reported in Figure 8 for the two cases L y = 0.67 cm (red curve) and L y = 5.36 cm (blue curve). As evident, both reveals the presence of the same dominant cyclotron harmonic k 1 = 3,650 m −1 , but using a larger simulation domain, new low-k modes appear whose amplitude are comparable to that of the dominant cyclotron harmonic. We have verified that the large wavelength mode is not associated with the particle re-injection method by observing that it is still present in the case without re-injection (already at t = 2 µs before the electron energy grows beyond an acceptable value).
Finally, the instability is observed to stabilize when the amplitude is large enough to trap a significant portion of ions: the confirmation comes from the ion heating and rotation since the free energy dissipates by giving the ion distribution a long tail, rather than driving exponentially growing azimuthal electric FIGURE 8 | Amplitudes (normalized to the maximum) of E y spectral components k y using two different azimuthal domain size. Results correspond to input data reported in Table 1 and using the following numerical parameters N ppc = 500 and y field. Figure 9A shows the ion velocity distribution function (IVDF) projected along the azimuthal direction (red curve) and ( Figure 9B) the corresponding ion phase space at t = 12 µs. In Figure 9A are also reported for comparison the initial Maxwellian distribution at T i = 0.2 eV (blue curve) and the Maxwellian best fitting the computed IVDF (green curve). The latter gives a final temperature of T i = 10 eV, fifty times larger than the initial one. Moreover, the ion rotation occurs with the presence of different ion beams. In particular, two populations are distinguishable, one dominant rotating with a fluid velocity of about 13 km/s and another minor population rotating with 20 km/s. Figure 9B clearly shows the vortex formation typical characteristic of trapping in the azimuthal wave.

PARAMETRICAL STUDY AND SCALING LAWS
Since operating parameters ( Table 1) are very effective on the ExB EDI and as a consequence, on the global electron cross-field transport, we are interesting to conduct a study of the behavior of the electron cross-field mobility [calculated following Equation (8)] as a function of the most important input data used: neutral density, plasma density, magnetic field, accelerating electric field, and ion mass. In addition, the total cross-field mobility derived by electron momentum conserving equation taking into account azimuthal instabilities, Equation (3), is also used for comparison. The azimuthal correlation between the electron density and azimuthal electric field is computed in the PIC model through a double integral formula by using a time-averaged (over an interval T = 2 µs) value n e E y = Ly 0 T 0 dy Ly dt T n e (y, t)E y (y, t) (10) discretized using the Simpson's 1/3 rule [33]. Figure 10A shows the gas density dependence keeping the other physical parameters fixed to the values reported in Table 1. This study allows knowing the effect of collisions on the mobility since the contribution of instability does not depend on the electron-neutral collisional frequency. At very low gas density, the contribution to the electron cross-field mobility predominantly comes from the instability: its initial value is mostly determined by the second term on the right-hand side of Equation (3). However, we must remember that the particle re-injection technique adopted (see section Description of the model and numerical convergence study) corresponds to an artificial collision always present even in the case of very low gas density. As the gas density increases the collisional contribution becomes stronger (mixed regime) and the cross-field mobility increases reaching its maximum for n N = 8 × 10 20 m −3 . For much larger gas density the regime is collision-dominated [the first term on the right-hand side of Equation (3) is dominant]: the mobility change behavior decreasing as the gas density increases following the classical collisional dependence Equation (1) Figure 10B shows the electron cross-field mobility as a function of the magnetic field imposed. Both estimations follows the classical behavior Equation (1) µ ⊥ ∝ 1 B 2 (green curve), but the effective collisional frequency ν m,eff results 50 times larger than the real value reported in Table 1. This proof that in HT regime, the electron ExB drift instability does not induce full turbulence Bohm-like behavior µ ⊥, Bohm ∝ 1 B . Figure 10C shows the electron cross-field mobility as a function of the plasma density studying the pure effect of the instability on the electron transport (the gas density is kept fixed to n N = 2 × 10 18 m −3 ). In agreement with values found by previous PIC simulations [16,17], the mobility always increases with the plasma density starting from the pure collisional value assumed at very low plasma density (no instability) given by Equation (1) In Figure 10D the behavior of the electron mobility as a function of the applied axial electric field is reported. A larger E z0 gives larger ExB drift and therefore stronger electron cyclotron instability. As a consequence the mobility shows a continuous growth in the range analyzed. Figure 10E reports the effect of the ion mass, where the most relevant noble gases (He, Ne, Ar, Kr, and Xe) have been used. The different cases all use the same cross section in order to avoid possible contamination on the mobility due to the different cross sections, and measure in this way the pure ion mass effect. As evident, heavier ions lead to a strong reduction (larger than 10 times) of the electron cross-field mobility. The behavior is in agreement with the kinetic theoretical treatment of [11] that have shown as the onset conditions for instability are reduced and the unstable region is broader for lighter ions. Finally, the parametrical study has been used to understand the character of the instability after the non-linear evolution. The case study using physical parameters reported in Table 1 leads to very similar wavelengths, whatever is the character of the instability, λ ExB EDI = 1.76 mm and λ IAI = 1.62 mm. It is not useful to identify the nature of the instability at saturation occurred. However, by observing Equations (4) and (5), one can see that the ExB EDI wavelength scales with the axial electric field and radial magnetic field, while the ionacoustic wave scales with the electron density and temperature. In Figure 11A we have reported the k-spectrum of the electron density at t = 10 µs for four different electron density, n p0 = 2 × 10 16 , 4 × 10 16 , 5 × 10 16 , and 8 × 10 16 m −3 , while Figure 11B shows the dominant wavelength computed in the non-linear phase as a function of magnetic field B 0 (blue circles) and accelerating electric field E z0 (red circles) using as electron density n p0 = 1 × 10 17 m −3 . As evident, the dominant wavelength does not depend on the electron density (i.e., on the Debye length since the electron temperature computed does not change with the density) and perfectly scales according to the electron cyclotron drift instability Equation (4). However, it has to be pointed out that, the missing axial coordinate can play an important role in hindering the ion acoustic like transition since azimuthal-axial models [19,[26][27][28][29][30] have always shown a scaling according to ion acoustic like character of the instability Equation (5).

CONCLUSIONS
In this work, we have presented a numerical study, based on 1D Particle-in-Cell methodology, of the electron ExB drift instability in Hall thruster. It gives correlated fluctuations in electron density and electric field leading to enhanced electron cross-field mobility. The local representation is simplified since it neglects any effects of axial gradients, radial contributions, and assumes a constant axial electric field. Nevertheless, the description allows having a high statistical level, condition necessary to reduce the effect of numerical heating on the linear and non-linear evolution and on saturation of the instability. FIGURE 10 | Calculated steady state cross-field mobility as a function of (A) xenon gas density, (B) magnetic field, (C) plasma density, (D) axial electric field and (E) ion mass, using data reported in Table 1.
Results show that fluctuation spectra remain quantized (electron magnetization signature) predominantly retaining the fundamental cyclotron mode. In addition, using larger azimuthal simulation domains (forcing the periodicity length toward a realistic value), the development of long wavelength envelope comparable to the entire azimuthal extension, appears.
An accurate parametrical study has also been conducted to find the behavior of electron cross-field mobility as a function of the most important operative quantities: gas and plasma density, magnetic field, axial electric field, and ion mass. The electron cross-field transport shows a strong dependence on the plasma density and ion mass. The computed dominant wavelength scales FIGURE 11 | (A) Amplitudes (normalized to the maximum) of n e spectral components k y using four different electron density. (B) Dominant wave detected in the non-linear phase as a function of magnetic field B 0 (blue circles using E z0 = 200 Vcm-1) and accelerating electric field E z0 (red circles using B 0 = 200 G). Results correspond to input data reported in Table 1 and using the following numerical parameters N ppc = 500 and y λ D = 0.36.
according to the ExB EDI and no transition toward ion acoustic like instability is observed.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
FT developed the numerical method and led the manuscript preparation. ZA contributed to the code development and to the article preparation. FT, ZA, and MS contributed to the analysis and discussion of the results.

FUNDING
This work was supported by the Apulia Space Project (grant PON03PE_00067_6).