Plasma instability of magnetically enhanced vacuum arc thruster

A two-fluid flowing plasma model is applied to describe the plasma rotation and resulted instability evolution in magnetically enhanced vacuum arc thruster (MEVAT). Typical experimental parameters are employed, including plasma density, equilibrium magnetic field, ion and electron temperatures, cathode materials, axial streaming velocity, and azimuthal rotation frequency. It is found that the growth rate of plasma instability increases with growing rotation frequency and field strength, and with descending electron temperature and atomic weight, for which the underlying physics are explained. The radial structure of density fluctuation is compared with that of equilibrium density gradient, and the radial locations of their peak magnitudes are very close, showing an evidence of resistive drift mode driven by density gradient. Temporal evolution of perturbed mass flow in the cross section of plasma column is also presented, which behaves in form of clockwise rotation (direction of electron diamagnetic drift) at edge and anti-clockwise rotation (direction of ion diamagnetic drift) in the core, separated by a mode transition layer from $n=0$ to $n=1$. This work, to our best knowledge, is the first treatment of plasma instability caused by rotation and axial flow in MEVAT, and is also of great practical interest for other electric thrusters where rotating plasma is concerned for long-time stable operation and propulsion efficiency optimization.


I. INTRODUCTION
Plasma propulsion is generally caused by jetting plasma in the opposite direction, according to Newton's third law, and accelerated either by electric force or magnetic force or both of them. [1][2][3][4][5][6] External magnetic field is usually employed for efficient plasma generation and propulsion enhancement and control. [7][8][9][10][11][12] Due to the nonuniform configurations of equilibrium magnetic field, discharge area and plume, plasma rotation driven by Lorentz force commonly occurs in various electric thrusters. [13][14][15][16] A few analytical models have been developed or/and applied to describe this flowing phenomenon, 17-21 but little attention was given to the resulted plasma instability which, however, can effect the propulsion efficiency, precise control, durable reliability and life time significantly. 22 arc thruster (MEVAT), [24][25][26][27][28] as an example and studies the instability evolution caused by plasma rotation and axial flow in detail. Schematics of the typical coaxial-type and ring-type MEVATs are shown in Fig. 1, which illustrate the radial expansion of plasma plume across the confining field lines, leading to azimuthal rotation. Axial equilibrium field is usually employed to reduce the divergence of plasma plume and thereby increasing the propulsion efficiency, which is relatively low for unmagnetized vacuum arc thrusters, making the MEVAT a promising candidate to provide micro propulsion for small spacecrafts. 28 The findings achieved here are also applicable to other types of electric thrusters with the involvement of external magnetic field, as long as the plasma rotation and flow are concerned for stable, efficient and safe operation.
Specifically, a two-fluid flowing plasma model developed originally for interpreting wave oscillations in vacuum arc centrifuge, 29 which is a cylindrical, rapidly rotating, low temperature, and confined plasma column, 18 will be employed, together with a shooting numerical scheme for finding solutions. It will show that the instability strength peaks near the maximum of equilibrium density gradient, an evidence of resistive drift mode, and it becomes larger for higher rotation frequency, higher field strength, lower electron temperature and lower atomic weight. Temporal evolutions of perturbed mass flow in cross section will be also shown. The paper is organized as follows: Sec. II describes the theoretical model and steady-state and perturbed solutions, Sec. III presents the employed numerical scheme, dispersion relation, fluctuation structure and parameter dependence, and Sec. IV summarizes the whole paper and remarks the possible applications of this work to other electric thrusters.

A. Model assumptions
The employed two-fluid flowing plasma model is based on the following assumptions: 1. Ions of different charge can be treated as a single species with average charge Z.
2. Plasma is quasi-neutral so that n e = Zn i .
3. Steady-state plasma is azimuthally symmetric and has no axial structure.
4. The effects of plasma fluctuation on external magnetic field is negligible. 5. Finite Larmor radius and viscosity effects are not considered. 6. Electron inertia can be neglected for the range of frequency considered. 7. Ion and electron temperatures (T i and T e ) are uniform across the plasma column.
8. Steady-state ion density distribution is in form of n 0 = n i0 exp[−(r/R) 2 ], with n i0 the on-axis ion density and R the characteristic radius where the density is 1/e of its on-axis value.
9. Steady-state velocities of ions and electrons are in forms of v i = (0, ω i r, v iz ) and v e = [0, ω e (r)r, v ez (r)], respectively, with ω i the ion rigid rotor rotation frequency, v iz the ion uniform axial streaming velocity, ω e (r) the electron rotation frequency, and v ez (r) the electron streaming velocity.
10. Radial diffusion of both ions and electrons caused by electron-ion collision can be neglected. Moreover, length and time are normalized to R and 1/ω ic respectively, with ω ic = ZeB z /m i the ion cyclotron frequency, so that a normalized cylindrical coordinate system becomes (x, θ, ς) = (r/R, θ, z/R) and τ = ω ic t, with x and ς the normalized radial and axial positions, respectively.

B. Governing equations
The model consists of the momentum and continuity equations for ion and electron fluids: with terms defined as: The subscript i and e refer to ion and electron parameters respectively, ϕ represents the normalized radial velocity divided by x, Ω stands for the normalized rotation frequency, u ς shows the normalized axial velocity, λ labels the ratio between ion and electron temperatures, ψ behaves a convenient constant which for λ = 1 becomes the square of the normalized ion thermal velocity, χ gives a normalized electric potential φ, l i expresses the logarithm of the ratio of ion density n i to its on-axis value n i0 , and n s means the ratio of electron-ion collision frequency ν ei to its on-axis value ν ei0 . Moreover, δ phrases the normalized resistivity parallel to external magnetic field, with η L the electrical resistivity of Lorentz gas and γ E the ratio of the conductivity of charge state Z to that in Lorentz gas. 30
Here, χ c is an arbitrary reference potential. Because the axial current in this model is unconstrained, which is consistent with MEVAT boundary conditions, it can be arbitrarily set to zero (u iς0 = u eς0 ).

D. Perturbed solution
To solve the perturbed solution, a linear perturbation treatment with plasma parameters ζ is considered, namely Here, ε is the perturbation parameter, m is the azimuthal mode number, k ς is the axial wave number and ω is the angular frequency. To first order of ε, Eqs. (1)-(4) can be reduced to Here, ̟ = ω − mΩ i0 − k ς u ς0 is the frequency in the frame of ion fluid, C = 1 + 2Ω i0 , Ψ = (λ + Z)ψ, y = x 2 and a new dependent variable X 1 (y) has been introduced For large axial wavelength modes of the resistive plasma column, i. e. k 2 ς ≤ δ, Eq. (8) turns out to be a second order differential equation where and f (y) = F 2 e y with the normalized axial wave number F = k ς / √ δ. For odd m modes, the boundary conditions are g 1 (0) = 0 and g 1 (Y ) = 0 with the infinite radius Y representing the edge of plasma column. For even m, these conditions become g ′ 1 (0) = 0 and g 1 (Y ) = 0. We only consider unstable solutions for which ̟ i > 0. The solutions ̟ = ±C are stable and thereby discarded.

A. Numerical scheme and conditions
To solve Eq. (11) numerically for perturbed solutions, we make use of a shooting method, as did by Hole et al. 29 For m = 1 mode, the boundary conditions are g 1 (0) = g 1 (Y ) = 0. The gradient at edge g ′ 1 (Y ) is arbitrary because the differential equation is homogeneous. As a result, we set g ′ 1 (Y ) = 1. For given F , a trial ̟ is first chosen and then the solution is matched from edge to core. We adjust the complex frequency ̟ until the on-axis boundary condition is satisfied. The procedure starts from F = 0, for which an analytical solution for ̟ can be found from with n the number of radial nodes in the plasma column. This numerical scheme has been benchmarked by previous studies. 18,29 For the computational parameters, we refer to existing devices and choose published experimental data. Table I shows the typical parameters of MEVAT employed for the present computation, together with those of a plasma centrifuge (PCEN 29,31 ) for comparison. Although MEVAT has much higher electron and ion temperatures and lower rotation rate, the two-fluid model originally developed for PCEN can still well describe the flowing plasma in MEVAT. The plasma density of n i0 = 5 × 10 19 m −3 is close to the measurement by Keidar et al. 32 Trial electron and ion temperatures are 100 eV and 50 eV, respectively, for current in order of 1 kA 33 and with reference to a previous experiment. 32 The strength of external magnetic field B z = 0.1 T is typical for various MEVATs. 21,28 The axial stream velocity of ion fluid is set to 1.3×10 4 m s −1 according to the experiment using titanium (Ti) cathode, 34 which lies inside the typical range of ion velocity 10 − 50 km s −1 , 35-39 and meanwhile the ion species of Ti is chosen. The rotation frequency is assumed to be half the ion cyclotron frequency, which is close to that of plasma centrifuge 29,31 to see the effect of rotation on instability evolution more clearly, and rigid rotation is a reasonable assumption according to various jetting plasma devices. 18,31,40 Moreover, charge number of Z = 2 is considered throughout the paper as a common value for different cathode materials chosen here (see Sec. III D), and characteristic radius of R = 0.03 m is made use of referring to previous studies. 20,28 The radial density profile is shown in Fig. 2, a typical plasma distribution near the exit of MEVAT. 27, 28 We shall compute the dispersion relation and instability physics based on these experimental data in the following sections.

B. Dispersion relation
Before presenting the obtained dispersion curve, we should check the eigenfunction associated with the computation. A typical radial variation of solved eigenfunction for F = 0 is shown in Fig. 3, which clearly satisfies the boundary conditions of g 1 (0) = g 1 (Y ) = 0 for m = 1 mode. Moreover, it shows a radial mode of n = 1 with visible node located near y = 2. The computed dispersion curves in the range of F = 0 − 1 are illustrated in Fig. 4, compared with those for PCEN. It can be seen that the normalized growth rate of MEVAT peaks on axis and decreases monotonously with F , whereas the normalized growth rate of PCEN peaks off axis (F = 0.3) and exhibits a mode crossing near F = 0.55, which is caused by centrifugal instability. Relatively, the instability level (proportional to growth rate) is much more lower on MEVAT than that on PCEN for F > 0.3, due to significantly smaller resistivity along the external magnetic field (it will show later that the instability is a resistive drift mode). Further, the peak growth rate of ̟ i max = 0.324 corresponds to ̟ r = −0.205 for MEVAT, implying that the instability frequency is smaller than the sum of plasma rotation frequency and axial velocity, because ̟ = ω − mΩ i0 − k ς u ς0 is the frequency in the frame of ion fluid, and it propagates in the −θ direction (same to the direction of ion diamagnetic drift); whereas, the peak growth rate of ̟ i max = 0.39 corresponds to ̟ r = −0.007 for PCEN, stating that the instability is near stationary in the frame of ion fluid. Overall, the normalized frequency is lower on MEVAT for F > 0.4, which may be attributed to lower rotation frequency and ion temperature as revealed in a previous study. 18

C. Fluctuation structure
To show a cross sectional view of the mode structure, the vector field of linearly perturbed mass flow was calculated through m i (n i1 u i0 + n i0 u i1 ). We computed the perturbed velocity components u i1 = (xϕ i1 , xΩ i1 , u iς1 ),u e1 = (xϕ e1 , xΩ e1 , u eς1 ) and perturbed density  n i1 from the solution of g 1 (y) and following equations: 18 l i1 (y) = −g 1 (y)  Figure 5 displays the radial profiles of perturbed density n i1 (r) and equilibrium density gradient |n ′ i0 (r)|. We can see that the perturbed density peaks (6.89 × 10 15 m −3 ) at r = 1.68 cm, which is very close to the radial location r = 2.12 cm of maximum density gradient (1.43 × 10 21 m −3 ) in equilibrium state. This suggests that the observed instability shown in Fig. 4 may be a resistive drift mode, which is driven by plasma density gradient. Temporal evolution of the cross sectional mass flow is given in Fig. 6, for a period of t = 1 s. The time dependence was achieved by multiplying exp[i(2πt)] with m i (n i1 u i0 + n i0 u i1 ). Here, the external magnetic field and z point into the page, and coordinates x and y label the cross section of plasma column, namely x = r cos θ and y = r sin θ. We can see that there exists a radial layer inside which the mass flow rotates in the anti-clockwise direction (same to the direction of ion diamagnetic drift), while it rotates in the clockwise direction (same to the direction of electron diamagnetic drift) outside, indicating a circularly sheared flow near the layer. The radial location of this layer is around the mode transition radius of r = 4.24 cm shown in Fig. 5, from n = 0 mode to n = 1 mode. This new pattern of rotation is different from the symmetric rotation observed before, 18 , and may be attributed to the large difference between electron and ion temperatures and very low normalized resistivity along the equilibrium magnetic field. Moreover, the mass flow is largest around the peak density gradient at r = 2.12 cm, and drops to zero when approaching to the core and edge of plasma column, consistent with the boundary conditions of perturbed density (Fig. 5). Although only linear oscillatory response is considered here, the rotation may be damped by a similar nonlinear flow pattern as claimed before. 18

D. Parameter dependence
To guide the experimental design of an efficient MEVAT, this section is devoted to studying the effects of rotation frequency, field strength, electron temperature and cathode material on the growth rate and frequency of plasma instability. As shown in Fig. 7, the normalized growth rate and frequency decrease when the rotation frequency is reduced (till Ω = 0.2), consistent with a previous observation, 18 because the centrifugal force is descending. However, when the rotation frequency drops below Ω = 0.1, the instability becomes much stronger and the profile of dispersion curve changes remarkably, indicating the entrance into a different fluctuation mode. Figure 8 shows the influence of external magnetic field strength on the dispersion relation of MEVAT. Different from previous conclusion, however, the normalized growth rate increases with growing field strength, although the normalized frequency decreases as before. 18 Moreover, as shown in Fig. 9, the normalized growth rate also increases with descending electron temperature, different from previous observation, whereas the frequency drops slightly at expected. 18 We claim that this unusual variation of growth rate with field strength and electron temperature may be attributed to the extremely low normalized resistivity parallel to external field, namely δ = 7.8 × 10 −5 , which makes the plasma nearly an ideal MHD medium. Additionally, we varied the ion temperature from T i = 50 eV to T i = 10 eV, and found that the change in dispersion curves is even smaller and negligible. Finally, we compare the plasma instability for different cathode materials, which are commonly used in MEVAT, including Ti for previous sections, magnesium (Mg), cuprum (Cu), molybdenum (Mo) and wolfram (W). As the atomic weight is increased, the normalized growth rate drops and the normalized frequency increases slightly, which is perhaps comprehensible in the sense that ions become heavier so that move slower.

IV. CONCLUSIONS
To describe the plasma instability associated with plasma rotation and axial streaming, a two-fluid flowing plasma model developed originally for vacuum arc centrifuge is applied to MEVAT. Based on typical experimental data, the dispersion curve showing the growth rate and frequency of instability evolution is first computed via a numerical shooting method. The strongest instability occurs on axis for rotation frequency close to that of vacuum arc centrifuge, and its frequency is smaller than the sum of rotation frequency and axial velocity. Then the perturbed density is calculated through a linearized technique, which shows a radial mode transition from n = 0 to n = 1 at r = 4.24 m, and it peaks near the radial location of maximum equilibrium density gradient, suggesting that the observed instability is a resistive drift mode driven by density gradient. Moreover, the temporal visualization of perturbed mass flow in the cross section of plasma column shows an anti-clockwise rotation (same to the direction of ion diamagnetic drift) inside the mode transition layer, and clockwise rotation (same to the direction of electron diamagnetic drift) outside. This implies that the plasma instability is dominated by different particle species in the core (ion) and at edge (electron). Finally, parameter scan shows that the instability strength increases with growing rotation frequency and field strength, and decreases with growing electron temperature and atomic weight. Possible reasons are suggested. Further research may focus on the comparison between these computations and experimental measurements on a well diagnosed MEVAT, and extend the analysis to other plasma thrusters where instability caused by plasma rotation and axial flow is of particular concern. These relevant thrusters may include helicon double layer thruster, 41 Hall thruster, 42 magnetized ion engine, 43 radio-frequency plasma thruster, 44,45 micro-wave (electron cyclotron resonance) plasma thruster, 46,47 and variable specific impulse magnetoplasma rocket (VASIMR), 48 to name a few.