Plasma acceleration in a magnetic arch

When two magnetic nozzles with opposite polarity are placed side by side, a ‘magnetic arch’ (MA) is formed, which connects the field lines of each nozzle into a closed-line configuration. The plasma expansion and acceleration in this magnetic topology are relevant for clusters of electrodeless plasma thrusters, as well as novel, non-cylindrical thruster architectures. A collisionless, quasineutral, two-fluid model of the plasma expansion in a MA, is introduced. The plasma properties (density, electron temperature, electrostatic potential, ion velocity, electric currents) in the 2D planar and zero plasma-beta limit are analyzed, and the magnetic thrust density is discussed. It is shown that the ions coming out of the two nozzles meet on a shock-like structure to form a single beam that propagates beyond the closed lines of the applied magnetic field, generating magnetic thrust. A small magnetic drag contribution comes from the final part of the expansion. The plasma-induced magnetic field is then computed self-consistently for non-zero plasma-beta expansions, showing that it stretches the MA in the downstream direction and helps reduce that drag contribution. Finally, the limitations of the present model are discussed.

When the plasma consists of warm electrons and relatively cold ions, as in the devices listed above, the MN is termed 1 http://ep2.uc3m.es. * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 'electron-driven.' These MNs work by perpendicularly confining the expansion of the warm plasma electrons, which must be well magnetized. This confinement occurs thanks to the applied magnetic field B a and the diamagnetic azimuthal electron current density j θe that forms as a consequence of the existence of a perpendicular electron pressure gradient and the E × B drift. This current density gives rise to a magnetic force density in the plasma. Part of this force density is directed radially inward (j θe B az ), while the other part is axially outward (−j θe B ar ). The former balances the electron pressure (and the electric force) in the radial direction; the latter gives rise by reaction to the magnetic thrust, which is the force felt by the thruster magnetic circuit due to the magnetic field induced by the plasma. In turn, the parallel electron pressure is balanced by the self-consistent ambipolar electrostatic field that forms in the MN. This field confines electrons and accelerates ions, converting the electron thermal energy into directed kinetic ion energy [5].
Downstream, the plasma jet must eventually separate from the turning magnetic lines to prevent the increase of plume divergence and the cancellation of thrust [6]. It should be noted that, at least for hot-electron and cold-ion plasmas, ions do not need to be magnetized for the MN to operate as intended; indeed, a high ion magnetization is generally undesirable, as it makes plasma detachment occur farther downstream, increasing plume divergence angle, and promoting the appearance of a paramagnetic azimuthal ion current density j θi in the plasma that results in magnetic drag [5]. However, special devices, such as the variable specific impulse magnetoplasma rocket [18], rely on the expansion of hot ions, where ion magnetization is a necessity.
A single cylindrical EPT creates a magnetic dipole moment that may induce secular torques on the spacecraft in the presence of the geomagnetic field. Flying EPTs in pairs with opposite magnetic polarities, such that the net dipole moment cancels out, is a straightforward and natural way to avoid this issue. Also, the use of more than one thruster (known as 'clustering') is a simple way of scaling thrust levels for larger space missions. A pair of EPTs has the additional benefit that, if each unit can be throttled independently, some degree of thrust vector control can be achieved without moving parts. In this configuration, the two MNs interact and their lines connect, resulting in a new magnetic topology that here we term 'magnetic arch' (MA), sketched in figure 1(a).
Similarly, the MA is an intrinsic part of some novel EPT geometries, such as the magnetic arch thruster (MAT) concept, where the cylindrical discharge chamber of traditional EPTs is replaced by a 'C'-shaped chamber, enveloped by coils that create a magnetic field essentially parallel to the walls, as represented in figure 1(b) [19,20]. By removing the rear wall that exists in cylindrical EPTs and ensuring full magnetic shielding of the remaining walls, it is hypothesized that this geometry could bring advantages with regards to losses, while reducing the appearance of external magnetic torques on the spacecraft.
The plasma expansion in an MA is radically different to that in an axisymmetric MN: while in a single MN the plasma flux is roughly parallel to the applied field B a (at least before detachment is well under way), in an MA the flux is only parallel initially; downstream, where the lines of the two MNs connect, the plasma flux must necessarily traverse the applied field roughly perpendicularly. Also, while the plasma currents in the axisymmetric MN are predominantly diamagnetic (i.e. thrust producing), they are expected to be diamagnetic and paramagnetic in the upstream and downstream regions of the MA plasma expansion, respectively. Relatedly, while in a MN the plasma-induced magnetic field B p plays a secondary role in deforming the shape of the lines, increasing divergence minimally if the MN is well-designed [21], it can play a more important role in the MA, potentially changing the line topology of the total field, B = B a + B p , with respect to that of the applied one alone, B a . Finally, the interaction of the two plasma jets coming from each end of the device may lead to collisionless shock-like structures in the plume, not found in smooth MN plasma expansions [5].
The objective of this work is to present a first model of the MA plasma expansion to examine the viability of this magnetic topology for plasma acceleration and discuss its main physical mechanisms, in particular behind ion acceleration and , which also features an MA in the plasma expansion region. In red, the location of the ionization chambers; in black, the position of magnetic coils; and in blue, selected magnetic lines. magnetic thrust production. By examining the zero plasmabeta expansion, β = µ 0 nT e /B 2 a = 0, and expansions with small β ̸ = 0, the effects of the plasma-induced magnetic field on the shape of the MA and the generated thrust are discussed. The model is of application to clusters of two cylindrical EPTs and to the novel MAT configuration described above. Finally, we identify the main physics currently outside of the present model that should be included in the future. Nevertheless, the major limitations of the study can be already stated from the outset: firstly, we shall only study a 2D planar version of the MA, rather than a full 3D geometry. Secondly, we shall ignore plasma kinetics, and employ a collisionless multi-fluid plasma model with a simple polytropic closure for the electron pressure.
The rest of the document is structured as follows. Section 2 presents the mathematical model of the plasma expansion in the MA and describes the approach followed to integrate it numerically. Section 3 contains the results of the first MA simulation using this model in the β = 0 limit, including plasma density, ion velocity, electrostatic potential, plasma currents, and magnetic thrust. Section 3.2 then discusses the plasma-induced magnetic field for β ̸ = 0, and how its presence alters the expansion and magnetic thrust with respect to the β = 0 case. The limitations of the model and its results are reviewed in section 4. Finally, section 5 briefly summarizes the main points of this work. A preliminary version of this work was recently presented in [19].

Model
A two-dimensional, two-fluid (ions i and electrons e) model of the steady-state plasma flow in an MA is considered. The model takes the following assumptions: (1) Quasineutral, collisionless, fully-ionized plasma.
(2) Inertialess, quasi-Maxwellian, perfectly-magnetized electrons with a polytropic closure relation. (3) Cold, singly-charged ions, with arbitrary magnetization, emitted from each source exit. Moreover, ions are assumed to remain cold downstream, neglecting the effects of any shock-like discontinuities on ion temperature/distribution that may exist in the solution. (4) Planar-symmetric geometry, as an intermediate step toward the actual three-dimensional geometry of the device. We consider the meridian plane of the plume and assume an infinite plasma with uniform properties in the perpendicular direction.
To normalize the model, we select the ion mass m i , the electron charge q e , and the radius R of one of the plasma thruster exits. And, using the properties at the center of one of the two symmetric thruster outlets (where variables are marked with subindex 0), the electron temperature T e0 (in energy units) and the plasma density n 0 used for injection. Note that, even if flux coming from one outlet ever arrives at the other, T e0 and n 0 are defined from the single-beamlet injection properties. This point of the outlets is also chosen as the origin of the electrostatic potential, where ϕ 0 = 0. In the following, all symbols are already appropriately dimensionless. In particular, the dimensionless magnetic field strength at the center of the outlet, B 0 , is normalized with √ m i T e0 /(eR), and coincides numerically with the dimensionless ion gyrofrequency Ω i0 , which defines the (initial) ion magnetization degree. Figure 2 sketches the problem domain. We define a righthanded reference frame with the plane Oxy coincident with the exit plane of the plasma sources, and the Oz axis pointing downstream. The plane under study is the Oxz plane, and in the 2D expansion the plasma is infinite and uniform in the y direction. The plane Oyz is a symmetry plane, and thus only the upper half of the plane (x ⩾ 0, shown in the figure) will be simulated. Without loss of generality, B is taken to point axially downstream in this part of the MA. We introduce the Cartesian vector basis {1 x , 1 y , 1 z } and the magnetic vector basis {1 b , 1 ⊥ , 1 y }, with 1 b = B/B and 1 ⊥ = 1 y × 1 b . Both bases are right-handed and orthonormal.
The applied magnetic field B a is generated by a set of thin, infinite electric wires w, each carrying an electric current I w along the 1 y direction. The arrangement of wires and their electric currents is antisymmetric about the Oyz symmetry plane, The rest of the boundaries are free (supersonic for ions) outflow boundaries (green lines). The applied magnetic field Ba strength (colormap) and streamlines (black lines) are shown. Magnetic lines connecting with the source edges and center are shown as thicker lines. The symbols ⊗ and ⊙ are used to denote the location of electric wires generating the field, with electric current going into and out of the paper, respectively. The magnetic sepatratrix line of the applied field, given by ψ Ba = 0, is plotted in red. and the sum of the I w over all the wires equals zero. The magnetic stream function of a single wire w is given by where ρ w is the polar distance from the wire. Summing over the wire contributions we obtain the streamfunction ψ Ba of the applied field. The plasma-induced magnetic field B p has the stream function ψ Bp given by Ampére's equation, which reduces to a manifestly elliptic partial differential equation: where j y is the out-of-plane plasma electric current and β 0 = µ 0 /B 2 a0 is the β parameter at the centerpoint of the thruster outlet, already normalized with n 0 and T e0 .
The total magnetic field B is the sum of the applied and plasma-induced ones, with ψ B = ψ Ba + ψ Bp , and When β 0 = 0 the plasma-induced magnetic field B p is negligible with respect to the applied one B a , and the total field coincides with the latter. Then, equation (2) may be dropped from the model. This is the case analyzed in the first part of section 3. Note that 1 ⊥ ≡ ∇ψ B /B and, for any single-variable function f(ψ B ), The relevant collisionless fluid equations of electrons and ions are where we have already imposed plasma quasineutrality, The electrons are assumed polytropic with the law T e = n γ−1 on the whole domain, with fixed exponent γ. We note that where the equality holds for γ ̸ = 1. Observe that the relevant dimensionless sound speed is c s = √ γT e . Under the assumption of full electron magnetization, electron streamlines coincide with magnetic lines. Indeed, from equation (5) we infer that, since there are no pressure gradients nor electric fields in the (uniform) 1 y direction, there is no electron fluid velocity along 1 ⊥ . Therefore we write the electron fluid velocity as: With these premises, equation (5) becomes Integrating this equation along the magnetic lines, we find that the electron energy H e is conserved along them, The function −H e can also be understood as a thermalized potential for the electron dynamics. The out-of-plane electron velocity u ye can be computed from the map of ∇H e [5]: This u ye results from the sum of the diamagnetic (i.e. pressuredriven) and E × B drifts, which are the only first-order drifts in the problem (and indeed, they scale as 1/B). The function H e , its derivative H ′ e , and consequently u ye , can be computed from the boundary conditions at z = 0 on each magnetic line. This computation can be done a priori, i.e. before solving the rest of the plasma problem. Observe that only one value of H e may be imposed per magnetic line, and this restricts the set of valid boundary conditions elsewhere.
Lastly, we note that u ∥e does not appear in equations (5)- (7), and is effectively decoupled from the rest of the problem. Indeed, it can be computed from equation (4) and the boundary conditions a posteriori, after all other variables have been solved for. In the steady state, and for zero perpendicular electron velocity (u ⊥e = 0), this equation reduces to There are two different types of magnetic lines in the MA: inner lines that connect the two plasma sources through the symmetry plane, and outer lines that go around the upper part of the domain without intersecting it. In the β 0 = 0 case, the separatrix between these two behaviors corresponds with the magnetic line labeled by ψ B = 0 (see figure 2), inner lines have ψ B < 0, and outer lines have ψ B > 0.
In steady state, the electron current on inner lines must be zero due to the symmetry of the problem, and therefore u ∥e = 0 there. This sets an additional consistency requirement on the electron velocity boundary conditions on these lines. On the other hand, outer lines are allowed to carry electron current, and u ∥e ̸ = 0 on them in general. For a globally-current free MA, the total electron current leaving the plasma sources along these magnetic lines must equal the total ion current emitted by the sources. This aspect of the model is discussed in more detail in section 4.
The electron equations have therefore been reduced to (1) a conservation law for H e , (2) an algebraic expression for u ye , (3) a line-wise differential equation for u ∥e . Equation (11) may be regarded as the law that provides the electrostatic potential on each magnetic line as a function of the electron density and the magnetic streamline function: Introducing relation (14) into the ion momentum equation (7) to eliminate ϕ and using (12) to eliminate u ye results in the following set of differential equations for n, u zi , u xi , and u yi : In the steady state, each species admits a stream function ψ j such that ∇ψ j = −nu xj 1 z + nu zj 1 x , for j = e, i. For the magnetized electrons, ψ e is a function of ψ B . For ions, which are nonmagnetized or only partially-magnetized, streamlines may differ from magnetic lines.
The last ion equation (18) can be integrated to yield (see [5] for the analogous equation in the axisymmetric MN): where D(ψ i ) depends only on the ion stream function and can be determined from the boundary conditions at the thruster outlet.
Observe that, if u yi = 0 at injection, the ion out-of-plane velocity u yi develops only when ion streamlines separate from their initial magnetic field lines; this u yi is positive if the ion streamline detaches inwardly from the magnetic field lines (i.e. along −1 ⊥ ), and negative if separation is outwardly (i.e. along +1 ⊥ ). Nevertheless, when ion magnetization is weak (B 0 ⩽ O(1)), the last term in the ion momentum equations (the ion magnetic force) is typically small. Then, if u yi ≪ 1 initially, it remains so everywhere else, and the electron magnetic force dominates in the right hand side of equations (16) and (17).
Finally, we define the in-plane ion velocity asũ i = u zi 1 z + u xi 1 x , and the in-plane ion Mach number as M i =ũ i / √ γT e .

Numerical integration
The differential ion equations (15)- (18) are in conservative form, and can be formally written as where The equations are discretized using a discontinuous Galerkin (DG) method, which for zeroth-order polynomials coincides with the finite volume method. The main advantage of the DG approach is that it enables easily improving the accuracy of solution by refining the mesh size h and/or increasing the order of the polynomials p. After multiplying equation (20) by a test vector V, integrating in an element D k with boundary ∂D k , and using integration by parts, the following weak form is obtained: where dΩ is the area differential and 1 n dS is the outwardoriented area vector differential. Upon summation over all elements D k of the domain, the second integral must be substituted by the corresponding numerical flux integral on all internal boundaries, taking into account the jump conditions across neighboring elements: where Γ int are the internal facets of the discretization, superindices '+' and '−' indicate the values of a discontinuous variable on one side and the other side of an internal facet, with 1 n pointing toward the + side, andF is a numerical flux function.
In this work a local Lax-Friedrichs flux is chosen, given bỹ with α computed as the maximum of all eigenvalues of the normal flux Jacobian (∇F · 1 n ) evaluated in each side of the facet.
A similar treatment is applied on the external boundary facets, denoted by Γ ext , except that on those facets the + side corresponds to the weakly imposed boundary conditions. The external boundary is further decomposed into Γ in , Γ out , and Γ sym for supersonic inflow, supersonic outflow and symmetry plane boundaries respectively (see figure 2). At the inflow boundary, the Q + vector on the + side is determined by the desired inflow conditions. At the supersonic outflow boundary, the Q + vector is taken equal to Q − (i.e. the value of Q on the corresponding boundary element of the domain), and finally, at the symmetry plane the density and parallel flux components of Q + equal those of Q − , while zero perpendicular flux is imposed (i.e. nu xi = 0).
The discretized plasma problem is initially integrated in time using a third-order Strong Stability Preserving Runge-Kutta scheme, given in [22]. As initial conditions for the time integration, any gross approximation of the expected steady state flow can be used to speed up convergence. After a sufficient amount of time steps, the steady state version of the equations are solved for.
In β 0 ̸ = 0 cases, the plasma-induced magnetic field problem is integrated using the continuous Galerkin method using Lagrange elements on the same mesh as the plasma problem. The weak form of (2) iŝ where V is a test function The boundary conditions used are B pz = ∂ψ Bp /∂x = 0 at the symmetry plane x = 0 as indicated above, and B px = −∂ψ Bp /∂z = 0 at the thruster exit plane z = 0. On the rest of the boundary, and on the outside of the plasma domain shown in figure 2, a thin absorbing layer with artificial anisotropic magnetic permeability is defined following [23], a method that is equivalent to a coordinate stretching, to better approximate the transition to infinity of B p . Notwithstanding this, it must be noted that the electric currents inside the thruster discharge chambers and beyond the simulation domain (i.e. further downstream) also affect the value of B p ; as these currents are unknown, and to partially mitigate their influence on the results, the peripheral part of the domain is cut off from the shown results whenever β 0 ̸ = 0. Then, the self-consistent plasma and magnetic field solutions are determined using an iterative procedure: The plasma flow is first solved for the B = B a (i.e. ignoring the plasmainduced field). This yields a first approximation to the outof-plane plasma electric current j y , that is used to compute ψ Bp for the next iteration using (2). The plasma flow is then recomputed for the new total field B = B a + B p , and this process is repeated until plasma variables and ψ Bp vary less than a prescribed tolerance from iteration to iteration, at which point convergence is reached.
The numerical implementation of the model employs GMSH [24] and FENICS [25] as open-source building blocks. The code has been verified successfully by simulating two simple cases: (1) a plasma flowing in a straight, uniform magnetic field; and (2) a 2D planar MN, and comparison against the existing DIMAGNO code [5]. Mass and momentum are successfully conserved in the simulation. A convergence study with mesh size and polynomial order was also conducted and confirmed the correct behavior of the code.

Simulation results
The applied magnetic field used for the simulations presented in this section is generated by four identical wires contained in the Oxy plane, located at x = 3, 7, −3, −7, as shown in figure 2. The thruster outlet in this half of the MA is located on the Oxy plane and goes from x = 4 to x = 6, and the normalized magnetic field at the center point of the thruster outlet, (z, x) = (0, 5), is B 0 = 1 (mild initial ion magnetization).
The boundary conditions at the thruster outlet are modeled as follows: i.e. the plasma density profile is assumed Gaussian, centered on x = 5 falling three orders of magnitude at the edges of the outlet, and the axial velocity is given as a function of the local sound velocity such that the in-plane ions are sonic at the magnetic throats (M i0 = 1) [5]. In the last expression, U ∥e is a constant electron parallel velocity imposed on outer magnetic lines (i.e. those that can carry electron current), computed such that the net electric current emitted by the plasma source is zero.
The electron polytropic exponent is set to γ = 1.2, a value commonly found empirically in MNs [26]. The out-of-plane electron velocity at the outlet of the thruster is given by equation (12) evaluated at z = 0 and the conditions above, This choice of boundary conditions means that, at the outlet, there is no x-directed electric field, and the electron pressure and the electron magnetic force are in balance. For a smoother numerical solution and to prevent regions of zero plasma density, we extend these conditions all the way from x = 3 to x = 7, where the thin wires that generate the magnetic field are located. This completely determines the value of H e on the whole the domain.
A mesh with cell diameter h = 0.29 and elements of order p = 1 are used to obtain the plasma solutions shown below, while in section 3.2, elements of order p = 2 are used for the computation of the plasma-induced magnetic field.

Plasma expansion in the β 0 = 0 limit
We begin with the analysis of the plasma expansion when β 0 = 0, and the total magnetic field is B = B a .
The map of H e (ψ B ) plays a crucial role in the plasma response, as its derivative H ′ e fixes the out-of-plane electron velocity u ye , which in turn defines the electron magnetic force. The resulting profile of H e and the u ye that follows from the boundary conditions are plotted in figure 3. The direction of the gradient of H e causes the electron out-of-plane velocity u ye to be positive and negative below and above the magnetic centerline of the plasma outlet, respectively, resulting in a magnetic force that confines the expanding electrons to their respective magnetic tubes. This change of sign contrasts with the typical situation in an electron-driven MN, where the outof-plane electron velocity u ye has the same sign everywhere in the meridian plane [5].  Figure 4 displays the steady-state solutions for the plasma density n, electron temperature T e , electrostatic potential ϕ, in-plane ion velocityũ i , and in-plane ion Mach number M i . Several aspects of these results stand out. Firstly, the plasma expansion is initially guided by the magnetic field, and as the (essentially unmagnetized) ions accelerate, their streamlines do not adhere to the magnetic lines, separating inward with respect to B, as in the axisymmetric MN case [6]. The plasma density, electron temperature, and electrostatic potential all decrease axially as the plasma expands in this first fraction of the domain. Secondly, ion streamlines on the periphery of the MA become essentially straight. Inward detachment proceeds as in an MN for outer magnetic lines above the separatrix ψ B = 0, which curve back and around the upper part of the domain. However, for inner magnetic lines below ψ B = 0, which eventually curve downward and intersect the symmetry plane, ion trajectories must traverse magnetic lines in the outward direction. This changes ion detachment from being inward-directed to being outward-directed in part of the domain, in contrast to what occurs in an MN. Thirdly, closer to the symmetry plane, an oblique shock structure forms at the location where ion streamlines coming from the two thruster outlets meet. Ion streamlines are deflected at the shock, and plasma density, electron temperature, and electrostatic potential rise across it. In-plane ion velocity and Mach number, which increase in the first part of the expansion, fall through the oblique shock, although ions remain supersonic downstream of it.
A major conclusion arising from these results is that the unmagnetized ions are not confined by the MA magnetic field, but are able to form a jet that propagates beyond it to infinity. This observation is crucial to the validity of the MA concept and therefore for the operation of a cluster of two cylindrical EPTs with opposing magnetic polarities. Figure 5 displays the in-plane electric current density,ȷ = n(ũ i − u ∥e 1 b ), imposing a uniform distribution of electron macroscopic velocity on the outer magnetic lines (ψ B > 0) to yield a globally current-free solution. As the electron flux on inner magnetic lines (ψ B < 0) must be zero, theȷ in this region results solely from the ion current. Above the separatrix line ψ B = 0, the strong compensating electron current responsible for making the system globally-current-free dominates. Clearly, the location of the separatrix line ψ B = 0 with respect to the thruster outlet is a major defining aspect of the MA plasma expansion with regard to the in-plane electric currents. This aspect is further discussed in section 4. Figure 6 depicts the x and z magnetic force densities j y B z and −j y B x , where j y = n(u yi − u ye ) is the out-of-plane electric current density. We note that j y is dominated by the electron contribution everywhere in the domain, as ion magnetization is low. Observe that the product j ye B is essentially independent of the magnitude B 0 by virtue of equation (12). Indeed, this product depends essentially on the initial electron pressure gradient at the thruster outlets, which determines the profile of H e .
The two components of the magnetic force density are largest near the thruster exit plane. The x force density, essentially perpendicular to the magnetic lines in the first part of the domain, confines the plasma expansion laterally. As it can be observed from figure 6(a), this confining force points in the x > 0 direction in the innermost part of the arch (i.e. in the region between the two plasma sources), while it points along x < 0 everywhere else, helping reduce the divergence of the jet.
The z force density gives rise to magnetic thrust, and is seen to be large and positive at the beginning of the expansion, where n, T e and B are large. A small negative contribution exists downstream on inner magnetic lines, beginning at the point where B x = 0 and lines curve down toward the symmetry plane. This negative contribution is mostly noticeable in the region after the shock wave, where plasma density (and therefore the out-of-plane current density) increases locally again. These characteristics are consequential on the magnetic force density and the generation of magnetic thrust: while positive thrust is generated initially in the region where they resemble a traditional MN, the magnetic force generates drag in the downstream region where the magnetic lines of each thruster connect, lowering the net thrust of the device.
As follows from the sum of the electron and ion momentum equations (5) and (7), the magnetic thrust force generated by the plasma contained in a rectangular control volume Ω(z) that spans the domain from the initial plane z = 0 to a variable axial position z can be equivalently computed as where ∂Ω(z) is the full boundary of the control volume. The first integral is the volume integral of the axial magnetic force density in figure 6(b), while the second integral is the flux integral of total momentum on the boundaries of the integration domain. Observe that the relative importance of electron pressure thrust decreases to zero sufficiently far downstream, and that ion momentum dominates as the expansion converts electron thermal energy into ion kinetic energy. Figure 7 displays the thrust force F(z) normalized with F(0), the initial momentum flux of the plasma coming out of the sources (directed ion momentum plus electron thermal momentum). Positive magnetic thrust is produced initially, in the first part of the expansion. When the plasma approaches the bend in the magnetic lines and the shock, magnetic thrust plateaus, and thereafter, a minor contribution of negative thrust (i.e. magnetic drag) results by which F(z) decreases by a small amount. As indicated above, this is a natural consequence of the closed shape of the inner magnetic lines and the maps of n, u ye , which give rise to a negative axial magnetic force density in the second part of the expansion, as shown in figure 6. In the present simulation with β 0 = 0, the negative contribution decreases F(z)/F(0) about an 8% at z = 20 with respect to its maximum, which occurs at z ≃ 7. Figure 8 displays the normalized plasma-induced magnetic field B p /(β 0 B a0 ), computed from the j y current density corresponding to the β 0 → 0 limit. To mitigate the influence of the plasma currents beyond the simulation domain on the solution, the peripheral part of the results has been cut out from this and following plots.

Effect of the plasma-induced magnetic field
As noted before, j y B a is essentially independent of B a0 in the low ion magnetization regime under consideration. Consequently, by virtue of equation (2), the dimensionless group B p /(β 0 B a0 ) is also essentially independent of β 0 and B a0 . The direction of B p opposes the applied one in the proximity of the symmetry plane, and points roughly axially downstream far from it. Hence, the trend of B p is to stretch the MA downstream as β 0 increases. Figure 9 displays the self-consistent total magnetic field B = B a + B p and the streamlines for ψ B = ψ Ba + ψ Bp for different values of β 0 . As β 0 is increased from 0, B p gains relative importance. The main effect is the modification of the geometry of the central lines of the MA, which are stretched downstream. The separatrix line displaces downward, and some inner magnetic lines that intersected the symmetry plane for β 0 = 0 are converted into outer lines that go around along the periphery of the domain instead. As the fraction of inner lines shrinks and the fraction of outer lines grows, more magnetic lines can carry electron current away from the device.
For the larger values of β 0 shown, a region of very low magnetic field strength forms near the symmetry plane and the separatrix eventually intersects with it, forming an 'X' point at which B = 0, visible for β 0 = 0.08 in figure 9. This brings about a topology change of the MA, which now features a new magnetic region that forms beyond the 'X' point, whose magnetic lines are disconnected from the upstream plasma sources.
While the general characteristics of the plasma expansion are qualitatively similar to the β 0 = 0 case, the value of β 0 has a major effect on the generated magnetic thrust F(z)/F(0). Figure 7 displays the evolution of the thrust force as a function of β 0 . It is evident that, while the initial part of the curve roughly coincides for all cases, the stretching of the MA reduces the negative drag contribution that occurs downstream. Indeed, as β 0 increases, the magnetic thrust force generated within the domain rises. For β 0 = 0.04, F(z)/F(0) remains almost flat after a weak maximum, and for β 0 = 0.08, the local maximum disappears altogether, with the relative thrust gain F(z)/F(0) increasing by 10% at z = 20 with respect to β 0 = 0.
These results suggest that the plasma-induced field B p plays a central role in shaping the expansion and the propulsive performance of the device.

Discussion
The results from the previous section merit additional discussion. Firstly, the location of the separatrix line and its effect on electron currents deserves closer inspection, as the global current-free condition is an essential one that must be satisfied by any plasma thruster operating in space. It is possible to distinguish different types of MA, depending on the connectivity of the magnetic lines passing by the plasma source exit with the symmetry plane:  Interestingly, as shown in section 3.2, the downward displacement of separatrix line as β 0 is increased may cause a change of MA type. In particular, a closed arch is expected to become an open arch for a sufficiently high value of β 0 .
Secondly, even in the strict β 0 = 0 limit, collisions and an out-of-plane electric field E y are mechanisms outside of the present model that could relax the electron transport in the perpendicular direction, thus enabling u ⊥e ̸ = 0 and, as a side effect and as dictated by the continuity equation, allow u ∥e ̸ = 0 (and therefore electron current extraction from the sources) even on inner magnetic lines.
The main effect of non-zero collisions on the in-plane electron transport can be understood by including a new term in electron momentum equation (5), which now becomes where R e = nm e ν e u e is a simple representation of the collisional term. The y projection of this equation yields with χ = B/(m e ν e ) the local Hall parameter. Hence, a perpendicular electron flux arises, with u ⊥e ̸ = 0 pointing in the direction opposite to the confining −eu ye B1 ⊥ force. An electric field in the out-of-plane direction, E y 1 y , can also enable perpendicular electron flux. The E × B drift induced by this field generates a collisionless u ⊥e . This mechanism may play a role e.g. in 3D MA expansions, where E y may arise if the plasma undergoes lateral polarization, but is not present in the 2D planar geometry studied here.
Thirdly, another major aspect to be discussed is the validity of the electron model. The profile of H e that is defined at the upstream plane by the boundary conditions on n, T e , and ϕ fully determines H e in the rest of the domain, and therefore H ′ e , which dictates u ye and defines the magnetic force density in the plume. The map of H e , together with the map of n, also defines the electrostatic potential ϕ. As such, H e has a central role on the MA plasma dynamics. However, while it is reasonable to prescribe H e on the plasma-carrying magnetic lines that pass through the source, it is not evident what should be the condition on the external lines outside of this main magnetic tube, where plasma density is negligible. Here, in this first simulation, we have opted to define H e by setting n ≃ 0 and ϕ = 0 at the upstream plane for those lines. A similar problem arises in cases with β 0 > 0, if an 'X' point forms that bears a new magnetic regions beyond it, disconnected from the plasma sources, as discussed in section 3.2. At present, we have extended the value of H e on the last magnetic line before the 'X' point to this new region. Incidentally, note that if H e were constant everywhere (which can always be achieved with the right choice of ϕ upstream), there would be no u ye and hence no magnetic force on the electron fluid, and the magnetic guiding effect of the MA would disappear; this conclusion applies to traditional MNs too [5].
The electron model may need to be revisited to include inertia effects, finite Larmor radius effects, and/or more advanced closure relations based on a kinetic description, which may play a non-negligible role in some of the regions mentioned above. Altogether, these effects may modify the H e conservation law and the parallel and perpendicular transport of electrons. Similarly, the assumption of quasineutrality may need to be dropped in favor of integrating Poisson's equation, in low density regions.
Fourthly, the applicability of the 2D planar MA model to describe the actual 3D MA remains to be assessed. While it is currently expected that the planar model captures the essence of the mechanisms at play in the actual MA, adding bounds to the plasma in the third dimension can have additional effects, such as the possible set up of a polarization E y field that further changes the axial dynamics due to the E × B drift as discussed above. Bounds in the y direction also demand the closure of out-of-plane plasma currents j y , which will modify the plasma solution and the plasma-induced magnetic field with respect to the 2D planar ones. Additionally, the plasma expansion is stronger in 3D than in 2D, resulting in a faster-decreasing plasma density and electrostatic potential.
Fifthly, it is noted that the present model with cold ions neglects the ion temperature increase that is expected to occur across the shock structure seen in the solution. Including finite ion temperature and the ion energy equation in the model would be a necessary first step to study these aspects consistently. Alternatively, ions could be modeled as two distinct species, representing the ion beams coming from each plasma source, and overlapping in the interaction region. Ultimately, the correct treatment of the collisionless shock requires a kinetic model.
Lastly, we note that while the MA model remains to be validated against experiment directly, it relies on similar hypothesis to those of the MN model of [5]. That model has shown good agreement with existing laboratory measurements in the literature, in particular the plasma expansion and magnetic thrust generation [27,28], ion-inertia-driven plasma detachment [29], and the role of the plasma-induced magnetic field [30].

Conclusion
A MA is expected to form when the MNs of two EPTs with opposing polarities interact. It is also the topology of the magnetic accelerator in the novel MAT concept. A first model of the external expansion of a MA has been presented, which already describes much of its interesting plasma physics in spite of its simplifying assumptions.
The ions are seen to form a free jet that traverses the closed lines of the magnetic field, even if electrons are fully magnetized. An oblique shock structure forms when the beams coming out of the two thruster outlets meet. Electron equations reduce to algebraic relations in the inertialess, fully-magnetized, polytropic, collisionless limit of the model, and describe how the out-of-plane electron velocity is fully determined by the magnetic field map and the upstream boundary conditions.
Electron current can only be extracted along magnetic lines that do not connect the two sources through the symmetry plane (i.e. outer lines), in the limits of the model. For a given applied magnetic field map, these lines are delimited by a separatrix line; the location of this separatrix with respect to the thruster outlet determines the type of MA. In the open MA considered here, a globally-current free solution of the plasma expansion is possible. If all the lines passing through the outlet connected with the symmetry plane, no electron current could be extracted in the strict fully-magnetized, collisionless, planar, β 0 = 0 limit of the electron model, and other mechanisms would need to be included to enable it.
Net positive magnetic thrust is produced from the interaction of the out-of-plane plasma currents and the applied field. These currents are dominated by the electron contribution at low and mild ion magnetization strengths. It is observed that most of the positive contribution to thrust comes from the initial stages of the expansion, while a small negative (drag) contribution results from the region where the magnetic lines bend back to the device.
The plasma-induced magnetic field B p has been shown to stretch the MA downstream when β 0 ̸ = 0, and to increase the fraction of electron current-carrying magnetic lines that pass through the source. The stretched arch has a smaller negative drag contribution to thrust, and indeed a monotonically increasing magnetic thrust curve results already for moderate values of β 0 . Hence, modeling the plasma-induced magnetic field is essential for the correct description of the MA dynamics. This contrast with the case of an axisymmetric MN, where