XML Generation Notice This is an automatically generated XML version of the publication: "PIC/Fluid simulations of the plasma expansion in a planar magnetic arch - https://doi.org/10.34657/31251", originally published under the license CC BY 4.0 International. The XML version was generated on 28-08-2026 by Technische Informationsbibliothek (TIB). In spite of high standards in the development of the script, we cannot outrule the possibility of mistakes in the process of file conversion. The authors are not responsible for possibly existing mistakes or differences in this XML-version. The link above directs you to the authorised version of the publication. If you spot mistakes in the XML-version, please contact us at oarepository@tib.eu. This XML-Version of the publication "PIC/Fluid simulations of the plasma expansion in a planar magnetic arch - https://doi.org/10.34657/31251" is provided "as is", without any guarantee (“ausdrücklich oder stillschweigend”). TIB does not take responsibility for potential mistakes, missing passages and consequences resulting from this use. Hinweis zur XML-Generierung Dies ist eine durch die Technische Informationsbibliothek – TIB am 28-08-2026 automatisiert erstellte XML-Version der Publikation "PIC/Fluid simulations of the plasma expansion in a planar magnetic arch - https://doi.org/10.34657/31251", die unter der Lizenz CC BY 4.0 International veröffentlicht wurde. Trotz sorgfältiger Entwicklung des Skripts können Fehler im Konvertierungsprozess nicht ausgeschlossen werden. Die Autor*innen der Originalpublikation sind für etwaige Fehler/Differenzen in der XML-Version nicht verantwortlich. Die oben verlinkte Datei stellt die autoritative Version der Publikation dar. Hinweise auf Fehler on der XML-Version nehmen wir unter oarepository@tib.eu gern entgegen. Die vorliegende XML-Version der Publikation "https://doi.org/10.1088/1361-6595/adab8e" wird "wie es ist" bereitgestellt, ohne jegliche Gewährleistung oder Garantie, sei es ausdrücklich oder stillschweigend. Die TIB übernimmt keine Verantwortung für etwaige Fehler, Auslassungen oder Folgen, die aus der Nutzung dieser Datei resultieren. xml_created 28-08-2026 dc.identifier.uri https://doi.org/10.1088/1361-6595/adab8e dc.rights.license CC BY 4.0 International dc.relation.doi https://doi.org/10.34657/31251 PIC/Fluid simulations of the plasma expansion in a planar magnetic arch 29 01 2025 PIC/Fluid simulations of the plasma expansion in a planar magnetic arch, Guaita, MAhedo, EMerin, M, , (2025-01-29). Magnetic arches (MA) (i.e. the magnetic topology that emerges when placing two magnetic nozzles with opposite polarities side by side) are an attractive option for the clustering of multiple electrodeless plasma thrusters, as they are characterized by a zero magnetic dipole moment and thus allow a reduction of perturbing magnetic forces on the spacecraft. This work employs the hybrid code EP2PLUS to simulate and study the plasma expansion for such a magnetic topology in the planar limit. First, a reference simulation is used to analyze the leading physical mechanisms that govern the plume properties. Ions are thus found to be characterized by a double peaked velocity distribution function close to the symmetry plane, where the plasma beams emitted by the two thrusters merge, while the magnetic force acting on electrons is shown to shape both the lateral confinement of the plume, and the thrust profile provided. Second, a parametric sweep on the strength of the magnetic field shows that its influence on the propulsive properties and on the characteristics of the plume saturates for values of the Hall parameter larger than around 10. Beyond this value of the Hall parameter, only the in-plane electron currents are found to be particularly sensitive both to the magnetization levels and the boundary conditions employed, although they are also largely decoupled from the rest of plasma properties. Finally, background pressure effects were considered by including collisions with neutral atoms in the simulations, highlighting the relevance of neutral entrainment in the modification of the plume properties and in the propulsive performance of the MA. IntroductionElectrodeless plasma thrusters (EPTs) are a rapidly emerging technology in the electric propulsion community because of their simplicity, scalability, and absence of plasma-immersed electrodes [1,2]. This family of thrusters most notably 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.includes the helicon plasma thruster [3][4][5] and the electroncyclotron resonance thruster [6][7][8], in which the plasma generation is guaranteed by exciting the electron population through electro-magnetic waves. A key component of these architectures is the magnetic nozzle (MN) [9][10][11], used to externally accelerate the plasma and increase the provided thrust. A traditional MN consists of a diverging magnetic field B, applied by the thruster's magnetic circuit, which confines the plasma in the perpendicular direction and accelerates it axially. In fact, in the common case of an 'electron-pressure-driven' MN (where the electron temperature is much larger than the ions') the orthogonal electron pressure gradient, combined with the E × B drift, gives rise to an azimuthal and diamagnetic electron current j θ e . This generates a magnetic force density that pushes electrons axially, and by reaction contributes relevantly to the total thrust provided by the EPT [12,13]. To maintain a current free plume in the presence of the axial magnetic force acting on electrons, an ambipolar potential drop must also develop, confining electrons and generating an increase in the axial momentum of ions equivalent to the magnetic thrust provided [9,14].A single MN produces a non-zero magnetic dipole moment that can give rise, in presence of terrestrial or planetary magnetic fields, to disruptive torques for the control of a spacecraft's attitude. However, with the aim of easily scaling the generated thrust and providing a form of thrust control, the idea of clustering multiple EPTs has been recurrently proposed in the literature [15][16][17][18][19]. It becomes therefore intuitive to envision the pairing of EPTs with opposite magnetic polarities, so as to cancel out the total dipole moment [20]. In such a configuration the gradually diverging magnetic topology of a MN is modified into a magnetic arch (MA), in which a portion of the magnetic field lines emitted by one of the EPTs curves and connects with those of the other [15]. These field lines are termed closed magnetic lines, while the outer lines of each EPT, which do not connect with those of the neighbouring one and close back on themselves are termed open magnetic lines. These definitions also allow to identify a separatrix, the magnetic surface that lays on the boundary between open and closed lines, whose position and geometry depends on the details of the MA topology.The magnetic topology of a MA implies important variations in the plasma response with respect to the traditional MN. Amongst these the most relevant are the coupling between the two plasma sources, which now result directly connected by the arching magnetic lines, and the interaction and merging of the two plasma beams emitted by each EPT [17]. Additionally, in a MA the plasma must eventually move across some closed magnetic lines in order to expand downstream. While ions typically result unmagnetized because of their large mass, electrons must either move along the open magnetic lines present at the edges of the plume or undergo magnetic detachment from the closed ones to reach infinity [16] and grant a current-free plume. The failure to do so would eventually cause an electrostatic deceleration of ions and a reduction in the propulsive efficiency of the MA. This causes a design trade-off in the choice of the strength of the magnetic field. In fact, at lower magnetic flux densities the lateral confinement of the plasma decreases, and the propulsive force generated by the diamagnetic electron currents also drops. On the other hand, as the magnetic field is increased, electrons injected on one of the closed magnetic field lines will find it more and more complex to reach infinity, eventually decreasing the propulsive performance [21,22].Although these differences with respect to traditional MNs are important, their investigation and treatment in the literature remains for the moment very limited. As a first approach to understand the physics of a MA, Merino et al [15,16] explored an expansion in a planar MA (constituted by two infinitely elongated plasma sources) by means of a quasi-neutral, fullfluid approach, that considered perfectly magnetized and polytropic electrons and a single and cold ion fluid. This allowed a first description of the general plasma and propulsive properties of a MA, highlighting the presence of a magnetic drag force in the region where the magnetic lines become perpendicular to the plasma flow, and the consequent importance of detachment in the propulsive performance. Indeed, it was found that a relevant role was played in this direction by the magnetic field induced by the electron currents in the plume, that caused an opening of the magnetic lines and an effective detachment of the plasma from the applied magnetic field. This model, however, also had important restrictive assumptions in the physical description that it offered. First, the ions injected from the two EPT sources, were modeled by means of a single, cold fluid, which, as this work aims to show, misrepresents the complex double-peaked velocity distribution function (VDF) that develops at the intersection between the two plasma beams. Second, collisional events were fully neglected, so electrons were perfectly magnetized (i.e. the Hall parameter was infinity). Thus, electrons were forced to remain within the initial magnetic surface and only energetic electrons emitted on open magnetic lines could flow downstream and cancel the ion current to form a current-free plume. A non trivial point was that conditions imposed at the downstream boundary of the domain had to be compatible with that electron behavior.The case of a 3D MA, created by two circular plasma sources, has been first simulated by Di Fede [17] and Cescon [23], who have investigated multiple architectures for EPT clusters by means of a particle in cell (PIC) approach, highlighting for the case of the MA a strong exchange of electron currents between the magnetically linked EPT sources. With respect to the case of two MNs, this was found to decrease the total electron current reaching infinity, and consequently reduced the total ambipolar potential drop that developed in the plume. Despite the relevance of those works, various open points were left pending. First, collisions were neglected in the principal study cases (marginal collisions were commented in an appendix, having no effect on the solution), thus having electrons fully magnetized, exactly as in [16]. Second, in order to reduce the computational time, the vacuum permittivity was increased by a factor of 45 2 , thus increasing the plasma non-neutrality and inverting the hierarchy of the Debye length with the electron Larmor radius. Some studies have suggested that in the case of magnetized plasmas a strong artificial scaling of the vacuum permittivity, such as the one employed in these works, may alter the plasma response and reduce the reliability of the results [24,25].Finally, a first experimental investigation of a MA was carried out by Boyé et al [26,27], and proved the feasibility of the concept by detecting the flow of plasma across the closed magnetic lines as well as the theorized reduction of the plume divergence with respect to a MN configuration. This investigation also noted, in accordance with the previous numerical studies [16,17], a reduction in the ion velocity in the far plume when the MA configuration was employed instead of two distinct MNs.None of the publications mentioned so far have provided any insight into the effects of testing a MA inside experimental facilities, which in the case of MNs are responsible for relevant differences between the plasma response obtained in a vacuum chamber and in space. These effects are predominantly linked to the larger background pressure of neutral gas found in a testing facility, which can lead to: (i) the modification of the plasma density map and an increase of the plume divergence due to resonant charge exchange (CEX) collisions between ions and neutrals [28], (ii) an increased perceived mass flux due to the entrainment of ionized background neutrals in the plume expansion [29][30][31], (iii) a reduction of the ion acceleration due to collisional drag forces [32] and the electron cooling caused by inelastic collisions [29,30], and (iv) an increased diffusion of electrons across magnetic lines caused by collisional events [33]. Additionally, it has been suggested that the finite extension of the plume expansion permitted in a vacuum chamber, and the presence of conductive facility walls, may alter the map of electron currents [34] and in some cases strongly influence the value of the magnetic thrust provided by a MN [35].This work aims to investigate the open points identified in the literature with regards to the MA and its application for clusters of EPTs. Namely, we attempt here to further detail the physical processes governing the expansion of a MA plume, with particular attention on the interaction between the distinct plasma beams produced by each EPT, and the transport of plasma downstream across the closed magnetic lines. We then assess the influence of a variation in the magnetic field strength on the expansion properties, as well as facility effects linked to the vacuum chamber size and the collisional processes involving the neutral background population.For such a parametric study, a physically accurate, but computationally undemanding model was required. The planar version of the in-house, hybrid code EP2PLUS was therefore employed, allowing through the reduced dimensionality fast but meaningful simulations of a planar MA. The hybrid formulation presents three advantages for the present problem. First, the PIC formulation for heavy species constitutes a more adequate approach for ions with respect to the single, cold ion fluid model. Second, the hybrid formulation allows to satisfy quasi-neutrality directly, keeping the Debye length as the smallest length of the plasma and limiting non-neutral effects to the Debye sheaths around surfaces. Third, the electron fluid model offers a computationally-efficient and lownoise way of obtaining electron currents and maps of the electric potential in a weakly collisional regime. Still, a polytropic law is used here for the closure of the fluid equations, meaning that the electron thermodynamics can only be captured through an empirical and approximate relationship with the plasma density. Nevertheless, this simplified description of the electron thermodynamics and the use of a planar study case are well within the scope of this work, that aims to enrich the fundamental understanding of the physical processes and qualitative relationships that govern the plasma properties in a MA.The 3D version of EP2PLUS has been presented, described, verified and validated in numerous previous publications. Amongst these we point to: [36] for a general overview of the code and a detailed description of the particle-wise algorithms and collisional operators; to [37,38], respectively, for the derivation and numerical solution of the fluid equations for magnetized electrons; and to [39,40] for examples of the verification and validation of the code. The planar version of EP2PLUS is a trivial mathematical and numerical simplification of the 3D version.The remainder of this article is structured as follows. The main characteristics of EP2PLUS and the modeling strategy are first explained in section 2. Then the settings of the simulations presented in this work are detailed in section 3, and section 4 describes the results of a reference simulation, employed to showcase the general plasma response when expanding through a MA and the resulting propulsive properties. An analysis of the effect of the magnetic field strength on the plume expansion is then presented in section 5, focusing on how its influence on the electron population translates into variations in the observed propulsive performance. Next, in section 6 we discuss the implications of the choice of boundary conditions (BCs) and domain size on the plume characteristics, examining also the effectiveness of the BCs introduced in [16] in this more complex simulation scenario. Finally, a simulation case with a more complete description of the neutral population in the plume, and of the collisional processes that it entails, is described and analysed in section 7. A summary of the main conclusions of this work is carried out in section 8. ModelEP2PLUS simulates the electrons by means of a weaklycollisional drift-diffusion fluid model, and heavy species through a PIC description. At each time-step of the code, the heavy species and electron modules are solved sequentially, with the former providing fundamental inputs to the latter through the hypotheses of quasi-neutrality and of a current-free plume. The computational loop is then restarted by advancing the heavy particles in the newly obtained electric field. Each module can be configured with several options and alternative models. An overview of their main features are described in the following subsections, while section 3 states which modeling options are used in each of the simulations presented in this work. Heavy species moduleAt each time-step the heavy species macro-particles are advanced employing a standard series of PIC algorithms. First, macro-particles are injected into the domain from appropriate injection surfaces (see section 3.1.2). Second, a Boris-Cylrad Leapfrog scheme [41] is employed to integrate the particle's equation of motion according to the local electric and magnetic fields, linearly interpolated to the particle's position. Third, boundary effects are evaluated for those particles that have crossed important surfaces (i.e. domain boundaries, or material walls). Finally, the macro-particles are weighed to the nodes of the structured mesh, through a linear weighing scheme [42], yielding the densities, fluxes and temperatures of each population. These are then passed as inputs to the Electron Module. The heavy species module can be used to simulate only ions, but may also include neutrals either as a prescribed background map or as PIC macro-particles. Another optional feature of this module is the inclusion of collisional effects on heavy species such as ionization and CEX collisions, evaluated after the particle motion and before their weighting to the mesh nodes. As this option is only used in section 7, the details of the collisional models used are deferred to that point. Electron moduleFor a quasi-neutral plasma, and a quasi-stationary (∂/∂t = 0), inertialess electron fluid characterised by a diagonal and isotropic pressure tensor, the electric current continuity and the electron momentum equations respectively read∇ • j = 0,(1)and− ∇ (nT e ) + en∇ϕ + j e × B − R e = 0.(2)Here most of the notation is standard, so that n = n i = n e represents the quasi-neutral plasma density, j = j i + j e indicates the total electric current density obtained as the sum of the ion and electron currents, and T e the scalar electron temperature in energy units. Moreover, e is the elementary charge, ϕ the electric potential, B the applied magnetic field vector and R e the collisional force acting on the electron fluid. For all the simulation cases presented here, the electron fluid equations are closed at the level of the momentum through the polytropic lawT e = T e0 n n 0 γ−1 ,(3)where both the polytropic coefficient γ and the constants T e0 and n 0 are calibrated empirically. Although such an empirical closure represents an important simplification of the electron thermodynamics, multiple experimental [43][44][45] and kinetic [24,46,47] studies have shown polytropic parameters in the range γ = 1.1-1.24 to qualitatively capture the electron cooling that takes place in the bulk of MNs sufficiently well for the purposes of this work. Indeed, the detailed study of the electron thermodynamics in a MA is out of the scope of this paper, and we simply employ a constant polytropic index γ = 1.2 in all the simulations presented here. The polytropic law allows to define a thermalized potential [48]Φ = ϕ + γT e0 e (γ − 1) 1 − n n 0 γ−1 ,(4)such that the 'generalized electric force' en∇Φ groups the pressure and electric forces on electrons. Numerically, expressing the momentum equation in terms of the thermalized potential,en∇Φ + j e × B − R e = 0,(5)results advantageous as it allows to reduce numerical cancellation errors [49] between these two counter-acting forces. Incidentally, in the case of a collisionless plasma (R e = 0), equation (5) would yield a balance between the magnetic force j e × B and the generalized electric force en∇Φ, so that Φ would be constant along the magnetic lines [16]. Typically, the plasma expansions in MNs and MAs are weakly collisional (e.g. in the collisional simulations of section 7 the collisional force is two orders of magnitude smaller than the magnetic and generalized electric forces), so that we here retain the collisional term in equation ( 5), although we only expect mild deviations from the collisionless limit just described.The general model for R e implemented in EP2PLUS isR e = S s=1 ν es m e n (u e − u s ) = B χ (j e + j c ) ,(6)where m e is the electron mass, u e = −j e /en the bulk velocity of the electron fluid, the suffix s is used to discern between quantities relative to the different heavy species present in the simulation (be they modeled as macro-particles or as a fixed background), and ν es is the momentum collision frequency between electrons and the heavy species s. In the second equality of equation ( 6) we have also introduced the dimensionless Hall parameter χ = eB/(m e ν e ), with ν e = s ν es the effective electron momentum collision frequency, and an equivalent current density j c = (en/ν e ) s ν es u s that groups the collisional effects from the heavy species.With the definition of these variables, the electron momentum equation can be solved for the total electric current density yielding the generalized Ohm's law [37]j = −K • n eB ∇Φ + j c χ + j i ,(7)where, if written in a magnetically aligned reference frame {1 ⊥ , 1 y , 1 ∥ } (with 1 ∥ = B/B, 1 y the versor pointing out of the simulation plane, and1 ⊥ = 1 y × 1 ∥ ), K = χ 1 + χ 2   1 −χ 0 χ 1 0 0 0 1 + χ 2  (8)is the normalized conductivity tensor. Of course, in the practical implementation of the code, this conductivity tensor is appropriately rotated to the employed Cartesian frame{1 x , 1 y , 1 z }.A first option in the electron module is to simplify equations ( 6) and ( 7) by prescribing a map of the Hall parameter χ in the domain, and dropping the negligible effect of j c (in the more complete simulation of section 7 we obtain eBj c /(χ n|∇Φ|) = O(10 −3 ) throughout the domain) [36]. This effectively decouples the collisional effects on the electron fluid from the solution of the heavy species, and affords for cleaner parametric analyses. This option is indeed followed in most of the simulations presented in this work, as described in the next section. Alternatively, one can retain the full collisional model for electrons, and provide the necessary momentum collision frequencies for electron-neutral and electron-ion collisions. This is the approach employed in section 7, using a set of collision frequencies reported there.Once the appropriate BCs are provided (see section 3.1) the system of equations ( 1) and ( 7) is solved for j e and Φ with the staggered Finite Volume scheme presented in [38]. The new potential map, required to further advance the ion macroparticles, is then retrieved through equation (4). Simulation setupAs mentioned, the simulations presented in this work are structured so as to reproduce and expand on the problem analyzed in [16]. A total of seven different simulations have been run, with the following modeling choices and settings:1. Reference Simulation (H30): employed to showcase the general characteristics of the plasma expansion in a MA. It makes use of the simpler version of the model described: neutrals populations and collisional effects on ions are omitted, and the collisional force on electrons is modeled by dropping j c and assuming a fixed and uniform map of χ = 30. In fact, we find that it aids the discussion to first analyze the plasma response in this simplified physical scenario, where the complexity of the numerous collisional processes that take place in a plasma plume is reduced. 2. H3, H10, H150: These three simulations are used to evaluate the effects of the magnetic field strength on the plasma plume. They employ increasing magnitudes of the magnetic field and the simpler version of the model. From one simulation to the other the fixed and uniform maps of χ are varied proportionally with the variation of the magnetic field. 3. Big domain: used to elucidate the effects of the downstream BCs on the simulation results. It uses exactly the same simulation settings as the case H150, but with a doubled domain size (and cell count) in both dimensions. 4. In space: proposes and alternative set of BCs more suitable for infinite in-space expansions in the case of high electron magnetization. It uses settings identical to the H150 simulation, except for the downstream BCs for electron currents, which are detailed in section 6. 5. Background pressure (BGP): proposes a more accurate study of the collisional effects on the plume structure. This is done by employing the more complex version of the model described: two neutral populations (a prescribed background one and a PIC one) and their collisions with ions are included, and the more general expression of the electron collisional force of equation ( 6) is used. The details of the settings for neutrals and collisional events considered in this simulation are provided in section 7.Table 1 summarizes the general simulation settings assumed for each case, along with the uniform values of χ employed in the simulations making use of the simpler version of the model.The domain investigated in our simulations is defined by a right-handed reference frame with the Cartesian vector basis oriented so that the Oxz plane coincides with the 2D expansion plane, and the Oz axis points downstream. As a result, the Oy axis lies on the out-of-plane direction of the 2D expansion, along which the plasma properties are assumed infinitely uniform. Two plasma sources, emitting along 1 z , are located on the Oxy plane and span respectively from x = −6 cm to x = −4 cm and from x = 4 cm to x = 6 cm, so that Oyz represents the simulations' symmetry plane. The applied magnetic field B is generated by four infinite electric wires positioned respectively in x = −7, −3, 3, 7 cm, and z = 0, each carrying the same electric current I along 1 y . The signs of the electric currents flowing along the wires are anti-symmetric about the Oyz plane, so that the total magnetic moment of the system is null, and the magnetic field at the center of each source is B 0 .Figure 1 shows a sketch of the considered magnetic topology, computed by the analytical superposition of the Biot-Savart laws for each wire, and the simulation domain in the reference simulation H30, which corresponds to half of the physical domain described thanks to symmetry conditions. Note that the strength of the magnetic field, which is varied for some of the simulations of table 1, scales linearly with the current passing through the wires (e.g. for B 0 = 100 G I = 500 A, while for B 0 = 300 G I = 1500 A). Finally, the reference point for the polytropic closure is fixed at the center of the injection surface (x 0 = 5 cm, z 0 = 0), where we also impose Φ = ϕ = 0. The reference plasma density and electron temperature are n 0 = 10 18 m −3 and T e0 = 5 eV. Domain boundariesAs mentioned, the solution of the electron fluid equations ( 1) and ( 7) requires the imposition of BCs on all of the domain edges. The unknowns of these equations are the thermalized potential (which replaces the more physical but numerically less stable electric potential) and the electron currents. The BCs can therefore be applied either on the value of Φ (Dirichlet conditions) or, more commonly, on the local normal electric Additionally, particles crossing any domain boundary must also be correctly treated depending on the boundary type. Figure 1 details the type of boundaries employed in this work, and the consequent BCs imposed on the electron fluid and on the the particle populations for the simulation cases H3, H10, H30, H150 and Big Domain are described in the following subsections. The minor deviations that the In Space and BGP simulations undergo with respect to these settings are detailed in their corresponding Sections (respectively sections 6 and 7). Symmetry boundary.In order to impose symmetry conditions, a null perpendicular current (j n = 0) is imposed as BC for the fluid equations, while all particles crossing such boundaries are specularly reflected. Note that under these conditions j ni = j ne = 0, so that equation ( 7) yields 1 n • ∇Φ = 0 as expected on a symmetry plane. Injection boundary.The plasma injected into the domain by the EPT source is composed of singly charged Xenon ions and electrons. The ion particles are injected by sampling through an acceptance-rejection algorithm [50] the forward portion of a Maxwellian distribution function [36]. The ion Maxwellian VDF is characterized by a Gaussian density profile centered in x 0 = 5 cm, a uniform temperature profile, and a sonic velocity c s = (γT e + T i )/m i :n i (x, 0) = n 0 10 −3[(x−x0)/R] 2 , (9) u i (x, 0) = c s (x, 0) 1 z , (10) T i (x, 0) = 0.1 eV. (11)Note that, with these injection conditions, ions at the center of the injection surface have an average sonic Larmor radius r li = m i c s (x 0 , 0)/(eB 0 ) = 1.5 cm, which is expected to rapidly decrease as they accelerate downstream and the magnitude of the magnetic field drops. Taking the half-width of the injection surface (R = 1 cm) as a characteristic dimension of our problem, ions result only marginally magnetized near the source and rapidly demagnetize as they expand downstream. Finally, any particle returning to the injection surface from inside of the domain is simply removed from the simulation. Regarding the electron fluid, instead, it is assumed that in the quasi-neutral bulk of the EPT source an equilibrium between the lateral electron pressure and magnetic forces exists (i.e. ∂ϕ/∂x = 0) [51,52]. This, from equations ( 2) and ( 9), results in a Dirichelet condition to be imposed on the thermalized potential,Φ (x, 0) = γT e0 e (γ − 1) 1 − 10 −3(γ−1)[(x−x0)/R] 2 . (12)Note that, from equation ( 4) and in absence of numerical noise, this condition indeed recovers a flat electric potential profile across the injection surface,ϕ (x, 0) = 0. (13) Dielectric boundary.A dielectric condition is imposed along all of the upstream boundary, with the exception of the injection region. Since a neutral particle population is omitted in all simulations except the BGP case, and therefore ion recombination at material walls is not simulated, all particles crossing this surface are removed from the simulation. From the point of view of the electron fluid, a Neumann condition is imposed on Φ by locally requiring that j ne = −j ni , with the latter value of the normal ion current obtained from the heavy species module. Incidentally, note that to assure the consistency of our simulations a case in which this upstream boundary was modelled as a floating metallic object was also run, and the simulation results were found to be essentially insensitive to the BCs imposed at this surface. Chamber walls.In the commonly adopted experimental scenario in which the vacuum chamber and the thruster body are electrically decoupled (i.e. a floating thruster configuration) a zero net current flowing to the facility walls should be obtained at steady-state. In order to enforce this condition, let ϕ W be the potential of the metallic chamber, and let us call Q the edge of the Debye sheath that develops between the quasi-neutral boundary of the domain and the chamber walls. The sheath is assumed electron-confining, planar, collisionless, and characterized by an undrifting Maxwellian electron distribution. Then, according to classical sheath theory, the electron current density crossing the sheath is [53]j neQ (ϕ W ) = −e n exp e (ϕ W − ϕ) T e T e 2π m e Q ,(14)and imposing that the total electric current collected by the chamber walls is zero one has with I iW the total ion current flowing towards the wall, computed from the heavy species module. Equation ( 15) is implicit for the wall potential ϕ W , and is therefore solved with an iterative Newton-Raphson scheme at each time-step. Once ϕ W is known, equation ( 14) provides the electron current density to be imposed as BC in the electron module.I iW + ˆAout j neQ (ϕ W ) dA = 0,(15)Regarding the heavy species, all particles crossing these boundaries are removed from the simulation, analogously to what is done in the case of the dielectric boundaries. Numerical settingsThe domain of the reference simulation is shown in figure 1. This domain is discretized by means of a uniform mesh composed of cells of dimensions ∆x = ∆z = 0.83 mm. To guarantee the stability of the model the simulation time-step is required to respect a CFL-like condition ∆t < ∆z/(5c s (x 0 , 0)) [42]. Where a safety factor 5 has been considered to ensure the global compliance to the constraint. Regarding the weight of the ion macro-particles, a population control algorithm [36] was used to guarantee the presence of at least N mp,min = 10 ion macro-particles per cell in the less dense regions of the domain (meaning that more that 100 macro-particles per cell are present in the denser parts of the plume), leading at the end of the reference simulation to more than 2 × 10 6 ion macroparticles in the simulation domain. The steady-state condition of all simulations was confirmed by verifying that both the number of macro-particles in the domain, and the value of the vacuum chamber potential ϕ W varied by less than 1% over the last 700 steps. A convergence analysis was performed on all numerical parameters discussed here (e.g. cell size, particle number, time-step, and simulation length), and no appreciable changes were observed in the simulations by refining them beyond these values. The most relevant simulation settings, shared among all simulations in this work, are summarized in table 2. For specific settings, refer to table 1. Reference simulationThe panels of figure 2 show the steady-state maps for the reference simulation case (H30) of the thermalized potential Φ, electric potential ϕ, plasma density n and in-plane ion velocity ũi = u i − u yi 1 y (with u yi = u i • 1 y , note that hereafter we will employ a tilde to indicate vectors thus projected on the simulation plane). To aid with the ensuing discussion, figure 3 also displays the profiles of Φ, ϕ and n along the symmetry boundary of the simulation (x = 0) and the axis x = 5 cm of one of the injection surfaces. Figure 4, instead, presents the same plasma properties at axial sections z = 0, 10, and 20 cm. The curves corresponding to the reference simulation analyzed here are shown in green in the two latter plots. Incidentally, through the polytropic closure, the electron temperature is simply proportional to the plasma density and is therefore not shown in these figures; for completeness, however, we note that it ranges approximately from 5 eV at injection to 2.5 eV at the downstream boundary. The analysis of the steady-state properties of a strongly magnetized plasma expansion (χ ≳ 10) is most easily tackled by starting from the map of the thermalized potential. As mentioned, in fact, equation (5) shows that for weakly collisional plasmas the thermalized potential Φ remains roughly constant along the magnetic lines. This can be ascertained by observing the magnetic lines overlain in white in figure 2(a). As a consequence, in this regime, the value of Φ in the domain is mostly dependent on the chosen BCs, and is decoupled from the remaining plasma properties. The inverted Gaussian profile of Φ imposed at the injection surface through equation ( 12) (see figure 4(a)) is in fact mapped over the whole simulation domain by following the magnetic field lines, forming the curved valley observed in the 2D profile of Φ. A given profile of the thermalised potential imposes a strict relationship between the electric potential and the ion density through equation ( 4), which together with the ions' equation of motion determines the maps of ϕ, n and ũi .The two lobes of high electric potential at the edges of the injection surface in figure 2(b), together with the sharp axial decrease immediately in front of it (see also figure 3(d)) force ions to initially accelerate axially, following the magnetized electron fluid. As the magnetic lines begin to diverge, however, so does the plasma beam, although the large inertia of ions allows a slight diffusion of the plasma density across the field lines (see figure 2(c)). This is evidenced by the relative angles between the velocity streamlines and the magnetic lines in figure 2(d). The strong initial axial acceleration combined with the outward curving of the magnetic lines, in fact, causes ions to detach inwards as they advance downstream [10]. This portion of the domain is the near-plume, in which the two ion beams proceeding from the distinct plasma sources are yet to interact, and the expansion in front of each EPT is independent of the other. As a consequence, the plume's behaviour here is analogous to the one observed for a traditional MN, which is typically characterized by an inward ion detachment and a strong axial drop in the plasma density and electric potential [9] (see figure 3, second column).Starting from z ≈ 6 cm, however, the ion beam begins to merge at the symmetry plane with the one emitted from the specular source (here represented by macro particles reflected at x = 0). As a result the ion density, along with the thermalized and electric potential, rises when moving downstream along the symmetry axis (see figure 3, first column). In the lateral Plasma properties for the reference case (χ = 30). The magnetic field lines are shown in white, with the outermost ones corresponding to the lines intersecting the extremities of the injection surface, while the magnetic separatrix is highlighted in red. Note that, for visualization purposes, all plots are saturated at the extremities of their color-maps. direction, on the other hand, these quantities still display two local maxima (see figure 4, second column, and consider the presence of a symmetry plane in x = 0), indicating that the coalescence of the beams is yet under way. The plume's behaviour thus drifts away from the one of a MN, signaling the entrance into the transition-region. Here, the magnetic lines begin to curve back towards the symmetry plane to close into the specular EPT source, and figure 2(d) shows that the majority of ion streamlines are now forced to cross the magnetic lines with a negative detachment angle α = ∠(ũ i , B) (i.e. the angle between the ion velocity and the magnetic field vector) in order to continue their axial motion downstream. Nevertheless, a region with positive detachment angles persists beyond the magnetic separatrix [shown in red in figure 2(d)], which divides open and closed field lines in the top left portion of the domain.As the ion beams begin to merge at z ≈ 6 cm, the lateral velocity of the total ion population reduces to zero near the lower edge of the simulation domain. This results in the sharp drop observed in figure 2(d) in the magnitude of the ion velocity and the swift axial swerve of their streamlines around x = 0, z ≈ 6 cm. Note that, conversely to what is observed in traditional MN geometries or unmagnetized plumes, the axial potential profile from this point downstream increases rather than dropping, as highlighted by figure 3(c). Consequently, the ions moving downstream along the symmetry axis experience a slight deceleration. This is caused by two concurring effects. The first is the closed magnetic topology of the MA, characterized by magnetic lines orthogonal to the axial direction in the proximity of the symmetry plane, which have the effect of reducing the axial mobility of electrons. To maintain quasi-neutrality, therefore, the axial force acting on ions must also decrease. The second is the symmetrical geometry, which causes a rise in the plasma density where the ion beams meet (see figure 3(e)). Under the polytropic assumption this also implies an additional confining contribution of the pressure forces acting on the electron fluid, which must again be balanced by the electric field. Mathematically this effect is obtained by the inversion and axial derivation of equation ( 4) to yield the axial electric fieldE z = − ∂Φ ∂z − γ T e0 en 0 n n 0 γ−2 ∂n ∂z . (16)While in most conventional plumes the second term in this equation provides a positive contribution as the plasma expands, in this case it can have a negative effect due to the increase in density observed along the symmetry plane.Likewise, the axial profile of the thermalized potential displays an initial sharp drop followed by a gradual axial increase due to the mapping of the injection BCs to the points of the symmetry plane (refer to figure 3(a)). It is this later, gradual rise in Φ that, together with the corresponding increase in the plasma density, causes negative values of E z and a deceleration of the ion population. It is interesting that, although at the outer edges of the plume the ions continue to accelerate downstream, the effect just described acts on the denser portion of the plasma bulk and on the axial component of the ion momentum, and therefore has, as will be seen later, a negative impact on the provided thrust.Starting from the local maximum in the plasma density along the symmetry plane (at x = 0, z ≈ 12.5 cm) the two plasma beams can be considered fully merged into a unique far-plume, and all properties present a single maximum in the profiles along x (see the right-most panels of figure 4). At the outer edges of the plume the ion streamlines continue to diffuse with a positive detachment angle across the field lines, similarly to what happens in a MN. At the center, instead, the magnetic lines become roughly parallel to the x axis, forcing the ions to cross them orthogonally in order to expand downstream. Soon, however, the electric forces acting on the ions become small enough to be neglected, and their streamlines result essentially straight.The general plasma response in a planar MA detailed here provides results in qualitative accordance with the study of a 3D MA carried out in [17]: in both works one may identify the local minimum at the symmetry plane of the electric potential in the near-plume, the plateau in the ion velocity in the transition-region, and the complete merging of the plasma beams in the far-plume. Nevertheless the different dimensionality of the problem (together with the other distinctions in the models employed) causes some quantitative discrepancies between the two configurations. At an axial distance z/R ≈ 25, in fact, the plasma density has dropped roughly by a factor 10 2 in the present planar case and 10 3 in the 3D case. Similarly, the potential drop from the injection surface at z/R ≈ 25 is e∆ϕ/T e0 ≈ 2.5 for the planar MA and e∆ϕ/T e0 ≈ 4 in the 3D one. Interestingly, though, the total potential difference between the thruster and the vacuum chamber walls does not display marked differences between the two cases, with the planar case yielding eϕ W /T e0 = −4.56 while the 3D simulations give eϕ W /T e0 = −5.1. The total potential drop that develops in the plume, in fact, does not depend on the dimensionality of the expansion, but rather on the confining energy required to control the outflow of electrons and impose a current-free plume. As a consequence, ϕ W depends strongly on the gradual cooling of electrons, here modeled with a polytropic parameter γ = 1.2 in good accordance with the electron thermodynamics observed in [17].The differentiation of the thermalized potential map allows, through equation ( 7), the retrieval of the out-of-plane and inplane electron current densities (j ye = j e • 1 y and ȷe = j e − j ye 1 y respectively). These are displayed in figure 5, together with the in-plane total current density ȷ = ȷi + ȷe . Note that, considering the conductivity matrix written in a magnetically aligned reference frame in equation ( 8), and keeping in mind that all out-of-plane gradients are null, the out-of-plane currents result a factor χ larger than the in-plane ones. This is confirmed when comparing the scales of panels (a) and (b). Indeed, the first panel displays a large diamagnetic loop of electron currents (exiting the simulation plane in the bottom, red part of the plot, and entering it in the upper, blue portion), which, because of the negligible ion velocities in the out-of-plane direction, can be considered equal to the total out-of-plane current density j y ≃ j ye . Note that the curvature of the axis of the diamagnetic current loop (i.e. the line j ye = 0), turns towards the symmetry plane so that it is always aligned with the externally imposed magnetic field. This generates the in-plane magnetic forces which are fundamental for the lateral confinement of the plasma plume and, as will be discussed shortly, for the propulsive performance of a MA.Figure 5(b), on the other hand, shows clearly that the value of the Hall parameter χ = 30 considered here is insufficient to force the electron current to precisely follow the magnetic lines. Indeed, numerous electron streamlines form large angles with the magnetic field in order to reach the downstream domain boundary. Nevertheless, the effects of magnetization may be appreciated by noting that the magnetic lines that lead to the domain corners (highlighted in the plots by the thicker white lines) carry larger current values with respect to the ones found in their surroundings (e.g. consider the higher current intensity surrounding the magnetic line to the bottom-right domain corner). This is due to the larger mobility of electrons moving along magnetic lines rather than across them, and to the presence of the symmetry plane, that imposes a null outflow of currents through the bottom domain boundary. As a consequence, slightly reduced values of the electron currents are observed in the inner region of the MA, between the magnetic line to the bottom corner and the symmetry plane (labeled Region B in the plot), where electrons must diffuse across the magnetic field to leave the domain through the downstream boundary. The flow of electron currents is instead freer in the portion of the domain found between the two field lines to the domain corners (Region A), where electrons can run along the magnetic lines to reach the downstream boundary. This effect is even more prominent in the plot of the total current density. The total currents in Region A, in fact, flow mostly upstream, while in the rest of the simulation domain (Regions B and C) the total currents tend to flow downstream. Indeed, the increased mobility of electrons flowing parallel to the netic lines in Region A has the effect of imperceptibly increasing the electron flux inside it, and causes the upstreaming total currents observed there. In Regions B and C, on the other hand, the electron mobility along z is reduced by the magnetic field mainly directed along x, so that the electron current drops below the magnitude of the ion current-unaffected by the magnetic field and therefore flowing straight in the axial direction-and forms a region of downstreaming total current.As mentioned, the out-of-plane diamagnetic currents ultimately provide the magnetic force crucial for the propulsive performance of MNs. Indeed, the accumulated axial magnetic thrust (per unit of y-length) generated in a control area Ω by a planar MA isF m (z) = − ˆΩ j y B x dΩ. (17)Taking advantage of the domain symmetry to consider also the plume in front of the second EPT (i.e. located at x < 0 and not represented in the plots), Ω(z) is defined as the rectangular region spanning from x = −15 cm to x = 15 cm and from z = 0 to a given z coordinate. Figure 6(a) plots the accumulated magnetic thrust thus computed, expressed as a fraction of the total plasma momentum flux across the injection surfaces (F 0 = 2 ´x0+R x0−R n(x, 0) T e (x, 0) + m i u 2 zi (x, 0) dx = 12.5 mN m −1 ). Once again the results obtained for the reference case described here are shown in green. The magnetic propulsive force displays an initial sharp rise, corresponding to the first portion of the plume where the intensity of j ye is largest. As the field lines bend towards the symmetry plane, however, the total magnetic thrust provided drops slightly, showcasing a region of negative axial magnetic force j y B x (i.e. magnetic drag), and the resistance of the magnetic field to the axial expansion of the plasma plume. This is a natural consequence of the closed topology of a MA, and is also highlighted by the deceleration of the ion population close to the symmetry axis discussed above. Indeed the axial coordinate at which we observed the decrease in the ion velocity corresponds roughly to the one of the maximum cumulative magnetic thrust. These observations highlight the necessity of plasma detachment in the far plume of a MA, to avoid the null total magnetic thrust that would ensue from a perfect plasma magnetization and the continuous decrease of the accumulated thrust beyond z = 6 cm. Finally, note that the computed magnetic thrust displays small gradients a the downstream domain boundary (z = 25 cm), indicating that the domain dimension is adequate for the discussion of the propulsive properties of this MA, and that the regions beyond the domain boundaries only provide marginal contributions to the total value of F m . η div (z) = Ėkin,z (z) Ėkin (z) ,(18)computed as the ratio between the flux of axial ( Ėkin,z ) and total ( Ėkin ) kinetic energy through the lateral and downstream boundaries of the control area Ω(z). Here the subdivision in near-plume, transition region and far-plume is evident once again. The divergence efficiency in fact decreases in the nearplume, analogously to what happens in a MN because of the portion of lateral acceleration that ions are subjected to. However, in the transition region, η div rises again because of the symmetry conditions that cause an increase of the ion density close to the symmetry plane, where their bulk velocity becomes essentially axial. Finally, once in the far-plume, the merging of the ion beams has fully taken place, and their lateral expansion causes once more a decrease of η div . It is important to mention that the results obtained in [16] for a simulation case analogous to the one presented here displayed in all ion properties a shock-like structure that propagated from the symmetry plane. The main difference between that work and the present one is, indeed, the model adopted for the ion population, which is here described through the PIC approach as opposed to the single fluid formulation employed there. By simulating the ions with a kinetic approach, the formation of a double-peaked ion VDF is observed in the region close to the symmetry plane, where the two ion beamlets deriving from the distinct plasma sources merge. Each peak of the VDF represents the ions pertaining to one of the two sources. Indeed, figure 7 displays the histogram of the ion VDF (i.e. the number of particles per unit volume and area of phase-space) measured at x = 0, z = 6 cm, which presents two marked and symmetrical peaks at energies corresponding roughly to the potential drop that develops between the injection surface and the considered location ∆ϕ ≈ 20 V. As mentioned, this VDF presents two marked and symmetrical peaks at energies corresponding roughly to the potential drop that develops between the injection surface and the considered location ∆ϕ ≈ 20 V. Although, taken individually, each ion beam only has a temperature of 0.25 eV and can be roughly considered to be monoenergetic, their coalescence results in a VDF with a spread in velocities corresponding to a temperature T i ≈ 5 eV. Such a distribution cannot be represented with a single and cold ion fluid as done in [16], causing the shock-like artefact observed in that work. Effect of the magnetic field strengthHaving presented and analysed the general properties of a plasma expanding through the magnetic topology of a MA, we now investigate how this response is influenced by the strength of the applied magnetic field. We therefore consider three additional cases (H3, H10, and H150), identical to the reference one except for the intensity of the magnetic field map, which is modulated to yield B 0 = 30, 100, and 1500 G at the center of the injection surface. Under the assumption that the electron collisionality remains unvaried between simulations, the uniform map of the Hall parameter is also proportionally scaled to χ = 3, 10, and 150. Figures 3 and 4, along with figure 6, already displayed the axial and lateral profiles of the plasma properties and magnetic thrust obtained for these new values of χ. It results immediately evident that the cases with χ = 10, 30 and χ = 150 yield very similar results, and the Hall parameter has to be decreased to χ = 3 for relevant changes to develop. In fact, the value of χ affects the employed electron model only through the conductivity tensor introduced in equation ( 7). As the Hall parameter increases, the conductivity tensor has the effect of decoupling the electron currents perpendicular to the magnetic field from the gradients of Φ through a factor 1/(χ 2 + 1). This essentially inhibits the diffusion of electrons perpendicular to the magnetic field, and confines their motion along the magnetic lines, impeding the lateral expansion of the plume under the effect of the electron pressure gradients. The magnitudes of ȷe , however, are delimited throughout the domain by the BCs employed and the continuity equation. As a consequence, because of the large values of the parallel conductivity σ e = χ n/(eB) that characterize the weakly collisional plasma typical of a MA, the parallel gradients of Φ must result negligible in order to produce finite electron currents along the magnetic lines. As we approach the case of perfect magnetization and an infinite Hall parameter, then, the thermalized potential must remain constant along the magnetic lines, so that its map in the domain will only depend on the values imposed at the injection surface and on the magnetic topology [9]. As a consequence, when χ > 10 changes in the Hall parameter will only cause small variations in the map of Φ, and in the general plasma response. Keeping this in mind, and for the sake of conciseness, the present section employs mainly the comparison between the cases with χ = 3 and χ = 150 to showcase the effects linked to variations in the electron magnetization.Figure 8 compares the general plasma properties obtained for these two cases. The first thing that jumps to the eye is the noticeable differences between the maps of the thermalized potential. As for the reference case, in the case χ = 150 the profile of Φ imposed at the injection surface remains practically constant along the magnetic field lines, and is mapped over the domain by following them. On the other hand, when χ = 3 the electron currents may move more freely across the magnetic lines (allowing larger currents to form in the domain) and close at domain boundaries which are instead shielded by the magnetic field for χ = 150. Steeper gradients of Φ can therefore form along the magnetic lines, and the inverted Gaussian profile of Φ imposed at injection through equation ( 12) gradually dissipates along the domain (see figures 3(b) and 4(b)). This parallel smoothing of the gradients of Φ also allows a reduction of the cross field gradients of the thermalized potential, so that the map of Φ in the far plume results nearly uniform and flat. The first effect that the larger cross-field mobility of electrons obtained for χ = 3 has on the map of ϕ is the relevant reduction of the two lobes of high electric potential at the edges of the injection boundary, caused by the decay of the lateral gradients of Φ. These lobes, as already commented in section 4, had the effect of laterally confining the ion beam in order to follow the magnetized electron fluid. As a result, both the maps of the plasma density and ion velocity depict a less dense and more divergent plume for the case with χ = 3 than with χ = 150. The potential rise that was observed in the reference case after the merging of the ion beams is also present and enhanced for χ = 150, but is instead shifted upstream and flattened out when lowering the Hall parameter to 3 (see figure 3, first column). In fact, although the axial increase in density along the symmetry plane is still present in both simulations, it moves upstream for the less magnetized and focused plume. Moreover, the axial rise in the thermalized potential observed in the more magnetized cases from z = 8 cm downstream, and which contributed to the negative axial electric field in this region of the domain, is absent for the case with χ = 3 because of the smoothing of the map of Φ obtained at low magnetic fields. Observing the map of the ion velocity, this results in a monotonic axial acceleration of ions along the symmetry axis of the domain for χ = 3, and a more complex profile for χ = 150. In this latter case, in fact, ions are strongly accelerated by the sharp potential drop present just in front of the plasma injection, but are then mildly slowed down close to the symmetry plane (note the axial potential rise after z = 6 cm in figure 3(c)). The result is that at the downstream boundary of the domain ions have roughly reached the same outflow velocity for both cases, and noticeable variations in the ion velocity are only observed at the upper simulation boundary.The last panels of figure 8 show the out-of-plane electron currents, and the corresponding magnetic force density j ye B directed perpendicularly to the magnetic field (here positive or negative respectively if aligned or opposed to 1 ⊥ ) and responsible for the magnetic thrust and the lateral confinement of the plume. It is interesting to note here that the magnitude of j ye displays a non-monotonic variation as the strength of the magnetic field increases. Between χ = 3 and χ ≈ 10, in fact, the plume undergoes a noticeable magnetic focusing under the effect of the rising magnetic force, reducing the lateral expansion and increasing the density in the plume bulk. This larger density and the smaller divergence of the plume (compare the χ = 3 and χ = 10 lines in figures 3(f) and 4(h)) are responsible for the increase in the magnitude of j ye observed between figure 8(i) and figure 5(a). As the strength of the magnetic field further increases (i.e. χ > 10), however, the plume does not display further relevant focusing, as shown by the density plots in the last rows of figures ( 3) and ( 4). The confining magnetic force j ye B necessary to balance the expanding action of the electron pressure (shown in figures 8(k) and (l)) is therefore roughly constant for χ > 10, and can be obtained with smaller out-of-plane currents as B increases. The magnitude of j ye is therefore observed to decrease as χ moves from 30 to 150 (see figures 5(a) and 8(j)).Nevertheless, because of the higher plasma densities found at the symmetry plane for larger values of χ, the propulsive performance of the MA is expected to increase as the strength of the magnetic field rises. Figure 6(a), in fact, shows a noticeable improvement in the magnetic thrust when moving from the mildly magnetized case (χ = 3) to higher magnetizations. Note however, that the imparted magnetic force does not scale proportionately with the strength of the magnetic field, and that it appears that already at χ = 10 the maximum magnetic thrust of this MA has almost been reached. A similar scaling of the magnetic thrust with the magnetic field was obtained by Andrews et al in the PIC simulations of a traditional MN [24]. The trade off between the strength of the magnetic field forming the MA, and the thrust obtained from the thruster under design should therefore take into account the non-linear dependence between these two values.Interestingly, the initial slope of all the thrust curves in figure 6(a) results identical. This is caused by the BCs applied at the injection surface, as demonstrated by retrieving the axial derivative of F m by inserting the expression for j y obtained from equation (7) into equation (17):dF m dz z=0 ≈ − ˆLx 0 B x B χ 2 (1 + χ 2 ) en 0 dΦ (x, 0) dx dx,(19)which yields approximately the same value for χ >> 1. The trends of the magnetic thrust with χ are further confirmed by taking into account the plot of the divergence efficiency. Figure 6(b), in fact, shows how η div increases under the effect of a stronger magnetic confinement. In the case of a higher χ, in fact, the plasma density results higher close to the symmetry plane, where most of the ion kinetic energy flows axially; n instead drops steeply at higher values of x, where the ion kinetic energy has a more relevant x component (see panels (h) and (i) of figure 4). As a consequence the energy flux becomes mostly axial, yielding the large values of η div . Conversely, for lower values of χ, we observe in figure 4 smaller densities close to the symmetry plane and a less pronounced drop in n as one moves along x, justifying the decrease of the divergence efficiency.Finally, the first two plots in figure 9 show the electron inplane current densities obtained for χ = 3 and χ = 150. In the former case, the map of ȷe exhibits two current loops leaving the plasma source near its edges, and returning back to the source along its center. The existence of these loops raises the overall value of ȷe in this region, but does not contribute to the net electron current delivered by the source into the plume. This structure, which was already present in the case χ = 30 of figure 5(albeit much smaller in importance and extent), is explained as follows. First, we remark again that, for sufficiently high values of χ, the strong perpendicular gradients of Φ imposed at the source are successfully propagated downstream along the magnetic lines, with little diffusion or smoothing. However, as B is lowered and χ becomes sufficiently small, the damping of the perpendicular gradients of Φ as we progress from the source and into the domain becomes more marked. Eventually, this mechanism leads to a nearhomogeneous value of Φ early in the expansion. Second, this smoothing of the map of Φ means that a non-negligible variation of this quantity takes place along the magnetic lines for small values of χ. We observe that, in the central part of the source exit, the parallel component of ∇Φ points upstream, while in the periphery of the source exit it points downstream. Third, since the parallel electron conductivity σ e is large, the existence of these parallel gradients drives the strong electron currents into and out of the source, respectively. Our numerical experiments prove that this is a robust mechanism whenever σ e is sufficiently large, as in this regime Φ results essentially governed by χ alone and is unaffected by σ e , whereas the parallel electron currents are proportional to the parallel gradient of Φ and σ e .Analysing now the rest of the domain in figure 9, it is apparent that for the low magnetization case χ = 3, the electron fluid is well capable of crossing the magnetic lines freely, reaching the outer boundaries without being noticeably influenced by the magnetic lines. In the χ = 150 case, instead, the largest portion of the electrons is forced to follow the magnetic tube that reaches the downstream boundary. The guiding effect of the magnetic lines, already commented in section 4, is in fact further enhanced in figure 9(b), causing the streams of high electron current density running along the magnetic lines to the domain corners, highlighted once again by the thicker white lines. These streams of high current density are also clearly visible in the lateral section of the in-plane currents presented in figure 9(e), where it is evident that they coincide with the streamlines to the domain corners, located on the vertical dashed lines labeled 1 and 2. As will be shown in the next section this is an effect liked to the domain size and geometry, which however only influences the map of ȷe . The strong magnetization obtained with χ = 150 also causes the reduced current densities found in the inner region of the MA, between the symmetry plane and the magnetic line to the bottom domain corner, where electrons attempting to move downstream experience a mobility across the magnetic lines decreased by a factor (χ 2 + 1). Boundary effectsIt is clear from the plots of the in-plane electron currents described above that, at high values of the Hall parameter, the domain corners and in general the shape and size of the domain affect the solution for ȷe (which is only weakly coupled to the rest of the plasma solution). If the simulation domain is dictated by the dimensions of a vacuum chamber then the influence of the domain boundaries on the simulation results is physical: electron currents will preferentially follow those magnetic lines that intersect with the downstream walls, which may vary for differently shaped and sized facilities. However, if the simulation aims to reproduce the infinite expansion that takes place in space, the BCs imposed at the domain edges should not represent the plasma sheath that develops between the quasi-neutral bulk of the domain and the facility walls, but a global downstream matching layer [54] that expresses through a reduced model the plasma expansion between the domain boundaries and infinity. In this case, therefore, the applied (numerical) BCs should not influence electron currents. For this reason we dedicate the present section to the analysis of the effects that the domain boundaries have on the simulation results.Firstly, we note that for lower values of the Hall parameter (χ = 3) the chamber BCs employed up to here hold also for an infinite expansion. Indeed, because electrons in this case can freely cross the magnetic lines when expanding to infinity the expression for a global matching layer results analogous to the expansion in a Debye sheath. On the other hand, at higher χ, when the magnetic field has a larger influence in the electrons' motion towards infinity, the choice of appropriate BCs becomes more challenging. In order to assess the effects of the domain boundaries in this case, we repeat the case χ = 150 but with the domain size doubled in both directions in the 'Big Domain' simulation case. The resulting electron currents are shown in figure 9(c), where the dashed red box indicates the original domain size. It is clear that once again the electron currents follow the magnetic tube that intersects the downstream boundary, highlighting a dependency of the currents on the domain size and shape. Indeed, marked differences are evident when comparing the electron currents inside of the dashed red rectangle with those of figure 9(b). In the case of the larger simulation domain, in fact, the extended symmetry plane imposes a null current over a bigger portion of the magnetic lines connected with the injection surface, pushing the electron flux tube to those outer magnetic lines that arch beyond the bottom boundary. As stated, this is to be expected when simulating vacuum chambers of different dimensions, but represents a numerical artefact when attempting to reproduce an expansion to infinity.The outflow BCs introduced in [16] for the electron fluid represent a possible solution to this problem for specific case of a MA with large values of χ. As already mentioned, in fact, two categories of field lines may be identified in this magnetic topology: closed lines connecting the two plasma sources through the symmetry plane, and open lines circling back towards their original source in the opposite direction. For perfectly magnetized electrons (χ → ∞), the zero outflow of current imposed at the symmetry plane causes all closed field lines to carry null electron currents. As a consequence, a condition on the normal electron current j ne = 0 can be applied on all boundary faces lying on closed filed lines, even those that do not cross the symmetry plane inside of the simulation domain. In order to obtain a current free plume, therefore, all electron current must flow through the boundary crossed by the open field lines. As before, equation ( 15) is therefore used to recover the value of ϕ W (now representing the potential at infinity) that grants the current free condition, and equation ( 14) yields the Neumann BCs to be applied at the cell faces composing this surface. Note for clarity that the treatment of particles remains unvaried with respect to the previous simulation cases.The currents obtained with the BCs described above are plotted in figure 9(d), where we have also highlighted in red the magnetic separatrix, that divides open and closed field lines. To aid comparison, this same magnetic line has also been traced in red in all other panels. Notice how, as we increase the domain size from panel (b) to panel (c), the beam of electron current shifts closer to the separatrix. This is even more clearly observable in panel (e), where we have plotted the magnitude of ȷe along a lateral section with z = 10 cm, and where it is evident that the peak in the current density is located above streamline 1, which as we move from the smaller domain (red line) to the bigger one (pink line) shifts towards the magnetic separatrix, marked here with a vertical red line. It stands therefore to reason that, if we were to simulate an infinite domain, the electron currents would follow the magnetic lines surrounding the separatrix to reach infinity. Indeed, already with a doubled domain size, it is clear that the electron current density map within the red dashed rectangle in panel (c) approaches the one obtained with the modified BCs [see also panel (e)].We remark that, although the domain size or BCs have an effect on ȷe , they cause extremely marginal variations in all other plasma properties, including the out-of-plane currents. The maximum relative variations in n, j ye and total magnetic Table 3. Maximum relative variation of n, jye, and Fm computed for the Big Domain and In Space cases with respect to the H150 case. The variations have been normalized with the reference density n 0 , the total injected plasma momentum F 0 , and a reference out-of-plane current j y0 = ∆Φ en 0 /(B 0 R), with ∆Φ = 22 V the variation of thermalized potential across the injection boundary. Big domainIn thrust F m (15 cm) between the Big Domain and In Space cases with respect to the H150 case are reported in table 3. For all these quantities the maximum variations are negligible, never exceeding 5%. These small variations when changing or moving the domain boundaries are due to the large values of the electron conductivity parallel to the magnetic lines. Any small variation of the thermalized potential parallel to B will in fact cause relevant changes in the electron in-plane currents (see equation ( 7)). As a consequence, the solution to the system of fluid equations is capable of complying to any BC imposed on ȷe with negligible changes in the map of Φ, and therefore negligible influence on the other plasma properties. In general, in fact, as one approaches the collisionless plasma limit, one has σ e → ∞, forcing the gradients of Φ along the magnetic lines to be identically zero in equation (7). Therefore the map of Φ in the domain becomes solely dependant on the upstream BCs imposed at the injection surface, and j ye = (en/B) ∂Φ/∂1 ⊥ also results insensitive to the downstream BCs. On the other hand, because of perfect electron magnetization, the in-plane electron currents are forced to flow parallel to the magnetic field and result completely decoupled from the rest of the plasma properties in our drift-diffusion model, since with an infinite parallel conductivity their value becomes inconsequential for the determination of Φ. Although this decoupling between the electron currents and the general plasma response may cause numerical challenges in the accurate determination of ȷe when σ e assumes very large values, it also ensures that any inaccuracies in its computation will not relevantly affect the plume's response. Finally, note that the small variation in the rest of the plasma and propulsive properties when comparing the different sized domains in table 3 is an indication of the fact that the domain sizes and BCs employed in this work are large enough to confidently allow an extrapolation of the simulation results to larger domains and to more realistic dimensions of experimental facilities. This should not entail, in fact, any changes in the plasma properties and propulsive performance, if not for the map of ȷe . BGP effectsWith the knowledge provided by the previous Sections on the plasma response to the magnetic field strength and the applied BCs, we now set up and analyze a more complex simulation scenario (termed BGP case), that aims to reproduce a MA expansion inside the environment of a vacuum chamber, and thus subject to facility effects. To this end, the simplified formulation for the electron collisionality has been dropped, and two neutral populations have been added to the heavy species module. The first population of neutral particles is characterized by a prescribed density n n1 = 2 × 10 17 m −3 , a zero bulk velocity and a temperature T n1 = 0.05 eV ≈ 500 K. It represents the background neutral population present in an experimental facility, and corresponds to a BGP of 8 × 10 −9 bar [26,27]. In order to avoid simulating the long transient times linked to the generation of this population of background neutrals, we provide it as a constant and uniform background, postulating that it represents the steady-state equilibrium solution between collisional events, plasma and neutral interaction with the chamber walls, recombination, and extraction from the chamber pumps [55]. This allows us to maintain the simple outflow BCs for macro-particles crossing the domain edges presented in section 3.1, since the physical scattering and recombination processes that would take place at material walls are already included in the hypothesized equilibrium of this background population. The second neutral population (this time modeled as particles) represents the neutrals injected by the EPT sources. For its description we assume a propellant utilization efficiency in the source η u0 = 0.43 and a uniform injection of neutrals across the injection surface [56]. Keeping the total mass flow rate constant with respect to the previous simulation cases, we therefore inject neutrals with a density n n2 = 4.1 × 10 18 m −3 , an axial velocity u zn2 = 300 m s −1 , and a temperature T n2 = 0.05 eV. The ion flux also has to be accordingly adjusted, so that the reference density in equation ( 9) now results n 0 = 4.3 × 10 17 m −3 . The simulation settings for these new populations of neutral particles essentially mimic the ones for ions: outflow BCs at the various simulation surfaces as described in section 3.1; a controlled macro-particle weight for a minimum of 20 macroparticles per cell (corresponding to a total of 1 million macroparticles at steady-state); and the same time-step as the one employed to advance the ion macro-particles and the electron fluid. The magnetic field for this collisional case is identical to the one of the reference simulation, with a magnetic induction B 0 = 300 G at the center of each source. The collisional processes that affect the electron fluid are now evaluated explicitly, rather than fixing a uniform Hall parameter χ in the domain as in the simulation cases H3, H10, H30, H150, In Space and Big Domain. In particular, in the BGP case we include the effect of electron-neutral elastic, excitation and ionization collisions, and electron-ion Coulomb collisions, on the electron fluid. These collisional processes influence the model through the collision frequencies ν es of equation ( 2). In the case of electron-neutral collisions these frequencies ν es (T e ) are obtained as a function of the electron temperature by integrating the collision cross-sections σ es (v e ) over a non-drifted Maxwellian VDF with temperature T e [36]. The collision cross-sections as a function of the electron velocity for elastic, excitation and ionization collisions are taken respectively from the Hayashi [57], NIFS [58], and Biagi [59] databases available on the LXCat website [60]. The electron-ion Coulomb collision frequency is instead obtained as [36,61]ν ei = 2 1/2 n i e 4 ln Λ 12π 3/2 ε 2 0 m 1/2 e T 3/2 e ,(20)with the Coulomb logarithm taken to be ln Λ ≈ 10. Note here that, although the inclusion of inelastic collisions would in principle affect the cooling of the electron fluid, we maintain the polytropic parameter γ = 1.2 to aid the comparison with the previous cases.Regarding the effect of collisions on the heavy species, we consider two collisional phenomena: ionization and CEX. These two are evaluated respectively by means of MCC [36,62] and DSMC [50,63] algorithms. For CEX collisions, the collision cross-sections as a function of the relative collision velocity are taken from the Miller [64] database, while the ionization frequency is taken locally equal to the one employed in the electron fluid [36] provides a more detailed description of the collisional algorithms implemented in EP2PLUS.Figure 10 presents the plasma properties for this new simulation case. Panel (a) shows that the Hall parameter in the bulk of the plume decreases from the initial values of 150 at injection to roughly 20 at the downstream boundary of the domain. In light of this, the previously employed assumption of a uniform value of χ throughout the domain results inaccurate. Nevertheless, the magnetization of the electron fluid falls within the parametric range previously analyzed, and for which the plasma and propulsive properties displayed a reduced sensitivity to the value of χ. As a consequence, the thermalized potential map for this case, shown in panel (b), results nearly identical to the one obtained for the previous cases where a constant χ > 10 was assumed. Once again, in fact, the value of Φ remains practically constant along the magnetic lines because of the large scalar conductivity obtained for this weakly collisional plasma and the low cross-field transport of electrons.The plasma density for this case, on the other hand, is expected to display variations with respect to the previous simulations because of two distinct effects. The first is the moderate propellant utilization efficiency, which implies a decrease in the injected ion flux (and therefore in the plasma density) by a factor η u0 . The second is the inclusion of collisional processes (ionisation and CEX) that cause the generation of low energy ion particles downstream of the injection surface. To better distinguish these two effects we divide ions into two separate particle lists: the first comprised only of ions produced in the EPT source and injected sonically into the simulation domain (i.e. fast ions), the second generated by CEX and ionization collisions in the plume (i.e. slow ions). The densities of these two populations are respectively shown in panels (d) and (e) of figure 10, while panel (f) shows the total plasma density that results from their combination (and the assumption of quasi-neutrality). Interestingly, in the near plume and transition region the slow ion density is a relevant fraction of the total ion density, an observation in accordance with experimental studies of MNs, which have measured double-peaked total ion VDFs (one peak for the slow ions and one for the fast ones) in the main bulk of the plume [30]. The generation of slow ions only decays in the far plume, where the decrease in electron temperature and neutral density (see panel (g)) cause a strong reduction of the ionization and CEX collision frequencies. More specifically, this population is particularly present in two locations. The first is just in front of the injection surface, where the neutral density and the electron density and temperature are highest, resulting in the maximum collision frequency found in the domain. The second is at the symmetry plane where the two ion beams first meet, because of the combined effect of the higher densities found here, and the merging of the slow ion populations that are generated in front of each source. As a consequence of the relevant presence of slow ions, the total plasma density in this case does not display the reduction by a factor η u0 that would be expected with respect to the previous cases, and indeed it results of a similar magnitude throughout the domain. This observation is representative of the neutral entrainment detected in multiple experimental studies [29,30], which can cause an overestimation of the thrust when measured in on-ground facilities. Finally, note that the majority of the total ion density remains concentrated close to the symmetry axis even for large z coordinates, likely because of the generation of ions close to the symmetry plane through ionization collisions.The differences in the plasma density described above are also responsible for the variations in the map of the electric potential shown in panel (c), with respect to what was obtained for the previous simulation cases for large values of χ. While the lobes of high potential at the edges of the injection surface, responsible for the axial confinement of ions in the nearplume, and the sharp axial potential drop in charge of their acceleration are still present, in fact, the potential well previously observed in the transition region close to the symmetry axis is here strongly reduced. This is likely due to the relevant presence of slow ions slightly upstream of where the fast ion beams emitted by the sources merge, causing a reduction of the axial pressure gradient along the symmetry axis, and a smaller contribution of the second term of equation ( 16), linked to the plasma density. As a consequence of collisions, therefore, the potential and plasma density maps displays monotonic profiles along x shortly after the near-plume, effectively shortening the length of the transition region.The potential map just described is responsible for the acceleration of both the high and slow ions downstream, whose cumulative in-plane velocity is plotted in panel (h). It results immediately evident that the total ion velocity is much lower than the ones depicted for the previous simulation cases. This is due to the large fraction of slow ions, which are generated in the plume from neutrals with a low velocity, and accelerated downstream only by the potential drop that develops from their generation point to infinity. As a consequence these ions possess a much lower velocity, and heavily reduce the average speed of the total ion population [32]. Note, however, that the fast ions produced by the source maintain a velocity map very similar to the one shown for the reference simulation, omitted here for the sake of conciseness. In general, the net axial acceleration of both these ion populations, as well as the one of neutrals through CEX and momentum exchange collisions, is representative of the magnetic thrust provided by the MA. Another notable feature of Plasma Sources Sci. Technol. 34 (2025) 015007 M Guaita et al panel (h) are the two regions of increasing lateral velocity that develop at the edges of the injection surface (the first in the region x = 6-15 cm, z = 0-4 cm, and the second in the region x = 0-3 cm, z = 0-4 cm) and move mostly along x. These are once again caused by the slow ions generated just in front of the injection, and accelerated laterally by the electric field. Note however, by observing panel (e), that the number density of these ions is largely negligible.Finally, panel (i) shows the in-plane electron currents for this simulation scenario. As may be expected by observing the axially decreasing values of the Hall parameter throughout the domain, the behaviour of the electron current gradually transitions from the one observed for χ = 150 to the one of χ = 10 − 30. The electrons near the source, in fact, initially move mostly along the magnetic tube encompassing the downstream domain boundary (once again delimited by the thicker white magnetic field lines). As the Hall parameter decreases as they proceed downstream, however, electrons begin diffusing across the field lines, and the confined current arc that was observed for χ = 150 breaks down into a more homogeneously spread out electron flux. Finally, at the downstream boundary, the electron current streamlines turn upwards in order to comply with the fluxes there imposed as BCs.Regarding the propulsive performance for the BGP case, the obtained accumulated magnetic thrust and divergence efficiency have been included in figure 6 with a black dashed line. Note, when considering the plotted thrust, that the force provided without a MA and used to normalize the plotted values reduces to F 0 = 5.11 mN m −1 in this case, because of the effect of the non-unitary utilization efficiency. Relatively speaking, the MA in these conditions results much more effective, as it provides a total magnetic thrust F m = 0.6F 0 (see figure 6), although in absolute terms the net magnetic thrust provided is actually slightly lower than the one obtained in the cases with χ ⩾ 10. This increase in the propulsive effectiveness in spite of the moderate value of η u0 employed, is due to the plasma production that takes place in the plume. Indeed, the plasma densities of the current simulation are only slightly lower than the ones of the reference case, pointing to an effective propellant utilization efficiency (when considering also the plasma production in the plume) which is considerably larger than 0.43. When taking the ratio between the total ion mass flow at the downstream domain boundaries and the total mass flow of ions and neutrals injected by the thruster, we obtain a utilization efficiency η u = 0.53. It is important to state, however, that the production of plasma in the plume is, from a propulsive viewpoint, not as efficient as its generation and confinement in the source, since it shifts the plasma production downstream of the more intense acceleration region. Note also that, for the values of BGP and utilization efficiency employed here, the majority of the electron-neutral collisions involve neutral atoms emitted by the thruster, rather than the ones composing the chamber background population, and their effects would therefore be present also in space conditions. Finally, figure 6(b) shows that the divergence efficiency for this BGP case ultimately results the highest out of all of the considered scenarios. This confirms the improved confinement of the plasma density close to the symmetry plane already observed in the density plots. ConclusionsThe in-house, hybrid code EP2PLUS has been employed to simulate the expansion of the plasma produced by two EPTs in a planar MA, and to further study the effects linked to variations in the magnetic field strength and to facility conditions. The analysis has allowed the identification of three distinct regions: (i) a near plume, where the plasma in front of each thruster displays a behaviour analogous to that of a traditional MN topology, (ii) a transition region, where the merging of the two distinct plasma beams takes place, and (iii) a far plume, in which the lateral profiles of the plume's properties become single peaked, and the distinct plasma beams have fully merged. In particular, a detailed study of the transition region unveils the complex kinetics of ions, characterized in this area by a double peaked VDF, and the presence of a region of magnetic drag due to the closed magnetic lines that arch between the two EPTs, that mildly reduces the total provided magnetic thrust.A parametric sweep on the strength of the applied magnetic field shows that only small changes are visible for values of the Hall parameter χ ⩾ 10 both in the general plasma properties and in the propulsive performance of the MA. The only quantities that show a relevant sensitivity to χ in this magnetized regime are the in-plane electron currents. These are found to be only weakly coupled to the other plasma properties for the weakly collisional plasma analyzed, and progressively result constrained to the magnetic lines that intersect the domain's outflow boundaries as the magnetization level rises. At lower values of χ, instead, the confining effect of the magnetic field begins to break down, and the plasma is capable of diffusing across the magnetic lines, increasing the plume's divergence and reducing the magnetic thrust provided by the MA.A further study on the influence of the BCs demonstrates that different domain sizes (corresponding to different dimensions of the vacuum chamber) impose different electron current paths in a strongly magnetized plasma, although this does not appreciably alter the other plasma properties of the plume. In light of this dependency of the in-plane electron currents on the domain size, an alternative set of BCs [16] was also employed and tested, and resulted more suitable for the simulation of the infinite expansion that takes place in space.Then, an analysis of BGP effects has been carried out by including in our simulations the neutral atoms emitted by the EPTs, as well as a uniform population of neutrals representing the background density found even far from the thruster exhaust in vacuum chambers. The collisional processes involving these neutral populations, and which modify the response of ions and electrons were also included. Although a propellant utilization efficiency η u0 = 0.43 was prescribed for the plasma source in this case, the thrust provided results only slightly lower than the one obtained in the previous simulations (in which η u0 = 1), because of the entrainment of the neutral atoms caused by ionization and Plasma Sources Sci. Technol. 34 (2025) 015007 M Guaita et al CEX collisions, which implies an increase in the effective value of the utilization efficiency. Interestingly, a majority of these collisions involved neutrals emitted by the thrusters, and not those populating the chamber background, so that a similar behaviour is to be expected also in space conditions under low utilization conditions.While the present work has offered an extensive analysis of the main physics in a planar MA and the sensitivity to central aspects such as the Hall parameter and the BCs, the implementation of the electron energy equation instead of the polytropic law, will provide a more consistent response of electron thermodynamics, mainly when collisional and facility effects are important [48]. Finally, we acknowledge that a big step forward would be to simulate a 3D MA, as done in [17], but with finite values of the Hall parameter, and focusing on the 3D details of the plasma detachment, the topology of the electron current lines, and the induced magnetic field. Figure 1 . Figure 1. Sketch of the reference simulation domain and of the boundary conditions applied for a planar MA. The two counter-streaming current wires and the plasma source are infinite in the y direction. The green color-map and the streamlines display the intensity and direction of the magnetic field for the reference case. Figure 2 . Figure 2.Plasma properties for the reference case (χ = 30). The magnetic field lines are shown in white, with the outermost ones corresponding to the lines intersecting the extremities of the injection surface, while the magnetic separatrix is highlighted in red. Note that, for visualization purposes, all plots are saturated at the extremities of their color-maps. Figure 3 . Figure 3. Comparison of the axial profiles of Φ, ϕ, and n for all simulation cases. The first column shows the plasma properties along the symmetry plane, the second along the axis of the plasma source. Note the changes in scales in the vertical axes of the plots of the thermalized potential. Figure 4 . Figure 4. Comparison of the lateral profiles of Φ, ϕ, and n for all simulation cases at z = 0, 10, and 20 cm. Note the changes in the scales of the vertical axes between the first column of plots and the other two. Figure 5 . Figure 5. The out-of-plane electron currents (panel (a)), and the in-plane electron and total currents (panels (b) and (c)) for the reference case (χ = 30). Figure 6 . Figure 6. Comparison of the accumulated axial magnetic thrust and divergence efficiency for all simulation cases. The magnetic thrust is expressed as a percentage of the total plasma momentum flux through the exhaust of the two EPT sources. For all simulations except the BGP case F 0 = 12.5 mN m −1 . For this latter one F 0 = 5.26 mN m −1 . Figure 6 ( Figure 6(b), instead, displays the divergence efficiency, Figure 7 . Figure 7. Histogram of the ion VDF at x = 0 and z = 6 cm for the reference simulation, obtained by discretizing the in-plane phase-space in 60 × 60 velocity bins, and averaging the ion macro-particle count in these bins over 600 simulation time-steps. Figure 8 . Figure 8. Comparison of the main plasma properties for χ = 3 and χ = 150. Figure 9 . Figure 9. Comparison of the in-plane electron currents for different simulation cases. The magnetic separatrix is highlighted with a red line, the magnetic lines to the reference domain corners (lines 1 and 2) are shown with thicker white lines, and the dashed red rectangle in panel (c) indicates the reference domain size. Panel (e) shows a lateral section of the in-plane currents. Here the lines to the corners of the reference domain (1 and 2) and big domain (1' and 2') are shown with dashed and dotted lines. The separatrix is marked with a vertical red line. Figure 10 . Figure 10. The plasma response in the BGP case. Note that the suffices fi and si indicate respectively the properties of fast and slow ions. Likewise, i and n mark the quantities of the total ion and neutral populations. Table 1 . Table 1 . Summary of the simulation cases presented in this work and of their defining characteristics. The tss column provides the simulation times. Note that for the majority of the simulations collisional effects are modelled in a simplified way by imposing a fixed value of the Hall parameter χ, while in the BGP case χ is self-consistently obtained through collisional algorithms. The reference simulation (H30) is shown with a bold font. Domain Case B 0 (G) χ (-) size (cm 2 ) N. Cells (-) tss (µs) BCs H3 30 3 25 × 15 180 × 300 500 EF a H10 100 10 25 × 15 180 × 300 500 EF H30 300 30 25 × 15 180 × 300 500 EF (reference) H150 1500 150 25 × 15 180 × 300 500 EF Big domain 1500 150 50 × 30 360 × 600 1000 EF In space 1500 150 25 × 15 180 × 300 500 Sp b BGP 300 Comp c 25 × 15 180 × 300 900 EF Table 2 . Table 2 . Summary of the most relevant numerical settings shared by all simulation cases. Parameter Value Parameter Value ∆x 0.83 mm ∆z 0.83 mm ∆t 7 × 10 −2 µs N mp,min 10 T e0 5 eV γ 1.2 AcknowledgmentThe authors would like to thank Celian Boyé for his contributions in the initial definition of the simulation setup. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 950466). ORCID iDsM Guaita  https://orcid.org/0000-0002-6579-8629 E Ahedo  https://orcid.org/0000-0003-2148-4553 M Merino  https://orcid.org/0000-0001-7649-3663 Data availability statementThe data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.5281/ zenodo.12751281 [65]. Plasmas for space propulsion Plasma Phys. Control, Ahedo, EFusion, 53, 124037 (2011). 10.1088/0741-3335/53/12/124037 Electric propulsion for satellites and spacecraft: established technologies and novel approaches, Mazouffre, SPlasma Sources Sci. Technol, 25, 33002 (2016). 10.1088/0963-0252/25/3/033002 , Regulus: a propulsion platform to boost small satellite missions Acta Astronaut, 157, (2019). Helicon-type radiofrequency plasma thrusters and magnetic plasma nozzles Rev, Takahashi, KMod. Plasma Phys, 3, 3 (2019). 10.1007/s41614-019-0024-2 Performance characterization of a helicon double layer thruster using direct thrust measurements, Pottinger, SLappas, VBoswell, Charles, RJ. Phys. D: Appl. Phys, 44, 235201 (2011). 10.1088/0022-3727/44/23/235201 , Electron-cyclotron-resonance (ECR) plasma thruster research 24th Joint Propulsion Conf. p 2916, , (1988). Optimization of a coaxial electron cyclotron resonance plasma thruster with an analytical model, Cannat, FLafleur, TJarrige, JChabert, PPackan, Elias, DPhys. Plasmas, 22, 53503 (2015). 10.1063/1.4920966 , Coupled plasma transport and electromagnetic wave simulation of an ECR thruster Plasma Sources Sci, 30, 45005 (2021). Two-dimensional supersonic plasma acceleration in a magnetic nozzle Phys, Ahedo, EMerino, MPlasmas, 17, 73501 (2010). 10.1063/1.3442736 , Plasma detachment in a propulsive magnetic nozzle via ion demagnetization Plasma Sources Sci, 23, 32001 (2014). A magnetic nozzle calculation of the force on a plasma Phys, Fruchtman, ATakahashi, KBoswell, Charles, RPlasmas, 19, 33507 (2012). 10.1063/1.3691650 Electron diamagnetic effect on axial force in an expanding plasma: experiments and theory, Takahashi, KLafleur, TCharles, CBoswell, Alexander, RPhys. Rev. Lett, 107, 235001 (2011). 10.1103/PhysRevLett.107.235001 Kinetic insights into thrust generation and electron transport in a magnetic nozzle, Hu, YHuang, ZCao, YSun, QPlasma Sources Sci. Technol, 30, 75006 (2021). 10.1088/1361-6595/ac0a48 Ambipolar ion acceleration in an expanding magnetic nozzle, Longmier, BPlasma Sources Sci. Technol, 20, 15007 (2011). 10.1088/0963-0252/20/1/015007 Preliminary model of the plasma expansion in a magnetic arch thruster (and overview of the first prototype) 37th Int, Merino, MGarcía-Lahuerta, DBoyé, CNavarro-Cavallé, JAhedo, EElectric Propulsion Conf., No. IEPC-2022-423, , (2022-06-19). , Merino, MGarcía, -Lahuerta D and Ahedo E 2023 Plasma acceleration in a magnetic arch Plasma Sources Sci, 32, 65005 (). 10.1088/1361-6595/acd476 Magnetic nozzle performance in a cluster of helicon plasma thrusters, Fede, DiManente, SComunian, MMagarotto, P, MPlasma Sources Sci. Technol, 32, 65013 (2023). 10.1088/1361-6595/acdaf2 , Preliminary analysis of the magnetic arch plasma expansion in a cluster of two ECR plasma thrusters 10th EUCASS Conf, , (2023-07-09). , Recent advances in the clustering of high power helicon thrusters AIAA Propulsion and Energy 2019 Forum, , (2019). , Motor espacial de plasma sin electrodos con geometría en U PCT patent, , (2019). Magnetic nozzle and plasma detachment model for a steady-state flow Phys, Breizman, BTushentsov, MArefiev, APlasmas, 15, 57103 (2008). 10.1063/1.2903844 , Plasma detachment mechanisms in a magnetic nozzle 47th AIAA/ASME/ SAE/ASEE Joint Propulsion Conf. & Exhibit, No. AIAA-2011-5999, , (2011). Cluster of magnetic nozzles: a fully kinetic 3d pic simulation study 38th Int, Cescon, FFede, SMagarotto, Manente, MElectric Propulsion Conf., No. IEPC-2024-415, , (2024-06-23). Fully kinetic model of plasma expansion in a magnetic nozzle, Andrews, SFede, DiMagarotto, S, MPlasma Sources Sci. Technol, 31, 35022 (2022). 10.1088/1361-6595/ac56ec , Non-maxwellian electron effects on the macroscopic response of a hall thruster discharge from an axial-radial kinetic model Plasma Sources Sci, 33, 25008 (2024). Magnetic arch plasma expansion in a cluster of two ecr plasma sources, Boyé, CNavarro-Cavallé, JMerino, MJ. Electr. Propul, , (2024). 10.21203/rs.3.rs-4412193/v1 Determining ion velocity in the interconnected plume of a cluster of two ecrts using 2d lif 38th Int, Boyé, CNavarro-Cavallé, JMerino, MLecervoisier, AMazouffre, SElectric Propulsion Conf., No. IEPC-2024-717, , (2024-06-23). , Experimental and numerical examination of the BHT-200 Hall thruster plume 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conf. & Exhibit, No. AIAA 07-5305, , (2007). Effect of background pressure on ion dynamics in an electron cyclotron resonance thruster, Wachs, BJorns, BJoint Propulsion Conf, , (2018). Background pressure effects on ion dynamics in a low-power magnetic nozzle thruster, Wachs, BJorns, BPlasma Sources Sci. Technol, 29, 45002 (2020). 10.1088/1361-6595/ab74b6 Neutral ingestion effects on plume properties of a radio-frequency plasma discharge, Caruso, N RWalker, M LJ. Propul. Power, 34, (2018). 10.2514/1.B36404 , Fully kinetic study of facility pressure effects on rf-source magnetic nozzles Acta Astronaut, 215, (2024). On plasma detachment in propulsive magnetic nozzles, Ahedo, EMerino, MPhys. Plasmas, 18, 53504 (2011). 10.1063/1.3589268 , Impact of facility electrical boundary conditions on the performance of an electron cyclotron resonance magnetic nozzle thruster 37th Int. Electric Propulsion Conf, , (2022-06-19). Ecra thruster advances: 30w and 200w prototypes latest performances, Désangles, VPackan, DJarrige, JPeterschmitt, SDietz, PScharmann, SHolste, KKlar, P JJ. Electr. Propul, 2, 10 (2023). 10.1007/s44205-023-00046-x , Hybrid 3D model for the interaction of plasma thruster plumes with nearby objects Plasma Sources Sci, 26, 125008 (2017). 3D simulations of a magnetized Hall Effect thruster plume 36th Int, Cichocki, FDomínguez-Vázquez, AMerino, MAhedo, EElectric Propulsion Conf., No. IEPC-2019-460, , (2019). Numerical treatment of a magnetized electron fluid model in a 3d simulator of plasma thruster plumes Front, Modesti, ACichocki, FAhedo, EPhys, 11, 345 (2023). 10.3389/fphy.2023.1286345 , Cichocki, FDomínguez-Vázquez, AMerino, MFajardo, PAhedo, EActa Astronaut, 187, (2021). 10.1016/j.actaastro.2021.06.042 , Magnetic nozzle and rpa simulations vs, 10, 84 (2022). , Computer Simulation Using Particles, , (1988). Clouds-in-clouds, clouds-in-cells physics for many-body plasma simulation, Birdsall, CFuss, DJ. Comput. Phys, 135, (1997). 10.1006/jcph.1997.5723 Electron cooling in a magnetically expanding plasma, Little, JChoueiri, E ;Giono, GGudmundsson, JIvchenko, NMazouffre, SDannenmayer, KPopelier, LLoubere, DMerino, MOlentsenko, GNon-maxwellian electron energy probability functions in the plume of a SPT-100 Hall thruster Plasma Sources Sci, 117, 15006 (2016). 10.1088/1361-6595/aaa06b , Plasma beam characterization along the magnetic nozzle of an ECR thruster Plasma Sources Sci, 28, 95004 (2019). , Kinetic electron model for plasma thruster plumes Plasma Sources Sci, 27, 35013 (2018). Collisionless electron cooling in a plasma thruster plume: experimental validation of a kinetic model, Merino, MFajardo, PGiono, GGudmundsson, JIvchenko, NMazouffre, SLoubere, DDannenmayer, KPlasma Sources Sci. Technol, 29, 35029 (2020). 10.1088/1361-6595/ab7088 Using electron fluid models to analyze plasma thruster discharges, Ahedo, EJ. Electr. Propul, 2, 2 (2023). 10.1007/s44205-022-00035-6 , Modeling and simulation of the plasma discharge in a radiofrequency thruster, , (2021). , Molecular Gas Dynamics and the Direct Simulation of Gas Flows (The Oxford Engineering Science Series), , (1994). Magnetized fluid electron model within a two-dimensional hybrid simulation code for electrodeless plasma thrusters, Zhou, JDomínguez-Vázquez, AFajardo, PAhedo, EPlasma Sources Sci. Technol, 31, 45021 (2022). 10.1088/1361-6595/ac64bc Parametric analysis of a magnetized cylindrical plasma, Ahedo, EPhys. Plasmas, 16, 113503 (2009). 10.1063/1.3262529 The bohm criterion and sheath formation, Riemann, KJ. Phys. D: Appl. Phys, 24, (1991). 10.1088/0022-3727/24/4/001 Analysis of the electron downstream boundary conditions in a 2D hybrid code for Hall thrusters 37th Int, Domínguez-Vázquez, AZhou, JSevillano-González, AFajardo, PAhedo, EElectric Propulsion Conf., No. IEPC-2022-338, , (2022-06-19). Simulations of ion thruster plumes in ground facilities using adaptive mesh refinement, Korkut, BLevin, DTumuklu, OJournal of Propulsion and Power, 33, (2017). 10.2514/1.B35958 Plume characterization of a waveguide ecr thruster, Inchingolo, MMerino, MNavarro-Cavallé, JJ. Appl. Phys, 133, 113304 (2023). 10.1063/5.0138780 , Nifs database, , (2019-10-23). , Lxcat: an open-access, web-based platform for data needed for modeling low temperature plasmas Plasma Process Polym, 14, 1600098 (2017). , , , (1995). A two-dimensional hybrid model of the Hall thruster discharge, Parra, IAhedo, EFife, J, Martínez-Sánchez MJ. Appl. Phys, 100, 23304 (2006). 10.1063/1.2219165 The direct simulation Monte Carlo method, Garcia, Alexander, AComput. Phys, 11, (1997). 10.1063/1.168619 Xenon charge exchange cross sections for electrostatic thruster models, Miller, JPullins, SLevandier, DChiu, YDressler, RJ. Appl. Phys, 91, (2002). 10.1063/1.1426246 , Data from: PIC/Fluid Simulations of the Plasma Expansion in a Magnetic Arch Zenodo, , (2024).