Magnetic nozzle performance in a cluster of helicon plasma thrusters

A numerical study of the plasma dynamics in a Helicon Plasma Thrusters’ (HPT) cluster is presented. For the first time in the literature, the three-dimensional (3D) plasma dynamics occurring in the plume of a HPTs’ cluster is analyzed. The physical investigation relies on ProPic, a 3D particle-in-cell (PIC) code specifically designed to simulate the plasma dynamics in magnetic nozzles and in a non-axi-symmetric domain. The code has been validated against experiments reported in the literature and cross-validated with Starfish, an open-source two-dimensional PIC software. The physical investigation has revealed an interesting mutual influence between the thrusters that constitute the cluster. Three significant phenomena that affect the cluster’s performance have been identified. The first phenomenon is related to the effect that clustering has on the shape of the magnetic field lines and, in turn, on the divergence angle of the plume. The second phenomenon is related to electron currents flowing among different thrusters, which affect the potential drop across the plume. The third phenomenon is related to the effect that neighboring thrusters have on the plasma potential map and, in turn, on the expansion of the ions.


Introduction
Electric propulsion is the subject of intense research activity despite ion and Hall effect (HET) thrusters being mature technologies with a consolidated flight heritage. Indeed, these systems are too complex and high cost for several applications * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
in the field of small-satellites (i.e. spacecraft's mass below 200 kg) [1][2][3][4][5][6][7]. In recent years, there has been an increased focus on the development of Helicon Plasma Thrusters (HPT). This technology consists of a radio frequency (RF) plasma source, operating in the MHz range, coupled to a magnetic nozzle (MN). Plasma is produced within the source by coupling RF power to the electrons using Joule and/or Helicon wave heating processes [2,3,[8][9][10][11][12][13]. When applying the MN, electric fields arise spontaneously, which boost the propulsive performance converting the electrons' thermal energy into ions' kinetic energy [14][15][16][17][18]. Moreover, the energetic electrons that surpass the potential drop occurring within the plume can neutralize the accelerated ions [19][20][21]; thus, HPTs do not need a neutralizer [22,23]. In HPTs, no electrodes are exposed to the plasma that might be damaged by ion sputtering or thermal loads [24]. Therefore, HPTs are expected to be mildly affected by erosion issues while ensuring reduced complexity and costs. This makes HPTs a good option in the growing market of low-budget satellites [25]. The main drawback of HPTs is the relatively poor conversion efficiency from RF to thrust power, which was about 0.5%-5% in early experiments [26][27][28] and approached 30% in recent academic experiments [22,29]. Despite the remarkable progress in HPT technology, the efficiency of commercially available HPTs is still significantly lower compared to Ion and Hall effect thrusters that operate at similar power levels [24,[30][31][32][33].
A cluster of thrusters may offer several advantages over a single similarly powered device, despite the risk of lowering the overall thruster efficiency and increasing the dry mass [34,35]. For example, clustering enables to improve the system's reliability because of redundancy. Moreover, the propulsive performance can be adjusted by turning off one or more thrusters. This throttling allows the cluster to operate at lower power without running any individual thrusters at off-design conditions [36,37]. Clustering can be advantageous for missions where power or propulsive needs vary over time, or it is necessary to implement thrust vectoring [38]. Additionally, clustering offers a high degree of scalability. Indeed, once the technical issues of operating a cluster are understood, a single flight-qualified engine can support a wide range of missions requiring various power levels by employing the appropriate number of thrusters [39]. Since 2002, the Air Force Research Laboratory and the University of Michigan have conducted numerous experiments on clustering of HETs [39][40][41][42]. Measurements have shown that clustering can cause significant changes in the ion energy profiles [41], the ion current density differs from a simple linear superposition of individual contributions [40], and the magnetic field of each thruster is affected by the neighboring devices [40]. Despite these findings, there are still open questions, such as how the cluster affects the neutralization process when a cathode is present [35], and how the plasma plumes interact with each other and with the spacecraft [34]. Nonetheless, the total thrust generated by a cluster of HETs is observed to be only 3.6% higher than the sum of the individual contributions [42], suggesting a relatively mild interference between multiple devices.
Clustering HPTs is highly appealing, provided the reliability and low cost of this class of thrusters. Nonetheless, only a few studies are focused on this topic, provided the complexity of evaluating the propulsive performance both experimentally [43] and numerically [44]. Experiments performed at the University of Washington in 2017 suggest that clustering HPTs might result in an increment of the performance provided the increase of both plasma density and ion exit velocity [43,45]. To date, only two-dimensional (2D) preliminary analyses of magnetic configurations representative of an HPTs' cluster have been accomplished [44]. Nonetheless, asymmetric magnetic field configurations arising in an HPTs' cluster are expected to affect the plasma dynamics in the MN and, in turn, the overall propulsive performance [43][44][45]. The lack of accurate three-dimensional (3D) simulations is mainly due to their computational burden, which has become affordable only recently [46,47].
The MN's dynamics can be modeled following different approaches, ranging from multi-fluid models to fully kinetic ones [48,49]. 2D fluid models are powerful tools for understanding the main phenomena governing the MN [20], but their closure remains an elusive problem. One-dimensional (1D) stationary kinetic models have allowed the analysis of the ion and electron heat fluxes, and the response to non-Maxwellian features of the ion and electron velocity distribution functions (VDFs) [50,51]. However, except for 1D cases, solving the Boltzmann equation directly is often computationally intensive [49]. Moreover, fluid and kinetic continuum approaches must make assumptions regarding VDF, one of the main impact parameters in magnetized plasma expansions [51]. On the other hand, the fully kinetic particlein-cell with Monte-Carlo Collisions (PIC-MCC) method is the numerical strategy with the lowest level of assumptions. The trajectories of ensembles of particles, called macro-particles, are integrated, accounting for the self-consistently computed electric and magnetic fields and collisions [52]. Nonetheless, the computational burden associated with this approach might be intense [46]. 2D-axisymmetric PIC-MCC solvers are commonly employed to simulate magnetized plasma plumes [46,[53][54][55][56][57]. However, this approach suffers from some limitations: macro-particles might be unphysically accelerated near the symmetry axis [58], azimuthal oscillations and anomalous transport shall be assumed [46,59,60], and asymmetric magnetic field configurations or spacecraft geometries cannot be handled [44,61]. The number of 3D PIC-MCC simulations and benchmarks remains limited [46] due to the intensive computing resources required. Nonetheless, the adoption of efficient parallelization strategies [62] and numerical acceleration schemes [53] is allowing 3D PIC-MCC to be performed in a reasonable computational time (e.g. less than 24 h per simulation [47]). This is the first time a 3D numerical study of the MN dynamics in a HPTs' cluster is presented. The physical investigation relies on ProPic, a 3D PIC-MCC code specifically designed to simulate the plasma dynamics in magnetic nozzles [47]. Section 2 summarizes the numerical method used to simulate the MN. In section 3, the numerical approach has been validated against measurements of plasma density and thrust [63,64], results have also been verified against the 2D-axisymmetric PIC code Starfish [53]. In section 4, the validated approach has been exploited to investigate the plasma expansion in a HPTs' cluster. Conclusions are given in section 5.

Methodology
ProPic simulates the behavior of a magnetized plasma plume in a 3D space (x, y, z). An overview of the simulation domain  (I, II, III), and Robin to (IV). Ions and neutrals are absorbed on (I, II, VI), on (III) ions are absorbed and neutrals reflected. Electrons are absorbed on (I, II, III) and selectively reflected on (IV) based on an energy criterion. and the boundary surfaces is shown in figure 1. The domain is delimited by the external boundary (IV) consisting of a cylinder of height H B and diameter D B . The spacecraft (III) is a parallelepiped of dimensions H s , L s , W s . The plasma source has been omitted because the aim of this work is solely to simulate the expansion of the plume. Ions, electrons, and neutrals are injected in the computational domain through boundaries corresponding to the thruster outlet surfaces (I and II) of diameter d o and oriented perpendicularly to the z-axis. This work describes the modeling of a cluster made up of two thrusters, but the simulation strategy can be easily generalized for any number of thrusters.

PIC simulation
In the PIC method, the plasma is treated as an ensemble of macro-particles p = 1, . . . , N p with positions x n p = ⟨x n p , y n p , z n p ⟩ and velocities v n p = ⟨v n px , v n py , v n pz ⟩ [65][66][67]. The mass and charge of each species are referred to as m p and q p , respectively. The typical PIC cycle is illustrated in figure 2: at every time iteration, the charge density is distributed to the mesh nodes according to the particle positions using a linear interpolation scheme (step 1) [65]. Maxwell's equations have not been solved in full since usually less than 1% of the total RF power coupled to the plasma is deposited outside the source [68]. Furthermore, the currents (mainly azimuthal [2]) induced in the MN are low enough that the self-generated magnetic field can be disregarded (errors of less than 1% are expected [2]). Therefore, the magnetic field B is static and matches the one imposed by the MN [9]. As a result, the Maxwell's equations are simplified to the Poisson's equation: with ϕ the electrostatic potential and E = −∇ϕ. The charge density computed at step 1 is used to solve for the selfconsistent plasma potential ϕ according to Poisson's equation via a finite element method with a preconditioned conjugate algorithm [69] (step 2). Then, particle motion is solved explicitly using the standard leap-frog Boris algorithm (step 3) [70,71]: Collisions (step 4) have been implemented according to the MCC method described in [72]. Various collision types have been implemented, including electron-neutral elastic scattering [73], electron-neutral excitation [73], electron-neutral ionization [73], ion-neutral elastic scattering [73], ion-neutral charge exchange [74], neutral-neutral elastic scattering [73], electron-ion recombination [75]. Data on charge exchange cross-sections is obtained from [74,76] for xenon propellant, and from [73] for argon gas. This sequence of operations is repeated at each time iteration, self-consistently evolving the particles and the electric field states until steady state is reached.
The steady state of a simulation is assumed to be reached when the number of particles inside the domain is constant. This condition is met after a certain number of time steps (α) to avoid a false-positive result due to PIC noise. This is implemented by requiring [77]: where ξ is a user defined number that represents a steady state condition, usually around 10 −4 , h is the current time step and N t p is the number of macro-particles in the system at the time step t. As a general guideline, α is assumed to be around 50. The thrust T produced by the engine is calculated by summing the momentum exchange and pressure contributions of each species k on the external boundary surface IV [78]: with k the external normal to IV andẑ the unity vector parallel to the z axis.
To speed up the simulation, the free space vacuum permittivity ε 0 has been increased by a factor γ 2 : Because of equation (6), sheaths are thicker and plasma dynamics are slower, allowing to use of a coarser grid and a longer time-step [79,80]. A sensitivity analysis in the γ scaling can be found in [53].

Boundary condition on the Poisson's equation
The boundary condition that has been imposed for the closure of equation (1) are: • Dirichlet condition on the thruster outlet surfaces (I, II) and spacecraft surface (III), with ϕ t I , ϕ t II and ϕ t III the corresponding potential at the time step t. ProPic can handle potentials that vary in space on the thrusters' outlet or spacecraft surfaces, such as those caused by spacecraft charging in more complex designs [81,82]. But, for simplicity, it is assumed that there are no potential gradients in I, II, and III. The values of ϕ t s change over time, and their steady state value cannot be predicted beforehand and must be calculated as part of the solution, as explained in section 2.4. • Robin condition on the external boundary (IV) [83]: with r the vector going from the geometric center of the simulation domain to IV. Equation (7) derives from the assumption of a monopole decay 1/r of the potential and that at the infinity ϕ ∞ = 0 V [83,84]. The derivation of equation (7) is provided in [53].

Boundary conditions on particle kinetics
From now on, pedix i, e, and g will refer to ions, electrons, and neutrals. At each time step, ions, electrons, and neutrals are injected from the thruster outlet surfaces (I, II). For all species, a Maxwellian VDF is assumed [53,85,86] with reference temperatures T k,s for k = i, e, g and s = I, II. A drift velocity equal to v t = Mv b in the direction perpendicular to the thruster outlet surfaces is imposed for ions, where M is the magnetic Mach number and v b = k B T s,e /m i is the Bohm speed. In typical HPT configurations, it is usually assumed that the Mach number is close to one [11,20,53,73]. A diffusion drift velocity v g = k B T s,g /2π m g is imposed on the neutrals [8].
Generating a stable, steady-state plume requires an energybased approach [47,53,54,87]. Assuming a steady-state electric field, the total energy of each electron can be defined as where ϕ IV is the potential at the external boundary and E tot is a constant conserved quantity of the motion in the collisionless case. Two populations of electrons can be identified based on their total energy. Electrons with E tot < 0 do not have enough energy to cross the potential drop across the plume, so they are forced to return to the plasma source at a certain distance downstream. Electrons with E tot ⩾ 0 have enough energy to cross the potential drop and can escape to infinity. Thus, electrons that reach (IV) with a negative E tot are reflected; instead, electrons with a positive E tot are absorbed. The main assumptions in this energy-based boundary condition are that the plasma is collisionless beyond the domain and that the magnetization is not strong enough to cause reflection of highly energetic electrons [50]. Since ions are accelerated outward by the ambipolar electric field [88], they are absorbed in (IV). Also, neutrals are absorbed in (IV) since no special treatment is requested for non-charged particles. Ions and electrons that reach the surfaces of the spacecraft, including the thruster outlet surface (I, II, III), are absorbed. Neutrals are absorbed at (I, II) and reflected at (III).

Control loop
A control loop has been implemented to compute the values of ϕ t I , ϕ t II and ϕ t III in steady state, and to maintain a quasi-neutral and current-free plume. This approach is similar to other codes that model the interaction between a spacecraft and a plasma plume, where an electric circuit connects the two [48,83]. The net current I s,k of the species k at the surface s is defined as I s,k = I inj s,k − I abs s,k , where I inj s,k is the current that is injected into the computational domain from the surface s and I abs s,k the current that is absorbed and leave the computational domain at the same surface. For each surface s = I, II, III, the potential ϕ t s is updated according to [11,47,53]: with I t−1 s,e and I t−1 s,i respectively electron and ion net current from the previous time step and c s a positive control coefficient. It is worth noting that the value of c s must be carefully chosen to find a balance between the fast convergence of ϕ s and the stability of Poisson's solver. A dedicated sensitivity analysis can be found in [53]. This method ensures that the system evolves self-consistently and guarantees that, once it reaches steady state, the ion and electron net currents are equal (I s,e + I s,i = 0) for s = I, II, III.
In general, the density n s,k of the species k at the surface s = I, II can be computed as [47]: with Γ s,k the particles flux density through s and µ k the mean of v k,z at the same surface. In typical HPT configurations, the plasma is mesothermal (i.e. T e ≫ T i ), ions are accelerated outward by the ambipolar electric field and |µ i | ≈ v t [47,89]; thus, the injected ion flux can be considered constant. Specifically, I inj s,i reads: withṁ s the propellant mass flow rate. Instead, µ e is unknown a priori, thus the quasi-neutrality condition n s,i ≈ n s,e is obtained at the steady-state by adjusting the injected electron flux. Specifically, I inj s,e is updated each time step according to with h s a positive control coefficient. Equations (9) and (12) guarantee that quasi-neutrality and current-free conditions are respected at the steady-state. It should be noted that equations (9) and (12) are general and can be applied to any number of surfaces 's', regardless of whether they are considered the outlet of a thruster or the surface of the spacecraft. Thus, the number of thrusters, whether one, two, or more, does not affect the applicability of the proposed methodology. It is worth noting that some particles might flow among thrusters in a cluster. Nonetheless, ions are accelerated outward by the ambipolar electric field. Thus, it is expected that this phenomenon is negligible for ions. On the contrary, a non-negligible flow of electrons might be exchanged between neighboring thrusters [44]. In the following, the electron current that is injected from I and is absorbed in II is noted with I * I,e , the electron current that is injected from II and is absorbed in I is noted with I * II,e . As proven in appendix A, the latter phenomenon does not prevent the possibility to establish a current-free condition on each surface given the possibility to control I inj s,e (see equation (12)). It is important to choose carefully the initial values of the control loop to guarantee convergence to steady state solutions. ϕ 0 III is set to 0 V since, for a typical small-satellite mission, the plasma plume and the ambient plasma are not expected to charge the spacecraft surfaces excessively (i.e. in the order of few Volts) [30,90]. The initial values ϕ 0 I and ϕ 0 II are set according to a theoretical value calculated by assuming current-free conditions, no magnetic field, and that the electron energy is conserved [86]. ϕ 0 s are obtained by solving It is worth noting that the control loop described in this section is derived from previous studies [11,47,48,53,54,86,87,91,92] in which successful validations are presented.

Experimental validation
The results of ProPic are compared against the measurements presented in [63,64]. Numerical results are compared against the plasma density and temperature measured with a 3 mm-diameter planer Langmuir probe, whose uncertainty is about ±15%. The thrust, directly measured using a pendulum thrust balance, is also used as a benchmark; typical errors are within ±15%. The results of ProPic have also been compared against those obtained using the 2D-axisymmetric PIC code Starfish, which has been used previously to model Hall thruster channels [93], as well as the plumes of ion [94], magneticallyenhanced vacuum arc [95], and Helicon [11,47,53] thrusters. The methodology proposed in section 2 is similar to the one used in Starfish. Further details can be found in [53].  respectively. Results on the y-z plane are not presented due to the axisymmetry of the domain. The results provided by the two codes show similar behavior (difference smaller than 20%). Specifically, the Starfish's outputs are noisier near the z-axis than the ProPic ones. This is because the 3D code has fewer assumptions on the z-axis which is not treated as a boundary of the computational domain. This, in general, leads to more accurate results [46]. The axial plasma density is shown in figure 3(c). The difference between ProPic and the experimental data is within the reported uncertainty bands. Starfish's results are outside the uncertainty bands only for x < 5 cm, with n overestimated by approximately 20%. This bias could be attributed to the reflective boundary condition imposed along the z-axis in the 2D-axisymmetric code [53].
In MN expansions, electron cooling can be characterized by the polytropic relationship: T e /T e,I = (n e /n e,I ) γe−1 , where γ e denotes the polytropic index [96]. Analytical estimations of γ e are available in literature and, according to [88], its value is 1.27 for argon propellant. A linear regression between log 10 T e and log 10 n e can be used to extrapolate the value of γ e from both experimental values and simulation data. The experimental results yield an approximate value of γ e ≈ 1.18, while the ProPic results indicate γ e ≈ 1.27, in perfect agreement with [88]. The discrepancy between the value extrapolated experimentally and the numerical simulation, which is roughly 7%, falls within the uncertainty error range. Table 2 compares the measured thrust against predictions obtained with the PIC simulations. Both ProPic and Starfish predict the thrust with an error lower than 10% (for ProPic the disagreement is about 1%). A lower accuracy of the 2D-axisymmetric code was expected provided the unphysical macro-particle accelerations that might occur near the symmetry axis [58]. Notably, the plasma potential and ion speed fields were not measured, but the density and thrust experimental data showed good agreement with PIC results. Therefore, based on the momentum equation, it is reasonable to assume that the velocity field and the potential drop across the plume are reasonably estimated.
The main disadvantage of ProPic with respect to Starfish is the greater computational cost. ProPic required a machine with at least 32 GB of RAM and took approximately 24 h to run the simulation using 16 processors. Starfish required less than 8 GB of RAM and took approximately 8 h to run the simulation using 8 processors. As shown in table 1, a larger γ scaling has been used in ProPic to get results in a reasonable amount of time. Despite the larger scaling, the number of mesh elements  N m in the 3D simulation is more than 20 times the one in the 2D case. Also, the number of macro-particle N p at steady state is 2 times larger in ProPic than in Starfish.

Supplementary validation
A second validation has been added to ensure the reliability of ProPic in estimating the plasma potential and ions' speed. The numerical code is validated against the experiments conducted in a Piglet Helicon reactor [97]. This mimics the approach presented in [53] to validate the implementation in Starfish of the numerical methodology discussed in section 2. The Piglet reactor is filled with argon gas at a pressure of 0.04 Pa. The magnetic configuration is generated by an electromagnet producing a magnetic field intensity of B 0 = 4 G at the throat. A retarding field energy analyzer (RFEA) has been utilized to measure the local plasma potential and the ion energy distribution functions. RFEA measurements have a given uncertainty of ±5% [53]. The input parameters for ProPic are given in table 3. The axial plasma potential is shown in figure 4: the PIC model reproduces the experimental trend with accuracy within the stated uncertainty. The ions' speed obtained with ProPic on the z-axis at z = 100 mm equals 8300 m s −1 . This value perfectly agrees with the one obtained by the RFEA measurements (the difference is less than 2%).

Results
The numerical investigation of the MN dynamics and performance in a HPTs' cluster has been conducted accounting for two thrusters, as shown in figure 1.  [11]. The magnetic configuration of each HPT is created by a solenoid yielding a magnetic field strength of 920 G at the center of the corresponding outlet surface. Asymmetric magnetic configurations that form in a HPTs' cluster are expected to influence the plume expansion [43][44][45]. Thus, the purpose of this study is to examine, for the first time, how clustering impacts the plasma dynamics within the MN. To simplify the investigation, a uniform plasma density profile has been assumed on the thrusters' outlet despite the non-uniformity observed experimentally [2]. However, in appendix B, we have conducted a sensitivity analysis to assess the impact of neglecting plasma non-uniformity at the thruster outlet. Similarly, collisions have been ignored. Nonetheless, to provide a more comprehensive understanding of the cluster's dynamics, we have included appendix C where results obtained with and without collisions are compared.
From now on, the notation I and II will be used to refer to the thrusters associated with surface I and II, respectively. The same notation will be used for the corresponding solenoids. The numerical investigation is conducted by comparing four different configurations: • Configuration A, which is simply a single thruster being on, is the benchmark against which the performance will be compared. The study of configuration B might be representative of a case in which the magnetic field is generated with permanent magnets and not electromagnets [30]. Moreover, the condition of having one thruster on and the other off could be used, among other things, to meet the requirements of a particular attitude control maneuver or to limit power consumption. Configuration C is the condition when two identical HPTs operate simultaneously. The magnetic field generated by configurations A-C contributes to having a non-zero magnetic moment of the satellite, causing potential problems with attitude control [98]. This problem can be solved by configurations similar to the D one, which has a resulting magnetic moment close to zero. Configuration D is also representative of the MN of the new arch-shaped thruster concepts (MAT) [44]. Indeed, MAT thruster features a 'C'-shaped discharge chamber and a toroidal magnetic field made up of two magnets with opposite polarities. The magnetic field corresponding to the four configurations is shown in figure 5.

Plasma density
The ion plasma density n i on the x-z plane is shown in figure 6, with white contours representing the magnetic field lines. In configuration A, the plume follows the magnetic lines of the axisymmetric MN until detachment, which happens approximately at z = 0.35 m. This behavior is similar to what was observed in previous studies of the plasma expansion in the MN of a single thruster [47,53,54,99]. In configurations B and C, the magnetic field lines are no more parallel to the zaxis near the thruster outlet surface (I and II respectively). The plasma, following the magnetic lines, is deviated radially; a more quantitative discussion on the divergence angle will be provided in the following section 4.4. In configuration D the lines of the MN are connected and the two plasma plumes that exit from I and II, respectively, merge into one flux. It is interesting to observe that the results of configurations C and D, in which two thrusters are on, are not a linear combination of the results of configuration A, in which only one thruster is operated.

Charge density and plasma potential
The plasma potential ϕ on the x-z plane is shown in figure 7. The results of configuration A are similar to what was seen in previous studies, with a potential drop forming between the thruster outlet and the surrounding 'infinity' to maintain a current-free plume [2,47,53,54,99]. In configuration B, the non-axisymmetric field lines alter the potential shape. In configurations C and D, the potential field is significantly changed due to a second positive peak in correspondence with the thruster II. The presence of two potential peaks has a strong influence on the cluster's dynamics since the potential drop across the plume is equal to the steady state values of ϕ I and ϕ II obtained from the control loop (see equation (9)).
It is worth noting that these findings are not perfectly consistent with previous results [100][101][102] that highlighted the presence of a plasma potential well in correspondence with the plume periphery. This potential well structure can be interpreted as the consequence of charge separation that results from ions with sufficient transverse energy overshooting the attached electron population [103]. As demonstrated in [53], the phenomenon of charge separation leading to the potential well structure is primarily associated with collisions. Moreover, [104] has shown that this phenomenon decreases significantly as the ratio of ion temperature to electron temperature (T i /T e ) decreases. Figure 8 shows the relative charge density (n i − n e )/n e , with white contours representing the magnetic field lines. Results are similar of what seen in [53]: within the plume, plasma is substantially quasi-neutral and there is a charge separation across the plume periphery. The existence of a strongly negative space charge along the outermost magnetic layer and a prominent positive charge beyond this layer can create a radial potential well. However, due to the collisionless nature of our approach and the low T i /T e ratio in our simulation, the charge separation is not significant enough to generate a potential well. It is unlikely that this will have a significant impact on the physical analysis of the system, since the intensity of the potential well remains relatively low, on the order of a few volts [100].

Control loop
The steady state values of ϕ I and ϕ II depend on the currents flowing through each thruster. As demonstrated in appendix A, the control loop implemented ensures the current-free condition on each surface of the system but the contributions that determine I s,k shall be analysed to grasp the cluster's dynamics. Specifically, electrons that exit from thruster I are attracted by the positive potential peak formed in correspondence of thruster II, and vice-versa. Instead, this phenomenon is negligible for ions that are accelerated outward by the ambipolar electric field. If I * II,e increases, then I abs I,e will also increase, which results in a drop of the electron net current I I,e . Since the ion net current I I,i is almost constant, a lower potential drop ϕ I is needed to keep the plume current-free; the same happens for II. Provided that electrons are strongly magnetized in usual HPTs [47], the electron current that flows between the two thrusters depends heavily on the magnetic topology. Thus, two different behaviors occur depending on whether the solenoids have the same polarity (configuration C) or opposite polarities (configuration D). In configuration C, the electrons that exit the thruster outlet with a sufficiently large amount of kinetic energy can arrive far enough to detach from the magnetic field lines [105]. Once detachment occurs, they could be attracted by the positive potential of the other thruster and eventually be trapped by the magnetic field lines that lead to the second thruster outlet. It is expected that this current is relatively low, because only the most energetic electrons can detach from the magnetic field lines and, once detached, they could have enough energy to escape toward infinity. In configuration D, because the magnetic lines are closed, the electrons follow the magnetic field that connects the two thrusters and can get close to the second outlet, eventually being absorbed, without detachment. Therefore, it is expected that a significant number of electrons that exit from one thruster's outlet, even those with relatively low energy, will flow to the other one, resulting in a large current flow. Table 5 shows the steady state values of ϕ I , the absolute value of the electron currents that flow from II to I |I * II,e | and the injected ion current I inj I,i . Because of geometric symmetry, ϕ II , |I * I,e | and I inj II,i for configurations C and D are equal to ϕ I , |I * II,e | and I inj I,i , respectively. The ion's current flowing between the two thrusters is not reported because it is negligible. Results confirm what is expected: the electron current that flows between the two thrusters in configuration D is 157 times the one in configuration C. The potential drop ϕ II in steady state in configuration D is 17.4 V lower than in configuration C. It is worth noting that the electron current flowing between the two thrusters in configuration C, which is equal to 0.05 A, is not negligible when compared to the ion injected current whose value is 0.145 A. This determines a steady state  value of the potential drop across the plume that is 3 V lower than in configurations A and B. Figure 9 shows the normalized VDF of the electron population in the proximity of the surface I for configurations C and D. Electrons are sorted depending on the thruster from which they have been emitted. Their distribution function is denoted VDF e,I and VDF e,II , respectively. The notation VDF e indicates the distribution function of the entire electron population. Table 6 displays the average temperature T e,s and density n e,s of the electron population in the proximity of the surface I. The subscript s refers to the surface from which electrons are ejected. Results confirm that the magnetic topology has a major influence on the electron's dynamics. In configuration C, T e,II is almost twice T e,I and n e,II is two orders of magnitude lower than n e,I because only the most energetic electrons that exit from II can detach from the magnetic field lines and arrive in the proximity of I. It is interesting to note that n e,II is so small that VDF e and VDF e,I are almost identical. Table 7 shows the probabilities associated with VDF e,II . It is interesting to observe that VDF e,II is truncated for v > 4.6 × 10 6 m s −1 for configuration C. In fact, the electron population ejected from thruster II is composed completely of 'trapped' electrons with E tot < 0 since electrons with E tot > 0 have enough energy to reach infinity once detached. In configuration D, T e,II is only 1.9 eV higher than T e,I and n e,II is of the same order of magnitude of n e,I because electrons emitted from thruster II do not need to detach to get close to I.    Table 7. Probabilities associated with the VDF e,II of figure 9.
P(V > 4.0 · 10 6 m s −1 ) P(V > 4.6 · 10 6 m s −1 ) Therefore, a significant number of electrons flow from II to I, even those with relatively low energy. Interestingly, the number of electrons arriving from II is large enough to cause VDF e to be different from VDF e,I . It is also interesting to observe that VDF e,II is not truncated for configuration D.

Plasma acceleration and ions' flux
Since ions are accelerated outward by the ambipolar electric field, their speed is strongly related to the potential drop across the plume. The relation between the ion speed |v i | and plasma potential ϕ is shown in figure 10, that reports ϕ and |v i | along the line 0 < z < 0.8 m, y = 0 m, x = −0.05 m (i.e. axis of the thruster I). In configuration C, the reduction of 3 V in the potential drop causes a reduction of 2.1% of the ion speed at z = 0.8, y = 0 m, x = −0.05. In configuration D, a reduction of approximately 20 V in the potential drop causes the speed to decrease by 21.0%. The speed of the ion population injected from I |v i,I | on the x-z plane is shown in figure 11. Figure 12 shows the current density J i,I of the ion population injected from I on the z ref = 0.1 m plane. The red contour in figure 12 corresponds to the minimum area whose integral of J i,I is equal to 0.85 · I inj I,i = 0.123 A. The maximum distance between the red contour and the I center is denoted with r θ . The divergence angleθ I of the four configurations, reported in table 8, is computed according to:θ In configurations B and C, the plasma following the magnetic lines is deviated radially, generating a strongly non axysimmetric current density J i,I . Additionally, when thruster II is turned on, it acts as a positive potential barrier, causing an even more pronounced asymmetry in the ion flux. This is especially noticeable when comparing ion speed and the shape of the red contour in configurations B and C. It is interesting to observe that the divergence angle is similar in configurations A-C. Instead, in configuration D, the plasma follows the closed magnetic field topology, causing a reduction of the divergence angle of ≈15.0 • . The significant distance between the two thrusters ensures that the positive potential peak of thruster II affects only a small fraction of ions injected from I. This is because the current density at the coordinates (x = −0.05, y = 0, z = 0.1 m) is over ten times higher than the current density at the point (x = 0.05, y = 0, z = 0.1 m).

Discussion
The results presented in this section suggest that three different phenomena can influence the performance of a HPTs' cluster. The first phenomenon is related to how the HPTs' clustering affects the shape of the magnetic field lines. The plasma follows the magnetic lines; therefore, the magnetic topology affects the shape of the plume and, consequently, the direction of ion flux. The results show that the divergence angle decreases when the magnetic field has a closed line configuration, such as when thrusters with opposite polarities of the magnetic field are arranged as in configuration D. Notably, the decrease in the divergence angle in a closed line configuration has also been observed in previous simulations of a MAT thruster [44]. A-C exhibit nearly identical divergence angles, but their plume directions differ primarily due to the orientation of the magnetic field lines. It is worth to observe that in previous experiments discussed in [43,45], solenoids with the same polarity were tilted to have a magnetic field more parallel to the z axis in the thruster outlet proximity. This design has enabled better collimation of the ion flux and suggested an improvement in cluster performance compared to a single thruster.
The second phenomenon is the effect that a second thruster has on the potential drop across the plume since it affects the acceleration of the ions. This effect appears strongly connected  to the electrons that flow between thrusters, affecting the electron current and, consequently, the potential drop needed to maintain the plume current-free. The results show that the electron current flowing between thrusters is strongly dependent on the topology of the magnetic field lines and is significantly enhanced for a closed line configuration. The third phenomenon is the effect that a second thruster has on the potential shape. When the second thruster is turned on, it creates a second potential peak which acts as a barrier for ions that have exited from the first thruster, resulting in a non-symmetrical ion flux. The influence of this potential peak is expected to depend strongly on the distance between the thrusters. In our study, the thrusters seem sufficiently far enough that this phenomenon does not significantly affect the system's dynamics. Table 9 shows the thrust T and the specific impulse I sp obtained for the four configurations. The plasma profiles in configurations C and D are not the linear superposition of configurations A or B. Thus, the thrust produced by the cluster is not simply double the one obtained with a single thruster. To quantify the cluster's influence on the propulsive performance, the specific impulse obtained in configuration A is compared with the ones obtained in configurations B, C, and D. Being in configuration B, the thruster II turned off, the clustering can affect the thruster performance only by modifying the shape of the magnetic field lines (the first phenomenon previously identified). The potential drop in configuration B is approximately the same as in configuration A. The reduction of I sp by approximately 5% compared to configuration A is therefore caused by the deflection of the ion flux direction due to the first phenomena.
Configuration C has the same magnetic field as configuration B, but the potential drop in configuration C is approximately 3 V lower than in configuration A due to the second phenomenon. The thrusters in configuration C are far enough that the third phenomenon does not affect the divergence angle and, therefore, the performance. Consequently, configuration C has an I sp that is reduced by approximately 8.6% compared to configuration A due to the deflection of ion flux caused by the first phenomenon and a decrease in ion acceleration caused by the second phenomenon.
In configuration D, the first phenomenon causes the two plasma fluxes that exit from I and II to combine into one flux with a reduction in the divergence angle. The potential drop in configuration D is approximately 20 V lower than in configuration A due to the second phenomenon. The acceleration of the ion is much lower in configuration D than in other configurations due to the significant reduction of the potential drop. Consequently, despite the reduction in the divergence angle, the specific impulse I sp in configuration D is 23.3% lower than in configuration A.

Conclusions
In this work, a 3D numerical study of the plasma dynamics in a HPTs' cluster has been presented for the first time. The physical investigation was conducted using ProPic, a 3D PIC-MCC code designed to simulate the plasma dynamics in magnetic nozzles and in a non-axisymmetric domain. The code has been validated against experiments reported in the literature and cross-validated with Starfish, an open-source 2D PIC software. The physical investigation has revealed a non-trivial mutual influence of the thrusters in a HPTs' cluster. Three different phenomena that affect the cluster's performance have been identified.
(i) The first phenomenon is related to how the HPTs' clustering affects the shape of the magnetic field lines. Our results show, as previously observed in simulations of a MAT thruster [44], that the divergence angle decreases when the magnetic field has a closed line configuration. When the solenoids of the two thrusters have the same polarity, the direction of the ion flux deviates resulting in a slight decrease of the propulsive performance. It is worth noting that this phenomenon could increase the collimation of the ion flux and, consequently, the thruster performance with an appropriate design, such as inclining the direction of the solenoids as in the previous experiments presented in [43,45]. (ii) The second phenomenon is the effect that a second thruster has on the potential drop across the plume. This effect is linked to the electrons that flow between thrusters, reducing the net electron current and, consequently, the potential drop needed to keep the plume current-free. Our results show that the electron current flowing between thrusters is strongly dependent on the topology of the magnetic field lines and is enhanced when the solenoids of a cluster have a closed line configuration, resulting in a considerable reduction of the performance. (iii) The third phenomenon is the effect that a second thruster has on the potential shape. When the second thruster is turned on, it creates a second potential peak which acts as a potential barrier for ions that exited from the first thruster. In this study, the thrusters seem to be sufficiently far that the divergence angle is not significantly affected by this phenomenon. It is expected that the influence of the potential peak strongly depends on the distance between the thrusters and on their number. Additional investigation is required to determine if this phenomenon can be utilized to enhance the performance of clusters or change the direction of the ions' flux, or if it is generally negligible in practical designs.
Future work will focus on improving the performance of a HPTs' cluster. This will involve looking at different inclinations of the solenoids and changing the number of thrusters and their positions. Further investigation is needed into the physical role of collisions and the interactions with the plasma sources. The particle flux that flows between the thrusters is expected to affect the source performance, and this could be analyzed using a fluid code such as 3D-VIRTUS [11].

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files). and net current I s are reported for ions and electrons. For every surface the sum of the ion net current I i and the electron current I e is zero. Therefore the control loop described in section 2.4 ensures the current-free condition. In configuration D, there is a significant difference in the electron-injected current (I inj e ) compared to configurations A-C. This is attributed to the substantial electron current flowing between the two thrusters as indicated in table 5.

Appendix B
A sensitivity analysis has been conducted to examine the nonuniform radial plasma density profile at the thruster outlet and its impact on the propulsive performance. Configuration C has been studied considering that the plasma at the source outlet exhibits a self-similar normalized radial profile that reads [106]: Here, R S is the source tube radius, which has been considered equal to d o /2, ρ = r/R s represents the normalized radial coordinate, with s = 2 and t = 6. The parameter h r can be computed according to [106]: where r ci is the ion cyclotron radius and c = 0.68. In particular, h r ≈ 0.055 for configuration C. The results obtained with the non-uniform radial plasma density profile will be referred to as 'configuration E'. Figure 13 shows the current density J i,I associated with the particles injected from I on the z ref = 0.1 m plane. The corresponding divergence angle isθ I = 56.5 • , about 13 • lower than in configuration C. Indeed, J i,I at the outlet's center is 2.71 times larger in configuration E than in configuration C. Consequently, a larger portion of the overall ion flux stays closer to the thruster's axis, namely the line 0 < z < 0.8 m, y = 0 m, x = −0.05 m. In configuration E, the current |I * II,e | is equal to 0.038 A, indicating a slight decrease compared to configuration C (see table 5). This reduction is likely attributed to the diminished divergence angle. Due to the lower |I * II,e |, the steady state values of ϕ I increase to 70.0 V, in accordance with the physical analysis presented in section 4. A comparison between configurations C and E is shown in figure 14 in terms of ϕ and |v i | along the thruster's axis. The potential drop and, consequently, the ion acceleration's on the z line are very similar. The specific impulse is calculated to be 650 s, exhibiting an increment compared to the value obtained in configuration C (see table 9). This result can be attributed to the reduction in the divergence angle.
Based on the findings presented in this section, it is evident that the non-uniform radial plasma density distribution at the thruster outlet has a non-negligible impact on the trajectory of ions and consequently affects the performance of the thruster. Despite that, all results are in perfect accordance with the physical analysis presented in section 4.

Appendix C
Configuration C has been simulated accounting for the collisions described in section 2.1. The mass utilization efficiency has been assumed equal to η u =ṁ i /ṁ = 0.98 for both thrusters [92], whereṁ i the ion mass flow rate at the thruster outlet. The plasma potential ϕ is depicted in figure 15 and the ion plasma density n i in figure 16. Specifically, a comparison against the results obtained in section 4 for the collisionless case is reported along the z and x axes for the parameters ϕ, n i , and n e . Collisional and collisionless results are similar. Thus,