Generalized lower-hybrid mode with, density gradient, equilibrium ExB drift, collisions and finite electron Larmor radius: Numerical studies with MATLAB solver
Ivan Romadanov, Andrei Smolyakov, Winston Frias, Oleksandr Chapurin, Oleksandr Koshkarov
GGeneralized lower-hybrid mode with, density gradient,equilibrium ExB drift, collisions and finite electron Larmorradius: Numerical studies with MATLAB solver.
Ivan Romadanov, ∗ Andrei Smolyakov, † Winston Frias, ‡ Oleksandr Chapurin, § and Oleksandr Koshkarov ¶ Department of Physics and Engineering Physics,University of Saskatchewan, Saskatoon SK S7N 5E2, Canada
Abstract
MATLAB solver for general dispersion relation has been developed to study local instabilitiesin partially magnetized plasmas typical for E × B discharges. Several examples of the Simon-Hoh, lower-hybrid and ion-sound instabilities with the parameters typical for Penning dischargeare shown. This solver allows parametrically investigate the local dispersion relation for the widerange of plasma parameters and prepare reports. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] a r X i v : . [ phy s i c s . p l a s m - ph ] M a r . INTRODUCTION Solver, ”Hall Plasmas Discharge Solver” (or ”HPDSolver”) has been developed to investi-gate local instabilities in partially magnetized plasmas typical for E × B discharge plasmassuch as in Hall thrusters, Penning discharges, and magnetrons. This solver has user inter-face, which allows easily change various plasma parameters and obtain results. The formof dispersion relation can be changed; thus, different physical effects can be included or ex-cluded and their effects on the unstable modes can be investigated. In this paper, examplesfor the Simon-Hoh, lower-hybrid, and ion-sound instabilities are presented, together withthe description of the code, and its features. Detailed explanations of the physics models,equations and notations are given in the companion paper in Ref. 1. The source code andexecutable installer can be downloaded by link 2. Current version of the installer can beused on the MS Windows 8 and 10 (the only two versions tested so far) computers and doesnot require the installed MATLAB. II. GENERAL DISPERSION RELATION FOR PENNING DISCHARGE
Important part of instability studies is an experimental investigation of plasma fluctu-ations and validation of theoretical and numerical models. Hall thruster is not convenientfrom experimental point of view, especially when probe measurements are performed. Thus,more flexible and versatile setup is required. Experiments with Penning discharges have longbeen demonstrating many instabilities which seem to be similar to those in Hall thrusters.One could expect that there are essential common physics elements between these two sys-tems. The geometry of Penning trap is similar to the geometry of Hall thrusters especially inlatest designs with strong axial magnetic field. The particular experimental design at Prince-ton Plasma Physics Laboratory[3] has the axial beam of energetic electrons that serves as acathode and the coaxial outer electrode as an anode. The applied voltage creates the radialelectric field and the axial magnetic field is created by the external solenoid, thus creatingcrossed E and B field configuration.Because of above advantages, theoretical investigation of Penning discharge seems to bea reasonable part of general study of E × B discharges, especially for verification of theproposed theories. Therefore, in this paper we use parameters of Penning discharge systemto investigate some predictions of local model. Typical parameters of the Penning dischargedevice[3] are presented in Table I. 2able I: Typical Penning discharge parameters Parameter Xe T e eVB GE V /mn m − L n -5 · − mc s . · m/sv T e . · m/sv ∗ . · m/sv . · m/sω LH . · rad/sω ce . · rad/sω ci . · rad/sω pe . · rad/sρ e . · − m where T e − electron temperature, B − magnetic field, n − density, E − electric field, L n − density gradient length, c s − ion sound velocity, v s − diamagnetic drift velocity, v − E × B velocity, ω LH − lower hybrid frequency, ρ e − electron gyroradius.Local models allow to investigate some characteristic features of the plasma dischargeswithout involving computation expensive and time consuming simulations. General disper-sion relation, which includes effects of electron-neutral collisions, density gradients, electroninertia and gyroviscosity, was developed [1, 4, 5]. For this model, E is in x − direction, B isin z − direction; y is a poloidal coordinate: ω ∗ − ω D + ( ω − ω + iν en ) k ⊥ ρ e + ik (cid:107) D (cid:107) ω − ω − ω D + ( ω − ω + iν en ) k ⊥ ρ e + ik (cid:107) D (cid:107) = k ⊥ c s + k (cid:107) c s ( ω − ω i ) (1)where ω ∗ − diamagnetic drift frequency, ω − E × B drift frequency, ν en − electron-neutralcollision frequency, ω i = k x V i - ion flow frequency, D (cid:107) - diffusion coefficient in z direction, k ⊥ = k x + k y - wave vector in the plane perpendicular to B , k (cid:107) - wave vector, parallel to B .3 II. MATLAB SOLVERA. User Interface
Hall Plasmas Discharge solver can be conveniently used to study parametric dependenceof local instabilities for various plasma parameters and different physical effects. Generalview of the software interface is shown in Fig. 1. It consists of three main fields: 1 – plasmaparameters; 2 – type of the solved dispersion relation; 3 – control buttons to calculate, plot,and save results.Figure 1: General view of the software interface. Section 1 – plasma parameters control;section 2 – dispersion relation form; section 3 – control buttons.At the startup all parameters are field with some initial values, user can change themaccordingly. Regular or scientific notation can be used for input values. Checkboxes on theright side allow to set required parameter as a variable. An example for the magnetic fieldis shown in Fig. 2.Figure 2: Example of selecting the magnetic field as a variable parameter.When user selects the checkbox, two additional fields appear: field for the final value ofthe parameter and field where user can select the number of points within the interval. It is4mportant to notice that all other checkboxes, except for k x and k y become disabled; thus,user is prevented from selecting other plasma parameters as variables. System length cannotbe selected as a variable parameter.An example of setting k y wavenumber as a variable is shown in Fig. 3. Meaning of fieldsis the same, but the number of point within the interval is bigger, thus a smooth dependencycan be obtained. If both k x and k y are selected, all checkboxes for plasma parameters willbecome disabled, thus user cannot set anything else as a variable.Figure 3: Example of selecting the k y wavenumber as a variable parameter.Either k x or k y wavenumber is the mandatory parameter to select. In case a user hitsthe ”Calculate” button without selection of k x or k y a warning sign will appear. Therefore,the next options for the output are available: a single 2D graph or a set of 2D graphs (incase one of the plasma parameter is chosen as a variable), which represent the dependencyof real and imaginary part of ω as a function of wavenumber (in x or y direction) at chosenset of plasma parameters; another option is two 3D surfaces, which represent dependency ofreal and imaginary part of ω as a function of both wavenumbers in x and y directions.Checkboxes in section 2 are used to control the form of the solved dispersion relation.Corresponding equation is shown at the bottom of the section 2, and it will be changedaccording to the selected parameters. One has to remember, that it is possible to modifythe dispersion relation in a way that no instability will occur.Control button are grouped in section 3. Once ”Calculate” button is pressed, the pop-upwindow with process bar will appear. After calculations are completed the window withplasma parameters will appear. It shows all velocities, plasma frequencies, and maximumvalues of ω together with corresponding wavenumber, at which this frequency is reached.Example is presented in Fig. 4. 5igure 4: Plasma parameter window.Plot functions are executed by ”Plot” button. As it was mentioned before, it is possibleto plot 2D and 3D graphs, depending on chosen parameters. ”PlotGV” button plot groupvelocity for the unstable modes. Output can be in 2D and 3D form as well. Group velocityis calculated numerically by v = ∂ω∂k (2)Save button generates two files. First file named as ”Figure name” field + ” plasma.txt”,contains calculation parameter fields and plasma parameters. It is a tab delimiter file, whichlater can be read by user from other applications. Files are saved in the same directory whereexecutable file is located. Example of the structure of this file is shown in Fig. 5.Figure 5: Plasma parameter window.First column is the parameter name, second is the value. If one of the parameter was6hosen as a variable, additional columns are created. Maximum values of omegas and cor-responding wavenumbers are saved as well.Second file is named as ”Figure name” field + ” data.txt”, it contains calculation re-sults, which is organized in columns: ”Parameter”, ”KY”, ”KX”, ”ROOT1”, ”ROOT2”,”ROOT3”. Parameters column is the string which contains the name of the field which washose as a variable. If no field was chosen then it is empty. ”KY”, ”KX” are wavenumbers.”ROOT” columns all roots for the dispersion relation. There can be 2 or 3 roots, dependingon the form of the equation. B. Installation
To install the software, one has to download ”HPDSolver web.exe” by link 2. During theinstallation it will automatically download and install all necessary files. However, if youalready have MATLAB run-time engine, you can download only ”HPDSolver.exe”, ”con-stants.mat”, and ”splash.png” files and run it directly. Remember, that ”constants.mat” isrequired to run the solver, otherwise it will crash.
IV. EXAMPLES OF INSTABILITIES FROM GENERAL DISPERSION RELA-TION
Here we present several cases of the unstable modes which can be found with the solver.All used parameters are taken from Table I.
A. Simon-Hoh instability
When the stationary electron current due to the E × B drift is present, the anti-driftmode becomes unstable due to the phase shift between the potential and electron density.Resulting gradient drift instability is described by the dispersion equation − ω ci k ⊥ L n k y = ω ω − ω , (3)where ω = k · V E is the azimuthal (closed drift) flow of electrons in crossed E × B fieldsand E = E (cid:98) x . This is the reactive instability of negative energy perturbations with aphase velocity below the stationary V E = E × B /B velocity. This mode is referred here as7he collisionless Simon-Hoh instability [6–8]. The condition E · ∇ n > ω ω ∗ > k ⊥ c s . (4)For given parameters, dependency of growth rate and frequency on k y ρ e are presented inFig. 6. k y ρ e ω , s - × Re( ω )Im( ω )k y c s k y v Figure 6: Growth rate and frequency of the Simon-Hoh mode as a function of k y ρ e , and k x = 1 m − .The instability condition is most easily satisfied for k x →
0. However, the modes with lowvalues of k x are not the fastest growing modes. In general, there exist multiple eigenmodeswith different values of k x ; roughly, the value of k x may be associated with a number of nodesin the x direction. For a given k y , the dependence of the growth rate on k x is nonmonotonous,see Fig. 7. k x ρ e ω , s - × Re( ω )Im( ω )k x c s k x v Figure 7: Growth rate and frequency of the Simon-Hoh mode as a function of k x ρ e , and k y = 1 m − .It follows from (4) that for high values of k x , the mode becomes stable. The maximumvalue of k x for unstable modes increases with k y . The growth rate is maximized when thefollowing condition is satisfied k y + k x = 2 ω ∗ ω c s = 2 E B k y ω ci L n . (5)8sing (5) one finds from (3) that the modes with the largest growth rate have the localeigenmode frequency given by the expression ω = ω + iω . (6)Therefore, for any given azimuthal wavenumber, k y , the real part of the frequency andgrowth rate for the most unstable local mode are equal to the E × B frequency.Solver interface with parameters, which were used to calculate data for Fig. 6, 7, is shownin Fig. 8. As it is seen, interface allows to control the form of the equation, what makes itmore convenient.Figure 8: View of the solver interface with parameters, which were used to obtain Fig. 6, 7. B. Effect of electron inertia on Simon-Hoh instability
As was noted in [1], Simon-Hoh instability is, in fact, can have high frequency comparableto ω LH even for low k y . Therefore, electron inertia effect become important and should beincluded. Corresponding dispersion relation is given in Eq. 7 ω ∗ + k ⊥ ρ e ( ω − ω ) ω − ω = k ⊥ c s ω . (7)Typical mode behavior is shown in Fig. 9. Electron inertia results in the cut-off for theinstability at high k y values, while dependence on k x remains approximately the same. Notethat the eigenmode frequencies become lower, compared to the Simon-Hoh case.Group velocity was calculated as well, results are shown in Fig. 10. It is seen that in y direction group velocity tend to be equal to v velocity.9 k y ρ e ω , s - × Re( ω )Im( ω )k y c s k y v k x ρ e ω , s - × Re( ω )Im( ω )k x c s k x v Figure 9: Frequency and growth rate of the Simon-Hoh instability with electron inertiaeffect as functions of k y ρ e at k x = 1 (left) and k x ρ e at k y = 1 (right). k y ρ e G r oup v e l o c it y , m / s × v g1 k y c s k y v k x ρ e G r oup v e l o c it y , m / s × v g1 k x c s k x v Figure 10: Group velocity in y (left) and x (right) directions.Solver interface with parameters, for Simon-Hoh instability with included inertia is shownin Fig. 11.Figure 11: View of the solver interface with parameters, which were used to obtain Fig. 9,10.10 . Effects of electron inertia and gyroviscosity on Simon-Hoh modes In general, modes with high k ⊥ ρ e (cid:29) ω ∗ + k ⊥ ρ e ( ω − ω ) ω − ω + k ⊥ ρ e ( ω − ω ) = k ⊥ c s ω . (8)Dependencies of unstable modes on k y ρ e and k x ρ e are shown in Fig. 12. In general,gyroviscosity shifts unstable region towards higher values of k y , and almost does not affectunstable region in x direction. k y ρ e ω , s - × Re( ω )Im( ω )k y c s k y v k x ρ e ω , s - × Re( ω )Im( ω )k x c s k x v Figure 12: Frequency and growth rate of the Simon-Hoh instability with electron inertiaan gyroviscosity effects as functions of k y ρ e at k x = 1 (left) and k x ρ e at k y = 1 (right).Group velocities in x and y directions are presented in Fig. 13. As before, maximumgroup velocity in y direction is defined by v velocity. k y ρ e G r oup v e l o c it y , m / s × v g1 k y c s k y v k x ρ e G r oup v e l o c it y , m / s × v g1 k x c s k x v Figure 13: Group velocity in y (left) and x (right) directions.Solver interface is presented in Fig. 14. 11igure 14: View of the solver interface with parameters, which were used to obtain Fig.12, 13. D. Full dispersion relation including effects of electron inertia, gyroviscosity, E × Bflow, and collisions
In systems with flows as typical for Hall discharge devices, there are modes with negativeenergy which become unstable due to the dissipation effects and collisions have been con-sidered as an additional destabilization mechanism. The dispersion relation with effects ofcollisions has the form ω ∗ + k ⊥ ρ e ( ω − ω + iν e ) ω − ω + k ⊥ ρ e ( ω − ω + iν e ) = k ⊥ c s ω . (9)Frequencies and growth rates as functions of k y ρ e and k x ρ e are shown in Fig. 15. Collisionsdo not affect values of maximum growth rates and frequencies a lot; however, they extendunstable region to higher wavenumbers. k y ρ e ω , s - × Re( ω )Im( ω )k y c s k y v k x ρ e -2-101234567 ω , s - × Re( ω )Im( ω )k x c s k x v Figure 15: Frequency and growth rate of the Simon-Hoh instability with electron inertiaan gyroviscosity effects as functions of k y ρ e at k x = 1 (left) and k x ρ e at k y = 1 (right).12roup velocities in x and y directions are presented in Fig. 16. Due to more smoothfunction of ω ( k y ), group velocity doe s not have an abrupt changes. The magnitude of groupvelocity is close to v . k y ρ e -4000-20000200040006000800010000 G r oup v e l o c it y , m / s v g1 k y c s k y v k x ρ e G r oup v e l o c it y , m / s × v g1 k x c s k x v Figure 16: Group velocity in y (left) and x (right) directions.Solver interface is presented in Fig. 17.Figure 17: View of the solver interface with parameters, which were used to obtain Fig.15, 16. E. Effect of magnetic field gradient
There is no magnetic field gradients in Penning discharge, however this solver can beused with Hall thruster parameters, where magnetic field gradients play an important role.Thus, here we show the effect of the magnetic field gradient on Simon-Hoh modes and onfull dispersion relation. Results are shown in Fig. 18.13 k y ρ e ω , s - × Re( ω )Im( ω )k y c s k y v Figure 18: Frequency and growth rate from full dispersion relation with magnetic fieldgradient as functions of k y ρ e at k x = 1.Magnetic field gradient stabilizes the instability. For selected parameter, there is noinstability for Simin-Hoh modes. Calculation for full dispersion relation shows that themagnitude of the maximum growth rate was significantly reduced. For the dependency on k x wavenumber it is too small and is not shown here.Solver interface is presented in Fig. 19.Figure 19: View of the solver interface with parameters, which were used to obtain Fig. 15. F. Search for growth rate and frequency as function of wavenumbers
It is possible to plot a dependency of frequency and growth rate for unstable mode asa function of both k y ρ e and k y ρ e . The results is two surfaces with represent frequency andgrowth rate. Example of such surfaces is given in Fig. 20.14igure 20: Frequency and growth rate.Group velocity can bu calculated as well. Result is shown in Fig.Figure 21: Group velocity.Solver interface is presented in Fig. 22.Figure 22: View of the solver interface with parameters, which were used to obtain Fig. 20.15 . CONCLUSION Simple and convenient solver for general dispersion relation, which describes plasma incrossed E and B field with unmagnetized ions was developed. It allows to make estimatesfor frequencies and growth rates of the unstable modes under various plasma parameters.Such estimations can be used for comparison with fluid or kinetic simulations, or in realexperiments.Several cases of unstable modes with different physical mechanisms were considered. Aninteresting property of the local modes is that the real part of the frequency and growthrate of the most unstable mode are equal to the local E × B frequency, ω = ω + iω forarbitrary value of the k y . [1] A. Smolyakov, O. Chapurin. W. Frias, O. Koshkarov, I. Romadanov, T. Tang, M. Umansky,Y. Raitses, I. Kaganovich, V. Lakhin, Plasma Physics and Controlled Fusion , 014041 (2017).[2] https://bitbucket.org/ivr509/hall-plasmas-discharge-solver/downloads [3] Y. Raitses, I. Kaganovich, and A. Smolyakov, ”Effects of the gas pressure on low frequencyoscillations in ExB discharges,” in Proceedings of Joint Conference of 30th ISTS, 34th IEPCand 6th NSAT, IEPC-2015-307.[4] W. Frias, A. Smolyakov, I. D. Kaganovich, and Y. Raitses, Physics of Plasmas , 011527(2013).[5] W. Frias, A. Smolyakov, I. D. Kaganovich, and Y. Raitses, Physics of Plasmas , 042105(2014).[6] A. Simon, Physics of Fluids , 382 (1963).[7] F. C. Hoh, Physics of Fluids , 1184 (1963).[8] Y. Q. Tao, R. W. Conn, L. Schmitz, and G. Tynan, Physics of Plasmas , 3193 (1994).[9] Y. Sakawa, C. Joshi, P. K. Kaw, F. F. Chen, and V. K. Jain, Physics of Fluids B-PlasmaPhysics , 1681 (1993)., 1681 (1993).