Fluid-kinetic model of a propulsive magnetic nozzle

A kinetic-electron, fluid-ion model is used to study the 2D plasma expansion in an axisymmetric magnetic nozzle in the fully-magnetized, cold-ion, collisionless limit. Electrons are found to be subdivided into free, reflected, and doubly-trapped sub-populations. The net charge current and the electrostatic potential fall on each magnetic line are related by the kinetic electron response, and together with the initial profiles of electrostatic potential and electron temperature, determine the electron thermodynamics in the expansion. Results include the evolution of the density, temperature, and anisotropy ratio of each electron sub-population. The different contributions of ions and electrons to the generation of magnetic thrust are analyzed for upstream conditions representative of different thruster types. Equivalent polytropic models with the same total potential fall are seen to result in a slower expansion rate, and therefore to underpredict thrust generated up to a fixed section of the magnetic nozzle.


Introduction
Magnetic nozzles [1,2] (MNs) act as the main plasma acceleration stage of electrodeless thrusters such as the helicon plasma thruster [3][4][5] (HPT) and the electron-cyclotron-resonance thruster [6][7][8] (ECRT), but is it also an essential part of the applied-field magnetoplasmadynamic thruster [9] (AFMPDT), the variable specific impulse magnetoplasma rocket [10] (VASIMR), and other devices [11]. Understanding the plasma expansion in the guiding magnetic field of the MN is crucial to develop predictive models of the performance of these thrusters, as well as to determine the electric potential map and assess the energy and amount of plasma backflow to the spacecraft. * 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.
The focus of this work is on steady-state, near-collisionless MN plasma flows composed of hot, magnetized electrons and comparatively cold ions, which are relevant to (at least) HPTs and ECRTs. These flows are characterized by a monotonically decreasing electrostatic potential φ along the plasma expansion that accelerates ions downstream and confines the majority of the electrons, except for the most energetic ones which are free to escape downstream with the ions (see figure 1). This electrostatic field maintains plasma quasineutrality and global charge current balance, while it converts the electron internal energy into directed ion kinetic energy.
The potential map φ(z, r) depends strongly on the electron thermodynamics [12], which arise from electron kinetics. Given the low number of collisions in the plasma, electrons are typically away from local thermodynamic equilibrium, which hinders an accurate representation of the electron species by means of otherwise convenient fluid models, since a consistent closure relation for the fluid equation hierarchy requires to account for the full kinetic electron response. Indeed, in the collisionless limit, electrons develop temperature anisotropy and cool down as a result of the existence of effective potential barriers in phase space, arising from the interplay between electrostatic forces and magnetic mirror forces, which creates empty regions in the electron velocity distribution function (EVDF), as well as isolated regions populated with doublytrapped electrons that do not connect neither with the plasma source nor with infinity downstream. These results have been found analytically [13,14], numerically [15], and experimentally [16,17]. The resulting electron behavior is far from the commonly-used isothermal or polytropic models, whose theoretical justification fails in a collisionless plasma. Accounting for the correct electron response is crucial for the study of the plasma expansion and the determination of the propulsive performance of the MN.
Previous works on the kinetics of the near-collisionless electrons in a MN have relied on 1D models. Arefiev and Breizman [18] consider the combined effect of the electric field and the magnetic mirror effect on individual electrons in the MN, and distinguish between 'coupled' electrons (reflected electrons in the nomenclature of the present work) and 'uncoupled' electrons (doubly-trapped electrons here). By invoking a transient rarefaction wave at the front of the expansion they are able to compute the EVDF of doublytrapped electrons, assuming the axial bouncing of electrons at the rarefaction wave has an associated adiabatic invariant. While that model introduces many of the necessary concepts to understand the kinetic expansion of electrons in a MN, it seems to ignore the essential role of free electrons and empty regions in phase space, and thus it is unable to correctly determine the electric potential fall in the expansion, φ 0 − φ ∞ , and its relation to the current-free condition that a MN must satisfy. The 1D kinetic converging-diverging MN model of Martínez-Sánchez et al [13] takes free electrons into account, which eliminates the problems listed above. This model establishes the basic framework that has been used and extended in [19,20] and constitutes an essential ingredient of the model used in the present work. By assuming Maxwellianlike doubly-trapped electrons, these previous articles investigate the electron and ion macroscopic moments in the expansion, the asymptotic behavior far downstream, and in particular discuss the collisionless cooling of electrons observed in the expansion as a function of the plasma properties far upstream. Recently, Sánchez-Arriaga, Zhou, et al also introduced a 1D time-dependent Vlasov model to study the filling process of doubly-trapped regions due to the transient set-up of the plume in the MN and low collisionality [14,21]. Their findings show that these regions are not empty, and that the degree of filling depends on the effective electron collisionality. Their model however is unable of following the expansion to infinity, and a numerical plasma sheath forms at the end of the finite simulation domain, which may affect some of the behavior observed downstream.
With the exception of some full-PIC efforts applied to limited-domain magnetized expansions [15,22], the majority of 2D/3D existing models of MNs treat the electrons as a fluid species, mainly with an isothermal or polytropic closure relation. The 2D studies show a rich two-dimensional behavior of the plasma properties, electric currents, and magnetic forces that cannot be reproduced with simpler 1D models. Examples relevant to the present work are due to Ahedo and Merino in [2], where ions are modeled as cold and weakly-magnetized while electrons are hot and fully-magnetized, and in [23], where the fully-magnetized limit of electrons was approached. The main drawback of fluid electron models is the inability to compute the potential fall φ 0 − φ ∞ and the thermodynamic properties of the near-collisionless electrons self-consistently, for which a kinetic electron model is needed. An incorrect value of φ 0 − φ ∞ leads to erroneous estimation of the energyper-ion downstream, which affects the specific impulse and thrust computations.
This paper analyzes the steady-state cold-ion and hotelectron expansion in a 2D-axisymmetric divergent MN in the collisionless, full-magnetization limit. Ions are modeled as a simple fluid, while electrons are treated fully kinetically and constitute the main object of study. The goal is to investigate the 2D kinetic response of the electrons, and its consequences on the operation of the MN. In particular, the axial magnetic force density on the ions and electrons is determined and discussed under different globally-(but not locally-) current-free configurations. Along the way, we also revisit briefly the limitations of a global polytropic model, even when based on the same potential fall φ 0 − φ ∞ of the kinetic model. The problem is approached by combining two models previously presented, AKILES [24] and FUMAGNO [23], both of which have been open-sourced and are available to the community [25,26]. AKILES is a 1D steady-state kinetic plume code valid in the unmagnetized limit, which is adapted here to the fully-magnetized electron limit consistent with the formulation used in [13]. It is shown that the electron equations are mathematically analogous in the two limits due to the existence of alternative adiabatic invariants that play the same role in each case. FUMAGNO is a two-fluid code that describes the fully-magnetized expansion of a plasma in a MN. The code is extended here to accept the tabular kinetic results from AKILES, to yield the 2D, axisymmetric plasma response using the iterative approach proposed in [19]. This enables analyzing the kinetic cooling and anisotropization of electrons in the plasma domain, and extracting conclusions on the operation of the device.
The rest of this document is organized as follows. Section 2 summarizes the combined 2D fluid-kinetic MN model. Section 3 presents and discusses the simulation results, first for a single magnetic line and then for the whole MN. Finally, conclusions are gathered in section 4. Appendix A details the relation between the unmagnetized and fully-magnetized kinetic electron models. A first version of this research was presented in an international conference [27].

Fluid-kinetic plasma model
The plasma/MN model presented here consists of a fluid ion model coupled with a kinetic electron model. The selfconsistent electric potential map φ(z, r) and the plasma properties in the MN plume are found iteratively by combining both models. The applied magnetic field B is known and forms diverging, axisymmetric, nested flux tubes labeled by the magnetic streamfunction ψ, which satisfies Without loss of generality, B is assumed to point in the direction of the expansion. We define the right-handed orthonormal vector basis {1 , 1 ⊥ , 1 θ } with vectors parallel and perpendicular to B in the meridian plane and in the azimuthal direction, respectively; we also define the magnetic field angle with the MN axis, α = arctan(B r /B z ). The necessary plasma properties at the magnetic throat (section of maximum magnetic field and thus minimum magnetic tube cross-sectional area) are also known.

Fluid ion model
The fully-magnetized, collisionless plasma model of reference [23] is adapted here to solve only for the ion flow in the MN. A full derivation of the model equations can be found in that reference. In the following, e, m e , T e , m i are the electron charge, mass, and average temperature, and the ion mass, while L is the characteristic macroscopic length of the plume (e.g. the plasma beam radius, the gradient length, or the curvature radius of the magnetic field). We shall call s the meridian arc-length coordinate along each magnetic tube or line. A subindex '0' denotes the value of a variable on a given magnetic line at the magnetic throat, e.g. B 0 . These values therefore depend on the initial radius r 0 of that magnetic line at the throat. A double subindex '00' denotes the value at the origin of coordinates, i.e. at the center of the magnetic throat.
In the collisionless, steady-state, fully-magnetized limit, the ion gyrofrequency satisfies [23,28] eBL where M i = u i /c s is the ion Mach number, with c s T e /m i the sonic velocity. As long as this condition is met, ions will follow magnetic lines and inward ion detachment due to inertia will be small [29]. Observe that, as M i = O(1), this condition is roughly equivalent toΩ i 1.
Additionally, all perpendicular drifts are much smaller than the parallel ion velocity: u ⊥i u i ∼ c s . Hence, to first approximation, the ion fluid velocity is parallel to the applied magnetic field: u i = u i 1 + O(1/Ω i ). The cold, singly-charged ions follow the magnetic lines, along which their continuity and (parallel) momentum equations take the simple integral forms: where G i (ψ) and H i (ψ) are integration constants on each magnetic tube, and all other symbols are conventional. Observe that, in equation (3), 1/B plays the role of the cross sectional area of each diverging magnetic tube. Once boundary conditions n i = n i0 , u i = u i0 are given at the magnetic throat s = 0 of each magnetic tube, these expressions suffice to determine the ion properties n i and u i as a function of B(s)/B 0 , the local magnetic field strength normalized with its value upstream. Finally, while the zeroth-order ion velocity u i is parallel to B and j θi = en i u θi = O(1/Ω i ), the magnetic force density ( j θi B1 ⊥ ) may be not small and indeed it enters in the zerothorder momentum equation. This force term can be computed from the 1 ⊥ projection of the momentum equation of ions [23], where κ is the curvature of the magnetic lines defined by κ1 ⊥ = ∇ 1 , with ∇ ⊥ = 1 ⊥ · ∇ and ∇ = 1 · ∇, and the identity has been used. The applied magnetic field satisfies κ = ∇ ⊥ ln B everywhere except inside the magnetic coils.
Observe that ( j θi B) can be calculated a posteriori with equation (5), once the rest of the model has been solved.

Kinetic electron model
The reader is directed to references [13,24] for a full account of the 1D kinetic electron model in the magnetized and the unmagnetized limits. Below, only a summary of the main aspects of the model and its generalization to 2D are given for completeness. Appendix A presents the analogies between these two limits. Consider the steady-state, collisionless electron flow in a 2D, axisymmetric, slowly-diverging MN. In the limit of large electron gyrofrequencŷ electrons are well-magnetized and all perpendicular drifts are small. Note that this condition is automatically satisfied if (2) is true. Then, to first approximation, the electron macroscopic velocity is parallel to the magnetic field, u e = u e 1 + O(1/Ω e ), and, in the collisionless limit, electrons remain on their respective magnetic lines. On each magnetic tube, the distribution functions of downstream-marching and upstreammarching electrons, f + e and f − e , are assumed to be uniform in the gyrophase angle β and the azimuthal angle θ, and to depend only on the coordinate s along the magnetic tube, the mechanical energy E = m e (v 2 + v 2 ⊥ )/2 − eφ, and the magnetic moment μ = m e v 2 ⊥ /(2B), The mechanical energy E is a conserved quantity of electron motion, and μ is an adiabatic invariant whose gyro-average is conserved to second order [30] in 1/Ω e . Note that this approximate conservation fails if/where relation (7) is not met, introducing a dependency of f e on the gyrophase angle β. This may occur far downstream as B decreases, especially in weaker MNs; however, note also that in practical devices, the characteristic macroscopic gradient length L increases downstream, (partially) compensating the decrease of B, suggesting that electron demagnetization will only be an issue at great distances from the thruster exit in many cases. Nevertheless, the extent of the validity of relation (7) should be checked for each particular instance of MN. Finally, within the validity of this assumption, the kinetic equation for f e can be written in first approximation as This equation states that f + e and f − e are propagated as constants along s for each E, μ, as long as v = 0. The condition v = 0 determines the turning points for the electron trajectories. Expressed in terms of E, μ and the electric potential and the magnetic field strength along the magnetic tube, φ(s) and B(s), this condition reads Equation ( With regards to the condition at s = ∞, this study does not consider backstreaming electrons, and therefore, regions of type 3 in the list above will be empty regions It should be noted that, since f − e is not known a priori at the throat of the MN, the full electron distribution function f e there cannot be determined until the full problem is solved. This means, in particular, that the density and average temperature at the throat n e0 and T e0 , which are in general different from n * e and T * e , are not known beforehand but must be computed as part of the solution. Moreover, even if we prescribe f + e to be a semi-Maxwellian distribution, f e need not be Maxwellian nor isotropic, as high energies (associated to free electrons) have f − e = 0. Nevertheless, as the simulations of next section will confirm, these differences and deviations at the throat are small in all cases of interest.
While the prescribed boundary conditions determine unambiguously the EVDF in regions 1-3, the isolated regions of type 4 can accommodate an arbitrary distribution function in this steady-state model. These regions can become populated during the transient thruster ignition [14] or by infrequent but non-zero collisions [21], or even by instabilities due to the gaps existing in the EVDF. Indeed, the existence of a non-Maxwellian distribution function in some regions of the plume, caused by the existence of empty regions, could lead to instabilities-such as two-stream instabilities between the upstream-marching and downstream-marching populations-that could induce transport among the different regions of phase space. To date, however, two-stream instabilities have not been observed in time dependent MN models [14,21], but other types of instabilities due to azimuthal oscillations in the MN have been experimentally identified [31]. Collisionality, on the other hand, suggests that in these regions an essentially-thermalized population (with itself and with neighboring regions of phase space) will exist after sufficiently long times. In the present study, these regions are considered to be populated with the same distribution function as upstream, i.e. a Maxwellian population function of reference density n * e and temperature T * e .
Once f + e (s) and f − e (s) are known everywhere, the v i v j ⊥ integral moment M i j of the distribution function, and in particular n e , u e , can be computed via direct integration. This integration can be carried out in (v , v ⊥ ) space, or alternatively, in (E, μ) space after applying the necessary coordinate transformations: Observe that only free electrons contribute to odd moments in v . Finally, the relevant magnitudes for this study, are defined below: where all symbols are conventional. The electron pressure tensor takes the form P e = n e [T ⊥e U + (T e − T ⊥e )1 1 ] where U is the identity tensor. Note that q e , q ⊥e are the parallel heat fluxes of parallel and perpendicular thermal energy, respectively, while q e is the parallel heat flux of total thermal energy. We also define the moments and variables above for each electron subpopulation by extending the integrals only to their corresponding region of phase space. These will be denoted with an additional subindex f , r, t for free, reflected, and (doubly)trapped electrons. For example, n e, f = M 00, f is the free electron density, and u e, f = M 10, f /n e, f the free electron velocity. Lastly, we also define the local cooling exponents γ , γ ⊥ and γ as γ = d ln(n e T e ) d ln n e ; γ ⊥ = d ln(n e T ⊥e ) d ln n e ; γ = d ln(n e T e ) d ln n e .
As for ions, the first electron fluid equations (which are just integrals of the Vlasov equation) have a particularly simple form. In the case of the continuity equation, the expression is analogous to that of ions: with G e constant along magnetic tubes. After neglecting the small electron inertia term [32], the electron momentum equation reads: The 1 -component of this equation reduces to [19]: This relation states the balance between pressure forces, including the (explicit, macroscopic) magnetic mirror force term along the parallel direction, and electrostatic forces. Finally, while j θe = −en e u θe is O(1/Ω e ), the magnetic force density term ( j θe B1 ⊥ ) in equation (15) is not small, and indeed is comparable to the perpendicular pressure term. This magnetic force term can be calculated from the 1 ⊥ -component of equation (15): As with ions, this force can be computed after solving the rest of electron variables. Incidentally, it is noted that the product of B with other quantities can be computed in the 2D MN in a similar manner. Such is the case of the perpendicular and azimuthal fluxes of thermal energy. It can be easily checked that the product of B with the perpendicular flux of thermal energy is zero to first order in 1/Ω e [33].

Iterative solution method
The ion and electron models are coupled through the electrostatic potential φ, which is part of the solution. Moreover, the quasineutrality condition, which is imposed at each point of each magnetic line, provides Additionally, either the value of the total potential fall along each magnetic tube ψ, φ 0 − φ ∞ , or the value of the (parallel) net electric current density j = j i + j e = en e0 (u i0 − u e0 ) at some point, has to be fixed. From equations (3) and (14) it follows that j = j 0 B/B 0 . For concreteness, we shall fix φ ∞ = 0 as our potential reference, and define the following magnetic tube variablesΦ As it will be shown later, these two variables are mutually dependent, i.e.Φ 0 =Φ 0 (Ĵ 0 ) for given plasma properties upstream. If the MN is to be globally current-free, then the net current across any z = const. cross-section must be equal to zero, 2π j cos αrdr = 0.
Lastly, it is assumed that ions undergo a sonic transition at the magnetic throat. Inspection of the plasma momentum equation in the parallel direction (which results from combining the differential version of (4) with (16) to eliminate φ) shows that the relevant sound velocity to define the ion Mach number, M i = u i /c s , involves the value of T e and its local variation with n e at the throat (determined by γ ), c s = γ T e /m i . The sonic condition can be written as: It turns out that the spatial dependence of the model on s along each diverging magnetic line can be expressed in terms of B by inverting the relation B = B(s). Moreover, the ion and electron models along a single magnetic line can be normalized with the values at the magnetic throat B 0 , n i0 , n * e and T * e , variables which are naturally line-dependent. This conveniently enables solving the coupled model to obtain φ and the rest of plasma properties on each magnetic line as a function of B/B 0 . The normalized line-by-line model depends only on the parametersΦ 0 orĴ 0 , and on the mass ratio m i /m e . The latter is fixed at 491.689, which is the corresponding value for xenon. As a side note, and to facilitate comparison with previous works, it is noted that these values of M i0 and m i /m e determine the parameter χ = u i0 / T * e /m e used in reference [24], which is roughly 0.002 in all simulations shown in next section.
An iterative approach to find the solution of the problem is then established as follows. First, a guess potential function φ is used on the ion and electron models to compute all plasma properties on each considered magnetic line. Note that n e0 , u e0 , γ 0 and T e0 are not known a priori as implied above, and are obtained as part of the solution. This affects the computation of M i0 and j 0 . Consequently, equation (18), one of (20), and (22) must be evaluated to define the iteration error. Finally, the information on the error of these equations is used to update the function φ, and the process is repeated until the error is below a chosen tolerance.

Results and discussion
In this section, the fluid-kinetic model is first applied to a single magnetic line, and the integral moments of each electron subpopulation are discussed. Then, 2D MN expansions are studied to discuss the generation of magnetic thrust.

Evolution along a single magnetic line
The normalized solution of the electric potential φ − φ 0 along a generic magnetic line is shown in figure 2. Most of the potential fall is localized in the near region of the MN: around 75% of it has already occurred by B/B 0 = 0.01. By B/B 0 = 0.001 (i.e. an area expansion ratio of 10 3 ), 93% of the expansion has been completed. The electron density in figure 3 shows that n e /n e0 goes roughly as B/B 0 , except at the beginning of the expansion, where ion acceleration is higher and n e drops at a higher rate. Free electron density, n e, f , decreases at a slightly lower rate than the overall density n e , which means that free electrons gradually gain in importance downstream. The electron bulk velocity u e increases downstream; in fact, since the electron and ion properties satisfy u i = u e in caseĴ 0 = 0 this curve also depicts how the ions gain speed in the MN as the electrostatic potential decreases. As advanced before, most of the plasma acceleration takes place at the beginning of the expansion. Observe that the velocity of free electrons satisfies u e, f = n e u e /n e, f , and consequently has a value three orders of magnitude larger than u e . u e, f is non-monotonic, decreasing in most of the domain except initially, coinciding with the region of high ion acceleration. This also corresponds, roughly, with the region where reflected electrons are the dominant subpopulation.
The electron parallel and perpendicular temperature components T e and T ⊥e and the average temperature T e of figure 4 are nearly constant in the initial part of the expansion, in agreement with the near-isothermal behavior expected close to the throat. Further downstream, electrons begin to cool down gradually. This electron cooling becomes apparent roughly from B/B 0 = 0.1 onward. Far downstream T e goes to a small but non-zero asymptotic value T e = 0.012, while T ⊥e goes to zero [13]. Interestingly, the development of noticeable electron anisotropy occurs far downstream only [20]: while T ⊥e /T e 1 initially, by B/B 0 = 10 −3 , T ⊥e /T e 0.94, and by B/B 0 = 10 −4 , T ⊥e /T e 0.82.
The evolution of the temperature components of each subspecies is markedly different. These temperatures, weighted with the corresponding subspecies density, account for the behavior of the overall electron temperature discussed before. After an initial drop, free electrons have increasing T e, f along the expansion. This residual thermal dispersion is responsible for the existence of the asymptote in T e and consequently in T e , as every other subspecies temperature components vanish downstream. Temperatures T e,r and T e,t decay with the same rate downstream. Regarding the perpendicular temperatures, T ⊥e goes as T ⊥e,t , since doubly-trapped electrons give the largest contribution to this temperature downstream. However, T ⊥e, f and T ⊥e,r drop at a faster rate (similar for both), and soon become negligible. The average temperature of each subspecies follows the trend of the dominant component in each case, T e or T ⊥e . Interestingly, the anisotropy ratio is also quite different for each subspecies. While doubly-trapped electrons have T ⊥e,t /T e,t 1 downstream (consistent with the asymptotic analysis of [19]), free and reflected electrons have ratios going to zero at different rates.
More insight on electron cooling can be obtained plotting the electron temperature components of each subpopulation versus the density in logarithmic scale as in the right panels of figure 4. The local slope of the lines in these diagrams corresponds with the local values of γ , γ ⊥ , γ respectively. The non-monotonicity of γ , f of free electrons is evident, as is the near-adiabaticity γ ,t 5/3 of doubly-trapped electrons after the initial stages of the expansion. (The initial behavior of the double-trapped electron curves in this plots is a result of their density being zero at the MN throat.) The perpendicular γ ⊥ of reflected and doubly-trapped electrons have a similar behavior and remains slightly below 5/3 downstream, while γ ⊥ f is above this value. The local exponent γ of the global population is seen to go from near-isothermal in the near region to near-adiabatic downstream. This suggests that a two-value polytropic model may be more apt than a simpler global polytropic model to capture some of the complexity of the electron expansion. The transition starts roughly when n e /n e0 = 10 −2 -10 −3 , and the behavior downstream is indeed dominated by the thermal energy of doubly-trapped electrons.
The flow along 1 of parallel and perpendicular electron energy is m e M 30 /2 and m e M 12 /2, respectively. Following the definitions given in section 2.2, one can write the 'energy fluxes per electron' Both energy fluxes per electron are small in the whole nozzle, although the parallel one goes to an asymptotic value while the perpendicular one vanishes downstream. The major contribution to the parallel energy flux per electron is q e /n e , which is essentially constant downstream and about 4.42T e0 T e0 /m i as shown in figure 5. When the analysis is restricted to the free electrons only, though, the convective terms (i.e. the first two terms in the right-hand side of (23)), and especially m e u 3 e, f /2, dominate over q e, f . The perpendicular energy flux per electron receives comparable contributions from u e T ⊥e and q ⊥e /n e , both of which vanish at infinity. However, the perpendicular energy flux of free electrons is dominated by the convective term u e, f T ⊥e, f . Both q ⊥e /n e and q ⊥e, f /n e, f change sign in the first part of the expansion, the former going from positive to negative, while the latter from negative to positive. Lastly, the total energy flux (i.e. the sum of the two) and the total heat flux are dominated by their parallel parts. In the free electrons, q e, f /n e, f inherits the change of sign of its perpendicular part, q ⊥e, f /n e, f . The stark differences between the behavior of the moments of the full electron population and those of its subpopulations highlights the importance of carefully distinguishing among them. In the particular case of the heat fluxes discussed above, these differences arise due to the disparate density, velocity, temperature T e , perpendicular temperature T ⊥e , and average temperature T e . Lines represent free (blue circles), reflected (green squares), doubly-trapped (red diamonds) electrons, and the whole electron population (solid black line). In the plots on the right, the slope of each line is γ , γ ⊥ , γ, respectively; the slope for γ = 5/3 is provided as reference. In the bottom plots, the dashed lines indicate the solution of the global polytropic model with sameΦ 0 for comparison (γ g = 1.155).
and temperature values of the whole population and those of free electrons (see the definitions of q e , q ⊥e , q e in section 2.2).
The overall behavior of the moments shown agrees with those previously reported in [13,19,34] when integration is solved from a far-upstream plasma reservoir rather than from the throat, and interestingly, are similar to those of unmagnetized plasma expansions as well [24].
We now turn our attention to the role ofĴ 0 on the electron kinetics. A zero-current magnetic line is characterized bŷ Φ 0 = 7.45. As figure 2 illustrates, a positive or negative net currentĴ 0 in the line results in a stronger or weaker total electrostatic potential dropΦ 0 , respectively. The maximum value ofĴ 0 is 1, which corresponds to zero electron current and Φ 0 → ∞, j i = j i0 = en e0 c s0 . The initial part of the expansion has roughly the same slope of φ − φ 0 independently ofĴ 0 , but further downstream differences arise. AsĴ 0 increases-and with itΦ 0 increases too-the amount of free electrons in the magnetic line is reduced. This leads to the expected decrease of n e, f and increase of u e, f (which is the centroid, in velocity space, of the high-energy tail of the electron distribution associated to the free electrons). Simultaneously, the electrostatic confinement of electrons becomes larger (i.e. reflected and doubly-trapped electrons gain even more relevance). Obviously, as the total potential fallΦ 0 becomes stronger, the final ion velocity increases.
The main aspects of interest in the behavior of the electron temperatures asĴ 0 andΦ 0 increase are discussed next. The larger variations of partial temperatures (not shown) happens in the confined populations, while changes on the free electrons are minor. T e,r , T e,t , and consequently T e (and likewise T e,r , T e,t and T e ) increase in the first part of the expansion and remain near-constant farther downstream. The effect on T e is shown in the first panel of figure 6. T e and T e still drop to a final asymptote far downstream, and since the asymptotic value of these temperatures is determined by the free electrons only, as these become less relevant, this value decreases. The change in perpendicular temperature of the free and reflected electrons is negligible, but the doubly-trapped T ⊥e,t increases with increasingĴ 0 , and this drives a similar behavior in the global T ⊥e . Regardless, this temperature component vanishes downstream at the same rate as before. As a result, the plasma anisotropy ratio T ⊥e /T e , shown in figure 6, is very dependent onĴ 0 : the higherĴ 0 (i.e. the strongerΦ 0 ), the longer the region where the plasma can be considered isotropic. This is consistent with the smaller fraction of free electrons in the population, and correspondingly, the smaller 'gap' that they produce in f − e . The behavior described so far suggests two limits in which the electron temperatures follow a particularly simple pattern: first, and of higher practical interest, asΦ 0 → ∞ and all electrons become confined, electron temperature becomes constant and isotropic everywhere; incidentally, reasoning along these lines is what justifies the isothermal approximation in some near-region MN models [2]. Second, asΦ 0 → 0 and all electrons become free, electrons reach a larger asymptotic value of T e and they do so earlier in the expansion, while T ⊥e follows a simple decay rule. Note that, in this second case, the electron anisotropy is large in most of the expansion.
The possibility of defining an 'effective', global polytropic electron model where T e ∝ n γ g −1 e in the MN has already discussed elsewhere [12,17,20], pointing out its limitations. To conclude this section, we make an additional comment regarding this topic in the light of the present results for each subpopulation. The polytropic exponent γ g that yields the same total potential fallΦ 0 along the magnetic line is We find γ g = 1.155 for theĴ 0 = 0 case, a value well in line with those commonly derived from experimental measurements in MN-based thrusters [16,17,35]. It is noted, however, that experimental calculations of γ g are quite device-and setup-dependent: for instance, near-adiabatic values have been found (assuming two-and three-degrees of freedom electrons) in other works [36,37]. These differences hint at the potential effects that vacuum chamber walls, background pressure, and measurement techniques can have on the experimental determination of the kinetic features of the electron population. Incidentally, it is also noted that most of the existing experimental data show cooling profiles with a non-constant local exponent γ, as in the present kinetic model, highlighting one of the drawbacks of using a global γ g . Figure 2 shows the electrostatic potential for this equivalent polytropic model as a dashed line for comparison with the kinetic solution. While early in the expansion the solutions approximately coincide (up to B/B 0 = 0.3 roughly), it is clear that the polytropic model approaches the downstream asymptote at a slower rate. Consequently, the polytropic model tends to underestimate the rate of ion acceleration in the MN, as evidenced by the profile of ion velocity u i in figure 3. This has an effect on the determination of the magnetic thrust produced up to a z = const. section, as will be discussed in section 3.2. The plasma density n e , on the other hand, is approximately the same in the polytropic and kinetic models, as this variable is mainly dictated by the area expansion ∝ B 0 /B. The global polytropic model also misses the rich behavior of electron temperature T e revealed by the kinetic model, as shown in figure 4. It neglects the development of anisotropy, and does not reproduce the quasi-isothermal behavior in the near region nor the near-adiabatic behavior farther downstream. Incidentally, however, the rate of cooling of reflected electrons, γ r , is found to coincide with γ g in most of the expansion. Nevertheless, as reflected electrons are a minor subpopulation except near the throat, is not indicative of the overall behavior of the electron species.
Thus, invoking a (global) polytropic approximation, even if chosen consistently with the final asymptotic behavior of the plasma plume, might have an important effect on the calculation of the MN performance figures like thrust, specific impulse, divergence and conversion efficiencies, etc. in practice, when only a finite length of nozzle is investigated.

Two-dimensional plasma expansions
In the fully-magnetized limit, the solution to the 2D MN can be constructed line by line using the 1D solution discussed above. Therefore, the basic behavior of the 2D expansion follows directly from the 1D electron kinetic model. However, there are two aspects of MN physics that cannot be analyzed in a 1D geometry, which are the focus of this section: the influence of the initial radial profiles at the exit of the plasma source on the electric currents and electron properties, and the generation of magnetic thrust. Central to this discussion is the fact that each individual magnetic line needs not be current-free, although the integral electric current through the MN must satisfy the global current-free condition of equation (21) for continuous operation in space. Figure 7 depicts the initial radial profile models chosen for this study, which intend to be representative of various MNs of interest without being a exact reproduction of any particular thruster implementation. In all simulation cases shown below, the plasma source exit has a radius R 0 , and the MN is generated by a single magnetic coil of radius R L = 10R 0 located on the z = 0 plane. Initially, ions are sonic (u i = c s ) and the plasma density profile is similar to that of [23], peaking on the centerline, and dropping 3 orders of magnitude radially; n e (0, r) = n e0 (r) = n e00 10 −3r 2 /R 2 0 (26) (As a reminder, the double 00 notation in a subindex denotes the value at the origin z = r = 0). The value n e00 is used to nor- Figure 7. Initial profiles of plasma density n e , electrostatic potential φ, and electron average temperature T e in the 2D simulations F, PHID, PHII, TD, and TI at z = 0. In all cases, the potential far downstream is set to zero φ ∞ = 0 on all magnetic lines. Quantities n e00 and T e00 are the density and electron temperature values at the origin and are used to normalize the graphs. malize all densities hereafter. For consistency, the electrostatic potential far downstream is the same for all magnetic lines and is fixed at φ ∞ = 0, meaning that φ 0 > 0 initially. The first three simulations have different profiles of φ 0 : simulation F corresponds to a flat initial φ 0 profile with zero charge current on each magnetic line, and is used as the reference in the comparison with other cases; simulation PHID has a radially-decreasing potential profile, and is expected in HPTs and ECRTs with large wall presheaths (e.g. if the magnetic confinement in the source is not large enough to keep a nearflat φ 0 profile); lastly, simulation PHII has a radially-increasing potential profile, which may be regarded as an approximation of AFMPDTs where the axial cathode drives down the potential locally: φ(0, r) = φ 00 + T e00 r 2 /e (case PHII).
These three simulations feature a constant radial temperature profile at the magnetic throat, T e (0, r) = T e00 , the value at the origin used for energy normalization.
The former is relevant for coaxial ECRTs, which exhibit a peak of electron temperature near the axis [38,39], while the latter may be a better model for some HPTs whose measured electron temperature rises toward the lateral walls [40]. In all five simulations, the value of eφ 00 /T e00 in each case is computed to satisfy the global current-free condition (21), i.e. j cos αrdr = 0: 7.45 for F, 7.59 for PHID, 7.31 for PHII, 7.20 for TD, and 7.74 for TI. Figure 8 exhibits the solution in the reference case F, while figure 9 displays the electrostatic potential φ, slip velocity u i − u e , and average electron temperature T e in simulations PHID, PHII, TD, and TI. The following observations can be made. Firstly, a non-uniform initial profile of either φ or T e affects the electrostatic potential only in the near region, as farther downstream φ → 0 along all lines regardless of the initial profile. The ion velocity in the meridian plane u i changes accordingly to the potential fall along each line. Secondly, these initial profiles determine the electric current map in the MN, which-while globally current-free-is in general not locally current-free. A radially-decreasing φ or a radiallyincreasing T e initial profile leads to a larger positive current at the axis and a higher electron flux near the periphery. Thirdly, it is noted that T e is also affected by either initial profile, and not only its own: a larger initial φ along a magnetic line (i.e. a larger total potential fall) extends the distance from the throat in which the electrons remain near-isothermal before cooling down, as discussed in the previous section, and as can be seen comparing simulations F, PHID and PHII. While the electrons start as an isotropic species in the cases under analysis, the anisotropy ratio T ⊥e /T e downstream (not shown) also depends on the initial profile: regions with magnetic lines on which the total potential fallΦ 0 is smaller have a larger fraction of free electrons, and consequently, a stronger anisotropy that manifests itself earlier in the expansion. In comparison, the density map n e is essentially dictated by the initial profile n e (0, r) and the infinitesimal area of each flux tube, B 0 /B, since the role of acceleration in the continuity equation is small and limited to the near region. Therefore, n e remains essentially unchanged with respect to the reference simulation F.
The last four panels in figure 8 display the force terms (per particle) that appear in the perpendicular momentum equations of ions and electrons, (5) and (17), and which give rise to the magnetic force densities ( j θi B) and ( j θe B). Observe that these two force densities can be combined to yield the total magnetic force density ( j θ B) = ( j θi B) + ( j θe B) in the MN: where the quasineutrality condition has been used to simplify the expression. This force density is responsible for the radial plasma confinement, which controls the shape of the plume expansion, and for the (axial) magnetic thrust generation [2]. Indeed, the magnetic thrust force is the reaction, felt on the permanent magnets or coils that generate the MN, to the integrated axial force in the plasma volume: Note that as α is small (B r B z ) in the near region, the axial magnetic force is always smaller than the radial one. A positive (perpendicular) magnetic force density as defined in (32) is paramagnetic and leads to a negative contribution to magnetic thrust, while a negative magnetic force density is diamagnetic and leads to positive magnetic thrust. A net diamagnetic ( j θ B) is required for the MN to be meaningful as an acceleration device. Note that the total thrust of a given thruster also contains a so-called pressure or source contribution, equal to the momentum flux present at the MN throat [2].
Inspecting figure 8, it is clear that the perpendicular electron pressure ∇ ⊥ (n e T ⊥e ) is the dominant contribution to F M , and that it is diamagnetic. This is consistent with the fact that pressure is the driving force of the expansion in a MN-based thruster, which from this viewpoint can be considered as an electrothermal device. On the other hand, the ion curvature term κn e u 2 i and the macroscopic mirror term κ(T e − T ⊥e ) are paramagnetic contributions, but they are always much smaller than the pressure term. Furthermore, the relaxation of the fullmagnetization assumption in the model is known to lower the magnitude of the ion curvature term, as ions gradually detach from the magnetic lines and their trajectories become straight [29], resulting in an even smaller paramagnetic contribution. This makes manifest that excessive ion magnetization can have a detrimental effect in magnetic thrust generation (and plume divergence).
Evidently, the contribution of the electrostatic field term to the magnetic thrust force is zero, as stated by equation (32). Nevertheless, the characterization of this field can shed additional light on the role that ions and electrons play in magnetic thrust generation, as it has an opposite effect on ( j θi B r ) and ( j θe B r ). The perpendicular electric field is strongly dependent on the initial profiles of φ and T e , and therefore is expected to vary among thruster types. In the fully-magnetized limit, in the flat-profile simulation F of figure 8, this field is small everywhere and points inward. This leads to equal ion paramagnetic and electron diamagnetic contributions to magnetic thrust that cancel out when summed together, and nonetheless negligible compared to the contribution of electron pressure. Figure 10 displays this field in the plume of the remaining simulations. Cases PHID and PHII, which start with a non-uniform φ profile at the throat, feature accordingly a perpendicular electric field that is three orders of magnitude larger than F's, This field is outward and inward pointing, respectively in PHID and PHII. In both simulations, the field peaks at the throat where the non-uniform profile is imposed, and decays downstream. Simulations TD and TI, with non-uniform initial T e profiles, also display a far stronger perpendicular electric field than the flat case F, albeit smaller than in the previous two cases; the field points inward in TD, and outward in TI. In these cases, the perpendicular electrostatic field arises due to the differential Figure 8. Reference simulation F. The first four panels present the electrostatic potential φ (with φ ∞ = 0 as potential reference), plasma density n e , ion velocity in the meridian plane u i , and average electron temperature T e . The last four panels display the terms in the perpendicular ion and electron momentum equation that give rise to the magnetic force. This simulation corresponds line to line with the 1D normalized solution with zero charge current of section 3.1.
expansion between magnetic lines, as dictated by the different values of the electron temperature upstream. Figure 11 presents the magnetic thrust density due to ions and electrons in the latter four simulation cases. In agreement with existing fluid model results [2], the contributions reach their maximum in a small domain at mid-radius in the near region throat, where the product n e T e B r is large. It is clear that electrons bear the major part of the magnetic force in the cases under analysis, while the ion contribution is always secondary. Indeed, magnetic thrust due to ions is only relevant when perpendicular electric fields are large, as this is the dominant term in equation (5): whenever this field points radially outward, this term is diamagnetic and contributes positively to magnetic thrust (PHID, TI); conversely, it is paramagnetic when the field points inward (PHII, TD). It is interesting to note that the socalled ion swirl thrust, associated to the azimuthal current of ions, is actually negative when the initial φ profile is minimum at the axis (case PHII), which is the expected situation in an AFMPDT with cold ions. Likewise, a plasma jet with hotter electrons in its core, as observed in some ECRTs (represented here by simulation TD), would lead to slightly paramagnetic ions.
It is remarked that, after cancelling out the opposing electric field terms in ions and electrons, the different simulation cases portrayed here show a roughly equal integrated magnetic thrust force F M after proper normalization. In all studied cases, roughly 90% of the maximum magnetic thrust is achieved already when the radius of the plasma beam is about 10R 0 . Observe also that MNs with non-negligible ion pressure may Figure 9. Electrostatic potential φ (with φ ∞ = 0 as potential reference), current j/(en e T e00 /m i ) = (u i − u e )/ T e00 /m i , and average electron temperature T e in simulations PHID, PHII, TD, and TI. Note that the scale of these plots differs from that of figure 8. behave differently, as it causes a diamagnetic contribution like that of electrons, and then it would be the sum of electron and ion pressures that would drive the expansion in the MN [12]. Ion pressure is expected to be important in several devices, for example the VASIMR.
Finally, to highlight some of the differences of the kinetic electron model with global polytropic models, figure 12 displays the integrated magnetic thrust generated up to an axial position z for case F and its polytropic approximation with the same asymptotic potential fallΦ 0 . As a consequence of the slower potential drop that was illustrated in figure 2, the polytropic model predicts a lower growth rate of F M . This may have practical consequences when measuring thrust with a MN in the laboratory, where the finite dimensions of the vacuum chamber will limit the expansion.
It is worth commenting on the validity of assumption Ω e 1 in equation (7), on which the conservation of μ of the electron model hinges. While B decreases downstream, this effect is partially counteracted as the characteristic gradient length L increases, and T e cools down. It has been checked that, for the reasonable values of B 00 = 2000 G and T e00 = 10eV, Ω e > 10 5 in the full domain of figure 8. To perform this computation, the smaller between the radius of curvature 1/κ and 1/|∇ ln B| was used as the (local) characteristic length L. Thus we conclude that this assumption of the electron model is applicable until very far downstream, in practical MNs. Additionally, it is noted that the ion model fails if/where conditionΩ i /M i 1 of equation (2) (which is a simplified version of the detachment condition in reference [29]) is not satisfied. This condition is indeed more restrictive than the corresponding electrons' due to the factor m i /m e , which is 500 for xenon propellant, and due to ion acceleration, which makes M i reach hypersonic values around 5-10 downstream. Failure to meet this condition will occur first downstream in the plume periphery, where inward ion detachment will become important first [10,29], leading to a modification of the ion streamline geometry and the solution to the perpendicular electric fields. For the numerical example of the previous paragraph, condition (2) is well-satisfied in the full simulation domain, except perhaps at the north-east corner of the plots in figure 8, whereΩ i /M i drops to 7. At any rate, ion detachment is not expected to modify substantially the kinetic aspects of the electron solution discussed here.

Conclusion
The open-source kinetic AKILES and fluid FUMAGNO models have been coupled together to provide a fluid/kinetic description of the 2D plasma expansion in an axisymmetric MN in the fully-magnetized, collisionless limit. The resulting model can be solved line by line, simplifying the iterative solution process.
Electrons, which can be divided into free, reflected, and doubly-trapped populations, do not follow any simple closure relation, develop a mild anisotropy that grows downstream, and exhibit non-trivial cooling rates and kinetic heat fluxes that are missed by simpler fluid models. Doubly-trapped electrons are seen to become the dominant population downstream. The net charge current in the magnetic line is determined solely by the fraction of free electrons, and is directly related with the total potential fall along it. The fraction of free electrons dictates the anisotropy ratio in the expansion and the extent of the region in which electrons can be approximated as isothermal.
Our two-dimensional study has covered simulation cases which, while globally current-free, feature differential charge currents among magnetic lines, as dictated by different initial profiles at the magnetic throat, relevant for different types of MN-based thrusters. The perpendicular electrostatic field, electron temperature, and electron anisotropy are affected. While perpendicular electron pressure is the origin of magnetic thrust in present simulations and the dominant contribution in all cases, ion streamline curvature and the electron mirror term produce small paramagnetic (i.e. drag) contributions. The perpendicular electric field in the plume varies widely with the initial profile, and determines the share of magnetic thrust generation of ions and electrons. A radially-decreasing electrostatic potential or a radially-increasing electron temperature profile at the throat result in diamagnetic, thrust-producing ions, while other choices can lead to a paramagnetic ion contribution. This discussion is also relevant in the context of vacuum-chamber experiments of MNs, where the electrical boundary conditions at the thruster and the chamber walls determine the charge current on each magnetic line. In fact, it should be possible to drive differential currents in a MN experiment using separate downstream electrodes at different potentials immersed in the plasma plume, a test that would help validate the present model. Note that in a practical (flight) device, the current on each magnetic tube cannot be controlled in this manner.
The main takeaways of the study are: (i) the richness and complexity of the kinetic response of electrons in 2D collisionless, magnetized expansions, which are missed by approximate fluid models; (ii) the link between the total potential fallΦ 0  assumption would require a different iterative approach to solve the model, as now the plasma properties on all magnetic lines must then be solved for simultaneously.
We observe that, due to this approximation, the present model is not adequate for the study of plasma detachment in the downstream region, a problem that has received substantial attention in the literature [10,29,[41][42][43], and in which demagnetization plays an important role. While this topic has been outside of the scope of this work, we quickly note that the detachment problem can be broken down into the detachment of ions and the detachment of electrons. As shown in [29], the former occurs early in the expansion of a practical MN, as the inertia of the accelerating ions wins over the magnetic and electric forces that try to deflect their trajectories to match the magnetic lines. The detachment of the massand momentum-carrying ions is deemed essential for the correct operation of the MN. Electron detachment, on the other hand, is expected to occur much farther downstream due to the higher electron magnetization degree, when finite Larmor radius effects become predominant [43]. The modeling of the demagnetization and detachment of electrons, which can be described as a warm, charge-neutralizing cloud with a global response, remains a still-open problem.
Future lines of research should also study the effect of a background plasma, whose electrons would backstream into the MN and affect the expansion. This study is straightforward with the current model. Additionally, it is also possible to adapt the 2D AKILES/FUMAGNO model to 3D MNs and study directional thrust, relevant for thrust-vector-control MNs [28].