Structure of nonlocal gradient-drift instabilities in Hall E × B discharges Ivan Romadanov, Andrei Smolyakov, Yevgeny Raitses, Igor Kaganovich, Tang Tian, and Sergei Ryzhkov Citation: Physics of Plasmas 23, 122111 (2016); View online: https://doi.org/10.1063/1.4971816 View Table of Contents: http://aip.scitation.org/toc/php/23/12 Published by the American Institute of Physics Articles you may be interested in Tutorial: Physics and modeling of Hall thrusters Journal of Applied Physics 121, 011101 (2017); 10.1063/1.4972269 Modelling nonlinear electrostatic oscillations in plasmas Physics of Plasmas 23, 122103 (2016); 10.1063/1.4968520 A technique to control cross-field diffusion of plasma across a transverse magnetic field Physics of Plasmas 23, 122105 (2016); 10.1063/1.4971222 Dynamo-driven plasmoid formation from a current-sheet instability Physics of Plasmas 23, 120705 (2016); 10.1063/1.4972218 Parallel electric fields are inefficient drivers of energetic electrons in magnetic reconnection Physics of Plasmas 23, 120704 (2016); 10.1063/1.4972082 Modulation and nonlinear evolution of multi-dimensional Langmuir wave envelopes in a relativistic plasma Physics of Plasmas 23, 122112 (2016); 10.1063/1.4971444 PHYSICS OF PLASMAS 23, 122111 (2016) Structure of nonlocal gradient-drift instabilities in Hall E × B discharges Ivan Romadanov,1,a) Andrei Smolyakov,1,b) Yevgeny Raitses,2,c) Igor Kaganovich,2,d) Tang Tian,3,e) and Sergei Ryzhkov4,f) 1 Department of Physics and Engineering Physics, University of Saskatchewan, Saskatoon, Saskatchewan S7N 5E2, Canada 2 Plasma Physics Laboratory, Princeton University, Princeton, New Jersey 08543, USA 3 University of Science and Technology of China, Hefei, Anhui 230026, People's Republic of China 4 Bauman Moscow State Technical University, Moscow 105005, Russia (Received 10 June 2016; accepted 18 November 2016; published online 15 December 2016) Gradient-drift (collisionless Simon-Hoh) instability is a robust instability often considered to be important for Hall plasma discharges supported by the electron current due to the E × B drift. Most of the previous studies of this mode were based on the local approximation. Here, we consider the nonlocal model which takes into account the electron inertia as well as the effects of the entire profiles of plasma parameters such as the electric, magnetic fields, and plasma density. Contrary to local models, nonlocal analysis predicts multiple unstable modes, which exist in the regions, where local instability criteria are not satisfied. This is especially pronounced for the long wavelength modes which provide larger contribution to the anomalous transport. Published by AIP Publishing. [http://dx.doi.org/10.1063/1.4971816] I. INTRODUCTION Plasmas with crossed electric E and magnetic B fields are often used as ion sources in various applications. For moderate values of the magnetic field, the discharge is maintained by the electron E × B drift, while allowing the extraction, separation, and acceleration of unmagnetized ions. Such regimes, generally referred here as Hall plasmas, are widely employed in material processing (magnetrons1–3) and electric propulsion (Hall thrusters4,5) devices. The electron current across the magnetic field in such systems often exceeds the collisional value by orders of magnitude. Presumably, this anomalously large electron current is a result of plasma instabilities which are present in these devices for a wide range of plasma parameters and operational regimes. Despite long history of Hall plasmas applications, the understanding of the nature and mechanisms of plasma instabilities and plasma turbulence is still lacking. Among different types of unstable modes, relevant to E × B Hall plasmas, gradient-drift modes have long been considered as a possible source of fluctuations.6,7,18 These modes are closely related to the so called anti-drift mode8 that exists in Hall plasmas with density gradient, ω = k_⊥^2 c_s^2 / ω_*; (1) where c_s^2 = T_e/m_i is the ion-sound velocity and ω_* = k_y v_* = - k_y T_e / (e B_0 L_n) is the electron diamagnetic drift frequency due to density gradient; L_n is the characteristic length scale of the density gradient, L_n^-1 = ∇ ln n_0(x); k_⊥^2 = k_y^2 + k_x^2 is the wave vector perpendicular to the magnetic field, B_0 = B_0 ẑ. Note that the standard drift waves do not exist in Hall plasmas with unmagnetized ions. The frequency of the anti-drift mode (1) in fact does not depend on the electron temperature and the dispersion relation can also be written in the form ω = -ω_ci k_⊥^2 L_n / k_y, where ω_ci = e B_0 / m_i is the ion cyclotron frequency. The condition ω > k_z v_Te, where k_z is along the magnetic field lines, is required so that the mode is propagating almost perpendicular to the magnetic field. When the latter condition is not satisfied, the mode converts into the ion sound mode, ω^2 = k^2 c_s^2, which may propagate at a finite angle with respect to the magnetic field.9 When the stationary electron current due to the E × B drift is present, the anti-drift mode becomes unstable due to coupling with the ballistic mode ω = k_y V_0.10,11 Resulting gradient drift instability is described by the dispersion equation, - ω_ci k_⊥^2 L_n / k_y = ω^2 / (ω - ω_0); (2) where ω_0 = k · V_0 is the azimuthal (closed drift) flow of electrons in crossed E × B fields and E_0 = E_0 x̂. This is the reactive instability of negative energy perturbations with a phase velocity below the stationary V_0 = E × B / B^2 velocity. This mode is referred here as the collisionless Simon-Hoh instability.10–12 The condition E · ∇n > 0 is required for the instability.13 The gradient-drift instability was identified in early Hall thruster experiments6,14 as related to violent large scale structures (spokes). It was shown7 that properly profiled magnetic field improves the stability of the Hall thruster. The modified condition for the instability was derived in the form E · ∇(n/B) > 0. a) Electronic mail: [email protected] b) Electronic mail: [email protected] c) Electronic mail: [email protected] d) Electronic mail: [email protected] e) Electronic mail: [email protected] f) Electronic mail: [email protected] The gradient drift instability in plasmas with inhomogeneous magnetic field was revisited in Refs. 15 and 16 where the full compressibility was included resulting in the modified dispersion relation (ω - ω_D) / (ω - ω_0 - ω_D) = k_⊥^2 c_s^2 / ω^2; (3) where ω_D = k_y v_D = - 2 k_y T_e / (e B_0 L_B) is the magnetic drift velocity due to the axial gradient of the magnetic field; L_B^-1 = ∇ ln B_0(x) is the characteristic length of such a gradient. One consequence of this modification was that the magnetic field gradient enters in the combination ∇(n/B^2) rather than ∇(n/B) in the original work.7 In a finite temperature plasma, there also exists an additional contribution due to magnetic field gradient because of the compressibility of the electron diamagnetic velocity18 (the term with ω_D in the denominator on the left hand side). It was also shown in Ref. 15 that finite temperature fluctuations are also important in situations with a strong magnetic field gradient (complete dispersion relation for this case is given in Ref. 16). The collisionless Simon-Hoh instability has been considered among the most important sources of plasma fluctuations, in particular, as a cause of large scale structures or spokes.19–22 Much of the previous work on gradient-drift instabilities,23–25 see also other works cited in Ref. 23, was done within the local approximation neglecting the effects of the density and electric field profiles. The authors in Ref. 26 have considered the nonlocal solution to compare the eigenmode structure and eigenmode frequency with the data from Hall thruster experiments.26 However, the theoretical model used in Ref. 26 (following to Refs. 27 and 28) does not apply to the collisionless Simon-Hoh instability in which ions are not magnetized. The global analysis has been employed in Refs. 23, 29, and 30 for linear stability studies of Hall thruster for rather general fluid model that included effects of ionization, collisions, and heat flow. This model involved several different unstable modes and destabilization mechanisms; therefore, it is somewhat difficult to directly compare the results of the local and nonlocal analysis. The standard collisionless Simon-Hoh instability is considered to be a low frequency mode with the eigenfrequency well below the E × B frequency.16 However, it was shown recently17 that depending on the value of the k_x wave vector (along the electric field direction), the mode frequency becomes comparable to the E × B frequency, ω_0 = k · V_0, and thus, for the short-wavelengths, may become comparable with the low hybrid frequency ω_LH = (ω_ce ω_ci)^1/2 which required to be accounted for electron inertia. The goal of this work is to develop a nonlocal model for the collisionless Simon-Hoh instability and related higher frequency lower-hybrid modes and investigate the role of density and electric field profiles on the eigenmode structure, real part of the frequency, and growth rate in conditions typical for Hall plasma experiments. In case of the general profiles of the electric field and plasma density, the non-local stability problem is reduced to the eigenvalue problem which is solved here by the spectral method with Chebyshev polynomials.35 The results from the spectral code have been verified by a finite-difference differentiation matrix, and shooting method. Both spectral and finite-difference methods show good convergence and precision for low mode numbers. However, higher modes require significant increase in number of polynomials or number of steps for finite difference method to achieve the accurate solution. For our simulations, we used parameters relevant to two Hall plasma devices: Penning discharge and Hall thruster. Typical plasma parameters and profiles for such devices were taken from the experiments at Princeton Plasma Physics Laboratory.31–34 Penning discharge is a device with cylindrical geometry, in which the applied magnetic field is along the z-axis, and the electric field is in the radial direction. Plasma is created by the electron beam injected along the z-axis, from either filament cathode or plasma cathode such as a hollow cathode or RF plasma cathode. Magnetic field is created by Helmholtz coils. The experiments with the Penning discharges demonstrate range of instabilities and structures (spokes) similar to the Hall thruster33,34 and provide easier access for plasma diagnostics. Typical experimental parameters are: electric field E = -50 V/m, magnetic field B = 50 G, v_0 = 1 × 10^4 m/s, L_n = -0.1, ..., -0.04 m, system length L = 0.1 m, ω_ci = 3.7 × 10^3 s^-1, ω_LH = (ω_ce ω_ci)^1/2 = 1.8 × 10^6 s^-1. For calculations in Sections III A and III B, the parameters of Penning discharge were considered in slab geometry, with the magnetic field along the z-axis, and electric potential and density gradient along the x-axis. Hall thruster is another device of interest in which the gradient-drift instabilities play an important role. Typical device has coaxial or cylindrical, for the Cylindrical Hall Thruster (CHT),31,32 geometry. The discharge in the Hall thruster is created in the channel between the anode, which can be used as a gas distributor, and the cathode-neutralizer. A gas propellant, typically argon or xenon, is supplied to the channel through anode and then is ionized by high energy electrons from the cathode and as well electrons heated in the discharge. Magnetic field is sufficiently strong to make the electron Larmor radius much smaller than the chamber characteristic length thus creating the strong azimuthal E × B current of electrons. Collisions allow the electrons move to the anode in the axial direction, across the magnetic field and in the direction of the stationary electric field. The azimuthal current due to the E × B drift can be tens times larger than the axial electron current. The axial field accelerates the ions towards the channel exhaust, where they are neutralized by electrons from the cathode neutralizer. The thrust is a reaction force to this electrostatic acceleration, applied to the magnetic circuit. Even though azimuthal current is responsible for the ion thrust, axial current is important as well, because it completes the electric circuit. The axial electron current in the Hall thruster is typically an order of magnitude larger than the collisional value. This anomalous current is created by fluctuations driven by the electric field and gradients of the magnetic field, plasma density, and electron temperature which are strongly inhomogeneous along the thruster axis. For the Hall thruster configuration, the electric field, magnetic field, and density gradients are in the x–(axial) direction; the y-direction is periodic (azimuthal) direction. Typical magnetic field, electric field, and density profiles will be presented in Section IV. II. BASIC EQUATIONS Here, we formulate the eigenmode problem for the slab geometry with uniform axial magnetic field B_0 = B_0 ẑ, and inhomogeneous plasma density n_0 = n_0(x) and equilibrium potential ϕ_0 = ϕ_0(x). We consider the nonlocal modes in the form ϕ̃ = ϕ̃(x) exp(-iωt + i k_y y). The y direction is assumed to be periodic, a finite length interval is considered in the radial direction x. This geometry is an approximation for the cylindrical geometry of the Penning discharge with the axial magnetic field as in Ref. 31, in which case y is an azimuthal direction and x is radial. The linearized ion momentum equation has the form -iω m_i v_i = -e ∇ϕ̃. (4) Here, we assume cold ions and neglect the effect of the magnetic field on ion motion. The linear continuity equation for ions is -iω n_i + ∇ · (n_i v_i) = 0. (5) After substitution of v_i into the ion continuity equation, one finds -iω n_i - (i e / (m_i ω)) [n_0 ∇^2 ϕ + ∇n_0 · ∇ϕ] = 0. (6) In neglect of the electron inertia, electron motion along the magnetic field, and electron temperature, the electron velocity is a sum of the E × B drifts due to the stationary electric field E_0 = -∇ϕ_0, and the perturbed potential ϕ̃, with velocity ṽ_E = ∇ϕ̃ × B. The resulting density perturbation is ñ_e = -i (ṽ_E · ∇n_0) / (ω - ω_0). (7) Using the quasineutrality, one gets the eigenvalue equation in slab geometry ∂^2ϕ / ∂x^2 - k_y^2 ϕ + (n_0' / n_0) (∂ϕ / ∂x) - (ω^2 / (ω - ω_0)) (k_y / ω_ci) (∇n_0 / n_0) ϕ = 0, (8) where ω is an eigenvalue problem for prescribed boundary conditions on ϕ. In general case, ω and ϕ are complex. This eigenvalue problem is solved with the following boundary conditions: ϕ̃(a) = ϕ̃(-a) = 0 or ϕ̃(a) = ϕ̃(0) = 0, depending on the geometry type. It is convenient to use the transformation ϕ(x) = (1 / (n_0(x))^1/2) ψ(x). (9) Then Equation (8) is reduced to the form ∂^2ψ / ∂x^2 - k_y^2 ψ + [(ω^2 / (ω - ω_0)) (k_y / ω_ci) (n_0' / n_0) + (1/4) (n_0'^2 / n_0^2) - (1/2) (n_0'' / n_0)] ψ = 0. (10) For the constant gradient profile, n_0(x) = n_0 exp(x / L_n), (11) one has the following equation: ∂^2ψ / ∂x^2 - k_y^2 ψ + [(ω^2 / (ω - ω_0)) (k_y / (L_n ω_ci)) - (1 / (4 L_n^2))] ψ = 0. (12) For constant L_n and ω_0, from (12) one gets the following local dispersion relation: [FIG. 1. The growth rate as a function of k_x, and k_y = 10 m^-1; (a) for L_n = -0.1 m, and (b) for L_n = -0.04 m. Bold lines—growth rates from Eq. (13); thin dashed lines—growth rates from Eq. (2); squares—growth rates from Eq. (8).] (k_⊥^2 + 1 / (4 L_n^2)) (c_s^2 / ω^2) = ω_* / (ω - ω_0). (13) This equation, basically, is a modified Eq. (2). From this equation, one can see that density gradient increases the effective value of k_⊥ modifying the spectrum of unstable modes. This is illustrated in Fig. 1(a) for L_n = -0.1 m, and L_n = -0.04 m in Fig. 1(b). Profile of the ω_0 is uniform and system size L = 0.1 m. Bold lines are growth rate calculated from the local Eq. (13) with the effect of L_n. Thin dashed lines indicate the growth rate from Eq. (2). Density gradient effect leads to a slight increase of growth rates at low k_x values. Square marks are the discrete unstable growth rates obtained from numerical solution of Eq. (8) for the system size L = 0.1 m. It is important to stress that k_x values are not arbitrary and are defined by the size of the system. Another thing to note is that stronger density gradients (smaller absolute values of L_n) significantly increase a number of unstable modes and range of unstable values of k_x; however, the maximal value of the growth rate is not changed. Dependency of growth rate on k_y wavenumber with and without the effect of L_n is shown in Fig. 2. Results were obtained for L_n = -0.04 m. Consideration of L_n term leads to the constant shift of growth rate. Therefore, for moderate gradient values, this effect does not lead to significant change of the maximum growth rate. One can see from Figs. 1 and 2 that for the considered plasma parameters, the effect of (n_0' / n_0) (∂ϕ/∂x) term is not essential. However, Eq. (8) does not take into account the effect of electron inertia. The important result presented in Fig. 1 is that the maximum of the eigenfrequency is close to the ω_0 frequency. In fact, one can show17 that for continuous k_x, the eigenmode with the maximum growth rate is ω = ω_0 + iω_0. Therefore, it may become comparable with the lower-hybrid frequency and electron inertia has to be included. The respective equation can be written as in Ref. 17 (similar equation was also given in Refs. 37 and 38) ∂^2ϕ / ∂x^2 - k_y^2 ϕ - F(ω) (1 / (ω - ω_0)) (k_y / ω_ci) (∇n_0 / n_0) ϕ = 0, (14) where F(ω) = ω^2 / (1 - ω^2 / ω_LH^2). The full spectra of unstable solutions for two values of k_y with and without the effect of electron inertia are shown in Fig. 3. These spectra were obtained by solving Eqs. (8) and (14) for constant ω_0 profile and L_n = -0.04 m. For low values of k_y (empty circles and squares), unstable frequencies are much lower than ω_LH, so inertia effect does not play a significant role; however, for the case of high k_y values (solid circles and squares) electron inertia leads to the significant decrease of unstable frequency and growth rate values. Another critical effect of the electron inertia is the stabilization of the instability at large values of the k_y. These results are shown in Figs. 4(a) and 4(b). The results for the fixed value of the k_x are shown in Fig. 4(a); k_x = 31 m^-1 is the value chosen as a minimal possible for this system. As was shown above, the value of the k_x at which the growth rate is maximal changes as a function of k_y. The dependence of the maximum growth rate on the k_y is shown in Fig. 4(b). The value of k_x is different at each k_y value. Dependency on k_x obtained from Equations (8) and (14) is presented in Fig. 5 for two different k_y and L_n values. As before, stronger density gradients lead to the increase in the number of unstable modes. For low k_y values, inertia does not bring significant effect, but for higher k_y values inertia increases the range of unstable modes and causes the appearance of unstable modes with much higher k_x values. Therefore, electron inertia causes qualitative change in the behavior of unstable solution with the growth of k_y, bringing the cutoff of the instability at high k_x values. Changes in the dependency on k_x are more quantitative, and appear mostly at high k_y values. [FIG. 2. The growth rate as a function of k_y, for k_x = 31 m^-1 and L_n = -0.04 m. Bold line—growth rates from Eq. (2); thin dashed line—growth rates from Eq. (13).] [FIG. 3. The full spectrum of unstable solutions with (circles) and without (squares) electron inertia effect for L_n = -0.04 m: empty markers—k_y = 20 m^-1; filled markers—k_y = 100 m^-1.] III. NONLOCAL EIGENMODES FOR SOME MODEL PROFILES Here, we analyze the properties of the eigenmodes for different model profiles of E × B velocity. These simple profiles illustrate some characteristic features of the eigenmodes which also persist in more realistic cases. A. Step-like profile of the E × B velocity Eigenvalue problem (8) with a step-like profile of the E × B drift frequency ω_0 is shown in Fig. 6. Solution inside of each region can be obtained in the form of the harmonic functions, ϕ ~ A_i cos(λ_i x) + B_i sin(λ_i x), and the dispersion relation is easily obtained from matching conditions. A formal dispersion relation is not very illuminating. Here, we just emphasize the main features of such solutions and compare the nonlocal modes with the predictions of the local theory. First, we consider a case with the low k_y wavenumber values. As was noted above, the k_x range of the unstable modes is decreasing for lower values of k_y. Since the allowed wavenumber k_x is discrete, the number of unstable eigenmodes decreases with decreasing k_y. The effect of density gradient L_n on nonlocal unstable eigenmodes is similar to the effect of k_y; thus, the density gradient is fixed and taken L_n = -0.04 m for one wavenumber k_y = 20 m^-1. A characteristic feature of the ω_0 profile, shown in Fig. 6, is that the local theory predicts the instability only in the region -5 < x < 0, while the region 0 < x < 5 should be stable. The nonlocal solutions may exist in the locally stable region. Full spectra of unstable solutions are shown in Fig. 7(a). Three characteristic modes are chosen to show the shape of eigenfunctions in Fig. 7(b): point (1) with the lowest value of the real part of the frequency ω = (0.01 + 0.04i) · ω_LH s^-1; point (2) with the largest growth rate ω = (0.11 + 0.1i) · ω_LH s^-1; and point (3) with smallest growth rate ω = (0.14 + 0.01i) · ω_LH s^-1. All solutions have the nonlocal structure across the whole domain -5 < x < 5 including the region 0 > x > 5 which is locally stable. The tendency toward the local theory increases for larger k_y, as is expected. The full spectrum of eigenvalues in phase space for step-like profile of ω_0 with k_y = 100 m^-1 is shown in Fig. 8(a). Some characteristic eigenfunctions are shown in Fig. 8(b); one can see that there are solutions extending into the locally stable region. [FIG. 4. The growth rate as a function of k_y at L_n = -0.04 m. Dashed line—growth rate without electron inertia effect; solid line—growth rate with electron inertia effect; squares and circles—growth rates from nonlocal model without and with electron inertia effect, respectively. (a) Fixed k_x = 31 m^-1; (b) maximal growth rate at each k_y value, k_x changes with k_y.] [FIG. 5. The growth rate as a function of k_x for k_y = 20 m^-1 (a), and k_y = 100 m^-1 (b). Circles—results with the inertia; squares—results without the inertia. Filled markers—results for L_n = -0.1 m; empty markers—results for L_n = -0.04 m.] [FIG. 6. Profile of ω_0.] General tendency is that the solutions with the lower growth rates are more deeply extended into the locally stable region. B. Linear profile of the E × B velocity Here, we investigate the eigenmode problem for the parabolic potential profile ϕ = α x^2 corresponding to the case of a constant shear of the E × B velocity and exponential density profile n = n_0 exp(-x/a) which makes the density gradient length scale constant, L_n^-1 = n_0^-1 ∂n_0/∂x = const. Such profiles may occur in Penning discharge configurations where ions are confined radially by the inward radial electric field.34 We assume that ω_0 = k_y v_0 x / L, where L is the width of the region in the radial direction x. Several unstable eigenfunctions for k_y = 20 m^-1 and L = 20 cm are shown in Fig. 9. The full spectrum of the unstable solutions is presented in Fig. 9(a). For the E × B velocity with a constant shear, there exists a discrete spectrum of multiple unstable modes. Similar to the local theory results, the ground state is not the most unstable solution. The eigenmode with the largest growth rate for k_y = 20 m^-1 is shown in Fig. 10(b) with respect to the profile of the E × B frequency, ω_0(x). Lines with circles and squares represent the growth rate obtained from local theory for k_x = 0 and k_x = 63 m^-1, respectively. As it is seen, with the increase of k_x value, unstable region shifts to the right. The eigenfunctions with the largest growth rate are approximately localized close to the resonance point Re(ω) ≃ ω_0; however, there is an asymmetrical shift due to the finite growth rate. Similar to the local theory, the growth rate increases with the effective value of the radial wave number k_x. [FIG. 7. Results for step-like profile of the ω_0 with k_y = 20 m^-1 and L_n = -0.04 m. (a) The full spectra of unstable eigenvalues. (b) The eigenfunctions are for eigenvalues (from top to bottom): 1 - ω = (0.01 + 0.04i) · ω_LH s^-1, 2 - ω = (0.11 + 0.1i) · ω_LH s^-1, and 3 - ω = (0.14 + 0.01i) · ω_LH s^-1.] [FIG. 8. (a) Full spectrum of the eigenvalues for k_y = 100 m^-1. (b) The eigenfunctions are for eigenvalues (from top to bottom): 1 - ω = (0.07 + 0.25i) · ω_LH s^-1, 2 - ω = (0.71 + 0.12i) · ω_LH s^-1, and 3 - ω = (0.71 + 0.004i) · ω_LH s^-1.] [FIG. 9. (a) Full spectrum of the unstable eigenvalues; (b) unstable eigenfunctions for the constant E × B shear profile; wavenumber k_y = 20 m^-1. The eigenfunction with the largest growth rate, ω = (0.15 + 0.11) · ω_LH s^-1 - 1 (red); the ground state unstable eigenfunction with ω = (0.002 + 0.014i) · ω_LH s^-1 - 2 (blue); the eigenfunction with largest real frequency ω = (0.27 + 0.05i) · ω_LH s^-1 - 3 (black).] This is illustrated by the eigenfunction for the same value of k_y = 20 m^-1 as in Fig. 10(b) but for the extended domain with the length L = 20 cm, see Fig. 10(d). In the extended domain, the most unstable eigenfunction has larger growth rate, higher effective k_x, and the localization region shifts to the right (toward higher local ω_0), compare Figs. 10(b) and 10(d). Again, there is very little resemblance between the results of the local theory and the nonlocal solution. The local growth rates for comparison are shown in Figs. 10(b) and 10(d) for two different values of the radial wave number k_x. The local theory predicts that the mode with k_x = 0 is unstable in the wide region. For the k_x = 84 m^-1, the instability exists locally only in the region for x > 6 cm. In fact, the nonlocal solution has the effective k_x which is higher than k_x = 84 m^-1 and for which the local theory predicts no instability in the whole domain 0 < x < 20 cm. [FIG. 10. Left side: the eigenfunction with the largest growth rate ω = (0.09 + 0.6i) · ω_LH s^-1 (b), with respect to the profile of the E × B frequency (a), k_y = 20 m^-1 and L = 10 cm. The growth rates from the local theory for the same parameters and k_x = 0—squares (red), k_x = 63 m^-1—circles (blue). Right side: the eigenfunction with the largest growth rate ω = (0.07 + 0.06i) · ω_LH s^-1 (d) in the extended domain: L = 20 cm, for k_y = 20 m^-1, with respect to the E × B frequency profile (c). The growth rates from the local theory for the same parameters and k_x = 0—squares (red), k_x = 84 m^-1—circles (blue).] [FIG. 11. The eigenfunction with largest growth rate ω = (0.29 + 0.24i) · ω_LH s^-1, k_y = 100 m^-1, L = 10 cm.] The growth rate increases almost linearly with the wavenumber k_y, and the localization region becomes narrower around the resonance point Re(ω) ≃ ω_0. The most unstable eigenfunction for k_y = 100 m^-1 and system length L = 10 cm is shown in Fig. 11(b) with respect to the profile of the ω_0(x). Several unstable eigenfunctions in the extended domain L = 20 cm with the same value of k_y = 100 m^-1 are shown in Fig. 12. One can see that these modes are relatively local in the sense that they are weakly dependent on the boundary conditions at x = 0 and x = L. The eigenfunction with ω = (0.29 + 0.24i) · ω_LH s^-1 and localized at x ≃ 5 cm, shown in Fig. 11, corresponds to almost the same frequency Re(ω) and localized at the same x ≃ 5 cm as mode (2) which is shown in Fig. 12(b). However, localization of mode (1) is far from resonant point for this eigenmode, see line (1) in Fig. 12. This is due to the effect of the electron inertia. Therefore, another property of nonlocal modes can be understood from the example with the linear profile of ω_0. For unstable eigenmodes with absolute values of frequencies and growth rates well below of the ω_LH, the localization region is approximately determined by the point of the resonance Re(ω) = ω_0. For higher frequencies this is not true, in general. IV. NONLOCAL EIGENMODES IN PLASMAS WITH TEMPERATURE, DENSITY, AND MAGNETIC FIELD PROFILES Inhomogeneous magnetic field plays an important role in operation of Hall thrusters4,5 and magnetrons.1–3 Modification of the instability criteria due to the magnetic field was noted in Refs. 4 and 5. An additional effect of the magnetic field gradient occurs due to the finite temperature.15,18 The appropriate eigenvalue equation36 is ∂^2ϕ / ∂x^2 - k_y^2 ϕ + (1 / (1 - ω^2 / ω_LH^2)) (ω^2 (ω_* - ω_D) / ((ω - ω_0 - ω_D) c_s^2)) ϕ = 0. (15) For simulations, we use typical parameters from experiments in a Hall thruster.39,40 The density n_0, electron temperature T_e, electric field E, and magnetic field B profiles along the thruster axis (x-direction) are shown in Fig. 13. [FIG. 12. Several eigenfunctions are shown for the extended domain L = 20 cm, k_y = 100 m^-1. The eigenfunction with the largest growth rate for the extended domain, ω = (0.56 + 0.48i) · ω_LH s^-1 - 1 (blue); ω = (0.24 + 0.29i) · ω_LH s^-1 - 2 (red). The eigenfunction with the eigenvalue ω = (0.24 + 0.29i) · ω_LH s^-1 - 2 (red) has approximately the same frequency and localization as mode on with Fig. 11, for L = 10 cm.] [FIG. 13. Initial profiles of electric, magnetic fields, electron density, and electron temperature for Hall thruster experiment.39] [FIG. 14. Characteristic frequencies profiles and unstable solutions for k_y = 8.1 m^-1.] The profiles of the drift frequencies ω_0, ω_*, ω_D, and ω_LH are shown in Fig. 14. Calculations were done for parameters from Fig. 13 and k_y = 8.1 m^-1. For the complex profiles of plasma parameter, the nonlocal solution in general requires high number of polynomials in the spectral method. To verify the solutions obtained with spectral and shooting methods, and confirm the convergence, we have also used the integral relations that follow from Eq. (15) [FIG. 15. Full spectra of unstable solutions for Hall thruster profiles at k_y = 8.1 m^-1. Three groups of unstable solutions are marked as (a)–(c). Eigenvalues for numbered dots are presented in Table I.] TABLE I. Eigenvalues for Fig. 15. Fig. 15(a) 1: ω = (1.2 + 0.9) × 10^6 s^-1 2: ω = (1.0 + 0.6) × 10^6 s^-1 3: ω = (0.6 + 0.1) × 10^6 s^-1 Fig. 15(b) 1: ω = (-0.6 + 0.5) × 10^6 s^-1 2: ω = (-0.5 + 0.2) × 10^6 s^-1 3: ω = (-0.6 + 0.3) × 10^6 s^-1 Fig. 15(c) 1: ω = (-3.4 + 0.8) × 10^6 s^-1 2: ω = (-2.1 + 0.9) × 10^6 s^-1 3: ω = (-4.0 + 0.1) × 10^6 s^-1 [FIG. 16. Unstable eigenfunctions from each group of unstable solutions at k_y = 8.1 m^-1.] ∫ (|ϕ'|^2 + k_y^2 |ϕ|^2) dx = ∫ Re [ (1 / (1 - ω^2 / ω_LH^2)) (ω^2 (ω_* - ω_D) / (ω - ω_0 - ω_D)) ] (1 / c_s^2) |ϕ|^2 dx (16) and 0 = ∫ Im [ (1 / (1 - ω^2 / ω_LH^2)) (ω^2 (ω_* - ω_D) / (ω - ω_0 - ω_D)) ] (1 / c_s^2) |ϕ|^2 dx. (17) The ratio of left hand side and right hand side integrals was tracked with the increase in the number of polynomials. Desired accuracy was 10^-5. Similar to general trend, for a given k_y, there exist multiple eigenmodes with different growth rates and different localization regions. However, full spectra of unstable solutions look very different from spectra for simple profiles. An example of such spectra for k_y = 8.1 m^-1 is presented in Fig. 15. Three distinct groups of unstable solutions can be identified. An important difference is the sign of the real part of the frequency. There exist unstable modes with positive and negative frequencies. Note that the negative sign of the real part of the frequency correspond to the rotation in the E × B direction. Corresponding eigenvalues and eigenfunctions for the numbered point in each group are given in Table I and Fig. 16, respectively. Unstable eigenfunctions for marked points are shown in Fig. 16. Eigenfunctions, corresponding to different groups, have different localization. Contrary to the case of simple profiles, narrowly localized modes appear even for low k_y value. Modes in group (c) have almost exactly the same localization, thus local theory should work very well in this region. However, modes in group (a) are located in the region, where the local theory does not predict any instability (see in Fig. 20). [FIG. 17. Full spectra of unstable solutions for Hall thruster profiles at k_y = 81 m^-1. Four distinct groups of the unstable solutions can be identified.] [FIG. 18. Each group of unstable solution from Fig. 17. Eigenvalues for numbered dots are presented in Table II.] For higher value of k_y, full spectra of unstable solutions is modified, there appeared fourth group of unstable solutions, see Fig. 17. Group (d) appeared in between of groups (a) and (b) from Fig. 15. A detailed view of each group is presented in Fig. 18. Corresponding eigenvalues and eigenfunctions for the numbered point in each group are given in Table II and Fig. 19, respectively. Unstable eigenfunctions are presented in Fig. 19. Modes for groups (a) and (b) become local, as it is expected for higher k_y values. However, unexpected result is that the modes in groups (c) and (d) become highly nonlocal. Absolute values of frequencies are comparable in these two groups, but growth rates are lower in the region farther from the anode (group (d)). Moreover, the real and imaginary parts of the most unstable eigenvalue almost do not change with the k_y wavenumber. Growth rate of the unstable mode mostly depend on mode localization. Real part of frequencies for fastest modes is comparable with low-hybrid frequency ω_LH = 3.4 × 10^6 s^-1. The nonlocal modes (with the least number of nodes) have the lower growth rates compared to the localized (with the higher effective higher k_x) modes. The nonlocal unstable modes, presented in Figs. 16 and 19, are present in the regions where local modes are stable (shown in Fig. 20). It is important to note that the fast instabilities are absent in the acceleration region where the electric field is large. TABLE II. Eigenvalues for Fig. 18. Fig. 19(a) 1: ω = (3.0 + 0.6) × 10^6 s^-1 2: ω = (3.1 + 0.4) × 10^6 s^-1 3: ω = (2.8 + 0.4) × 10^6 s^-1 Fig. 19(b) 1: ω = (-0.9 + 0.23) × 10^6 s^-1 2: ω = (-0.8 + 0.13) × 10^6 s^-1 3: ω = (-1.0 + 0.1) × 10^6 s^-1 Fig. 19(c) 1: ω = (-2.9 + 0.14) × 10^6 s^-1 2: ω = (-2.4 + 0.03) × 10^6 s^-1 3: ω = (-3.4 + 0.02) × 10^6 s^-1 Fig. 19(d) 1: ω = (1.2 + 0.01) × 10^6 s^-1 2: ω = (1.0 + 0.17) × 10^6 s^-1 3: ω = (0.8 + 0.02) × 10^6 s^-1 [FIG. 19. Unstable eigenfunctions from each group of unstable solutions at k_y = 81 m^-1.] [FIG. 20. The growth rate from local model are shown as a function of distance along the thruster axis for k_y = 8.1 m^-1, and k_x = 0 m^-1, k_x = 100 m^-1.] The integral forms, Eqs. (16) and (17), were used to track the convergence. An example of the convergence for one group is shown in Fig. 21. Each marker and color represents different number of polynomials. Green markers are solutions for low polynomials number. With the increase in the number of polynomials, eigenvalues reach asymptotically the limit values which are shown as solid black circles. V. CONCLUSION Linear, nonlocal model for gradient-drift modes including the collisionless Simon-Hoh and lower-hybrid instabilities was developed and used to study the nonlocal structure of unstable mode in Hall plasma with inhomogeneous profiles of plasma density, temperature, electric and magnetic fields. The eigenvalue problem was solved numerically by using spectral and shooting methods and verified by the integral relations. Our model includes the collisionless Simon-Hoh destabilization mechanism13 E · ∇n > 0 and electron inertia that couples this mode to the lower-hybrid mode. A number of factors were not included in the model, such as collisional,41 ionization23 effects, and additional terms that explicitly include the shear flow effects.17,41 In general, neglect of these effects may limit the applicability of our model to some practical Hall devices. Our emphasis, however, was on the higher frequency modes with frequencies much higher than some typical low frequency (kHz range) oscillations, such as breathing mode due to ionization process. We have considered the eigenmode structure in the axial direction of the Hall thruster, and the radial direction for the Penning discharge. One of the important results of our analysis is demonstration of the existence of multiple modes with the eigenmode frequencies comparable to the electron drift frequency.42 There are significant differences between predictions of the local and nonlocal models, as it has been noted earlier.23,30 The nonlocal model for the step-like ω_0 profiles shows that for the low wavenumbers k_y nonlocal solutions can propagate into the region of the domain where local instability criteria are not satisfied. For higher k_y modes, the results of nonlocal model become similar to the local model where the unstable mode is highly localized. However, there are still unstable modes that extend into the region of stability and grow there even for high k_y values. Simulations for the E × B flow with a constant shear and constant density gradient show that there exist multiple unstable eigenmodes with different growth rates and different localization regions. Mode localization becomes more evident with increase of k_y. Moreover, such modes are almost independent on boundary conditions at x = 0 and x = L. We have also studied the eigenmode structure for complex profiles of the magnetic field, electric field, and electron temperature relevant to Hall thruster experiments.39 In general, the results demonstrate characteristic features similar to the outlined above. There exist multiple eigenmodes with different growth rates and different localization regions. The most unstable modes tend to localize in the region with higher gradients; however, modes with lower growth rates are present in the whole domain. There are significant differences between the predictions of the local and nonlocal models. This discrepancy is especially important for the modes with low k_y(m) and low effective k_x (longer wavelength). Often, the modes with largest growth rate are thought to dominate the dynamics on the longer time scales. However, in case of complex profiles, different eigenmodes may be localized in different spatial regions, even in regions which are locally stable. For example, in case of Hall thruster profiles there are modes which reside at the middle of the domain, where the local model predicts no instability. Some modes propagate from locally unstable regions to locally stable. The anomalous current in these conditions is due to the electron convection in the direction of the equilibrium electric field, ~⟨ñ Ṽ_x⟩.43 Then one can estimate the diffusion transport using a mixing length formula, D_a ≃ γ/k_y^2, where γ is the mode growth rate, k_y is the azimuthal wave number. Hence, significance of multiple modes is especially important for low k_y modes which have the nonlocal character and whose stability criteria could be very different from the predictions of the local theory and which provide larger contribution to the anomalous transport (D_a ≃ γ/k_y^2). These features of nonlocal solutions will be important for nonlinear transport calculations. In view of different stability criteria for multiple modes, especially for low k_y values, the predictions of the local theory may be misleading. ACKNOWLEDGMENTS The authors thank Ivan Halzov and Winston Frias for the help with numerical methods and model formulation, and Edward Startsev for the fruitful discussion. This work was supported in part by NSERC of Canada and U.S. Air Force Office for Scientific Research FA9550-15-1-0226. S.R. was partially supported by the Russian Ministry of Science and Education (Minobrnauka), Project No. 13.79.2014/K. [FIG. 21. Convergence of group (b) in Fig. 15 with the increase in the number of polynomials. Filled black circles—the stage when convergence is reached.] 1 P. Kelly and R. Arnell, Vacuum 56, 159 (2000). 2 T. Ito, C. V. Young, and M. A. Cappelli, Appl. Phys. Lett. 106, 254104 (2015). 3 S. N. Abolmasov, Plasma Sources Sci. Technol. 21, 035006 (2012). 4 A. I. Morozov, Plasma Phys. Rep. 29, 235 (2001). 5 A. I. Morozov and V. V. Savelyev, "Fundamentals of stationary plasma thruster theory," in Reviews of Plasma Physics (Springer US, Boston, MA, 2000), pp. 203–391. 6 A. I. Morozov, Y. V. Esipchuk, A. M. Kapulkin, V. A. Nevrovskii, and V. A. Smirnov, Sov. Phys. Tech. Phys. 17, 482 (1972). 7 Y. V. Esipchuk, A. I. Morozov, G. N. Tilinin, and A. V. Trofimov, Sov. Phys. Tech. Phys. 43, 1466 (1973). 8 A. M. Fridman, Sov. Phys. Dokl. 9, 75 (1964). 9 A. I. Smolyakov, W. Frias, I. D. Kaganovich, and Y. Raitses, Phys. Rev. Lett. 111, 115002 (2013). 10 A. Simon, Phys. Fluids 6, 382 (1963). 11 F. C. Hoh, Phys. Fluids 6, 1184 (1963). 12 Y. Q. Tao, R. W. Conn, L. Schmitz, and G. Tynan, Phys. Plasmas 1, 3193 (1994). 13 Y. Sakawa, C. Joshi, P. K. Kaw, F. F. Chen, and V. K. Jain, Phys. Fluids B 5, 1681 (1993). 14 L. A. Artsimovich, I. M. Andronov, Y. V. Esipchuk, I. A. Bersukov, and K. N. Kozubskii, Kosm. Issled. 12, 451 (1974). 15 W. Frias, A. Smolyakov, I. D. Kaganovich, and Y. Raitses, Phys. Plasmas 20, 052108 (2013). 16 W. Frias, A. Smolyakov, I. D. Kaganovich, and Y. Raitses, Phys. Plasmas 21, 062113 (2014). 17 A. Smolyakov, O. Chapurin, W. Frias, O. Koshkarov, I. Romadanov, T. Tang, M. Umansky, Y. Raitses, I. Kaganovich, and V. Lakhin, Plasma Phys. Controlled Fusion 59, 1, 14041 (2016); see http://stacks.iop.org/0741-3335/59/i=1/a=014041. 18 A. Kapulkin and M. M. Guelman, IEEE Trans. Plasma Sci. 36, 2082 (2008). 19 J. P. Boeuf and B. Chaudhury, Phys. Rev. Lett. 111, 847 (2013). 20 M. J. Sekerak, B. W. Longmier, A. D. Gallimore, D. L. Brown, R. R. Hofer, and J. E. Polk, IEEE Trans. Plasma Sci. 43, 72 (2015). 21 C. L. Ellison, Y. Raitses, and N. J. Fisch, Phys. Plasmas 19, 013503 (2012). 22 A. Kapulkin and M. Guelman, "Lower-hybrid instability in Hall thruster," in Proceedings of the 29th International Electric Propulsion Conference (Electric Rocket Propulsion Society, 2005). 23 D. Escobar and E. Ahedo, Phys. Plasmas 22, 102114 (2015). 24 Y. V. Esipchuk and G. N. Tilinin, Sov. Phys. Tech. Phys. 21, 417 (1976). 25 A. Kapulkin and M. Guelman, "Low frequency instability and enhanced transfer of electrons in near-anode region of Hall thruster," in Proceedings of the 30th International Electric Propulsion Conference (Electric Rocket Propulsion Society, 2007). 26 B. A. Jorns and R. R. Hofer, Phys. Plasmas 21, 053512 (2014). 27 R. F. Ellis, E. Marden-Marshall, and R. Majeski, Plasma Phys. 22, 113 (1980). 28 F. F. Chen, Phys. Fluids 10, 1647 (1967). 29 S. Barral, Y. Jayet, S. Mazouffre, M. Dudeck, E. Vron, and P. Echegut, "Hall effect thruster with an AlN discharge channel," in Proceedings of the 29th International Electric Propulsion Conference (Electric Rocket Propulsion Society, 2005). 30 D. Escobar and E. Ahedo, IEEE Trans. Plasma Sci. 43, 149 (2015). 31 Y. Raitses, P. Baele, and V. M. Donnelly, in Proceedings of the 64th Gaseous Electronics Conference (2011). 32 A. Smirnov, Y. Raitses, and N. J. Fisch, J. Appl. Phys. 95, 2283 (2004). 33 Y. Raitses, J. B. Parker, E. Davis, and N. J. Fisch, "Background gas pressure effects in the cylindrical Hall thruster," in Proceedings of the 46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit (American Institute of Aeronautics and Astronautics, 2010), Paper No. AIAA-2010-6775. 34 Y. Raitses, I. Kaganovich, and A. Smolyakov, "Effects of the gas pressure on low frequency oscillations in E × B discharges," in Proceedings of the Joint Conference of 30th ISTS, 34th IEPC, and 6th NSAT (2015). 35 L. N. Trefethen, Spectral Methods in MATLAB (Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2000). 36 W. Frias, A. I. Smolyakov, I. D. Kaganovich, and Y. Raitses, Phys. Plasmas 19, 072112 (2012). 37 A. M. DuBois, T. Edward, W. Amatucci, and G. Ganguli, Phys. Plasmas 21, 062117 (2014). 38 G. Ganguli, Y. C. Lee, and P. J. Palmadesso, Phys. Plasmas 31, 823 (1988). 39 Y. Raitses, D. Staack, M. Keidar, and N. J. Fisch, Phys. Plasmas 12, 057104 (2005). 40 A. A. Litvak and N. J. Fisch, Phys. Plasmas 8, 648 (2001). 41 A. A. Litvak and N. J. Fisch, Phys. Plasmas 11, 1379 (2004). 42 A. Lazurenko, V. Krasnoselskikh, and A. Bouchoule, IEEE Trans. Plasma Sci. 36, 1977 (2008). 43 R. Spektor, "Quasi-linear analysis of anomalous electron mobility inside a Hall thruster," in Proceedings of the 30th International Electric Propulsion Conference (Electric Rocket Propulsion Society, 2009).