University of Southern California Dissertations and Theses

Kinetic studies of collisionless mesothermal plasma flow dynamics

(USC Thesis Other)

Kinetic studies of collisionless mesothermal plasma flow dynamics

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content KINETICSTUDIESOFCOLLISIONLESSMESOTHERMAL PLASMAFLOWDYNAMICS by Yuan Hu A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY (ASTRONAUTICAL ENGINEERING) August 2018 Copyright 2018 Yuan Hu Epigraph ii Epigraph If I have seen further, it is by standing on the shoulders of giants. – Isaac Newton, 15 February 1676 Dedication iii Dedication To my beloved... Acknowledgments iv Acknowledgments I am eternally grateful to the large number of mentors family, friends and colleagues who helped me along my path toward earning my Ph.D.. Only with their inspiration, advice, support and help was I capable of pushing forward and stretching my limits farther than I ever thought was possible. First, I would like to express my deepest gratitude to my advisor, Professor Joseph J. Wang, for his valuable guidance and advice throughout my entire Ph.D. life, for his generous support and constant encouragement, and for his trust that allows me to purse my own research ideas while providing constructive academic criticism at the same time. In addition to all the academic research support provided by him, I would also like to thank him for always ensuring that my long term career goals are a high priority. I would like to thank the other members of my dissertation committee for their useful help and suggestions. I am deeply grateful to Professor Aiichiro Nakano for all the knowledge of computational physics and parallel computing from his classes, for his insightful guidance and continuous support on my dielectric materials microwave heating research, and for stimulating me pursuing further on the way of academia with curiosities of life and science. I sincerely appreciate Professor Daniel A. Erwin for his very helpful comments on my dissertation and valuable suggestions on how to pursue the career in academia. I also wish to acknowledge ProfessorJosephKuncandProfessorMikeGruntmanfortheirgenuinequestionsandbeneficial insights on my ongoing and proposed researches during my qualifying exam. I would like to thank my colleagues, Dr. Daoru Han, Dr. Dayung Koh, Dr. Scott Hughes, Dr. Kevin Chou, Dr. William Yu, Mr. Yinjian Zhao, Mr. Martin Hilario, Mr. Daniel Depew, Acknowledgments v Mr. Chen Cui, and Mr. Robert Antypas and Mr. Jeffrey Asher. Their support kept me motivated and our discussions over the years brought to light ideas that would have otherwise remained in the shadows. I am grateful for all of the hard work of the staff of the Department of Astronautical Engineering, Dell Causon, Linda Ly, Marlyn Lat, Nicole Valdez, Norma Orduna, Marrietta Penoliar, and Ana Olivares. Only through their efforts can the department function and can the students thrive. The generous support by the Viterbi School of Engineering Ph.D. Fellowship from the University of Southern California (USC) is gratefully acknowledged. The Center for High- Performance Computing of USC is sincerely appreciated for providing the computational resources for the numerical simulations carried out in this dissertation. Finally, I would like to express my greatest appreciation to my father and mother for their unconditional and endless support to raise the first Dr. Hu in our family. I also wish to thank all other family members for their constant care throughout my whole life. Most importantly, there are no words that can express my gratitude for my wife, Wenjuan, for her love, understanding, and support, and for her extraordinary courage and efforts to be the greatest mother to raise our daughter, Sijia (Scarlett). I owe both of you a lot. Abstract vi Abstract The collisionless, mesothermal plasma flow dynamics is a generic plasma physics problem having many potential applications. In space engineering, the plasma-spacecraft (S/C) interactions and the expansion of a plume from an electric propulsion (EP) thruster are two subjects in which the collisionless, mesothermal plasma flow dynamics plays an very important role. To understand the electron kinetics and validate the commonly used electron fluid approximations in such problems, this dissertation carries out comprehensive fully kinetic particle-in-cell (PIC) simulations as well as hybrid fluid-electron/particle-ion PIC simulations. The plasma-S/C interactions are first studied by considering a collisionless, mesothermal plasma flow over a large unbiased conducting plate. The hybrid PIC modeling first employs the most widely used Boltzmann relation, in which the electron is approximated by a massless, isothermal fluid, to describe the electron dynamics. The validity of this model is investigated by directly comparing the results between the hybrid and full PIC simulations. The comparison shows that the plasma wake may be characterized into two regions based on electron behavior, a fluid expansion region and a kinetic electron expansion region. In the fluid electron expansion region, the approaches based on the electron fluid approximation and the fully kinetic model lead to a similar result. In the kinetic electron expansion region, the fluid approximation for electrons breaks down and the fully kinetic approach must be adopted. The results also show that the error of the hybrid PIC simulation strongly correlates with a new non-equilibrium parameter ( ) based on the weighted deviation of the actual electron velocity distribution from the one at the equilibrium. The capability of the hybrid PIC modeling that incorporates the more general polytropic electron fluid approximation is Abstract vii then examined. The results predicted by this model exhibit very little improvement over the Boltzmann electron based hybrid PIC method. The influence of surface charging on the plasma-plate interactions is investigated next. In particular, the plasma wake structures, ion impact and current collection on the plate are studied with the Boltzmann electron based hybrid PIC and the full PIC simulations. The hybrid PIC modeling is shown to accurately describe the ion dynamics in the vicinity of the plate, but to fail downstream. The hybrid PIC failure is primarily due to the breakdown of electron fluid approximation, confirmed by the local non-equilibrium degree of electron flow measured by the non-equilibrium parameter . The EP plume problem is studied first with a 2D plume injected from a pre-neutralized source by using both the hybrid and full PIC simulations. The full PIC simulations are carried out using the real ion-to-electron mass ratios of H + , Ar + and Xe + , while the hybrid PIC model employ the Boltzmann relation for the electrons. The simulations reveal that the plume structure exhibits four distinct expansion regions. These regions are identified as, from the plasma emission plane to the downtream along the plume direction, the unperturbed, quasi- steady expansion, self-similar expansion and electron front regions, respectively. The behavior of electrons is strongly anisotropic, causing considerably different expansion characteristics betweentheplumedirectionandthetransversedirection. The3Deffectsonthepre-neutralized plume expansion is further examined. The results show that the 2D plume is qualitatively similar to the 3D plume, but quantitative differences exist. The comparison between the hybrid and full PIC results shows that modeling electrons either as an isothermal fluid or as a polytropic fluid with an isotropic temperature is over-simplified and leads to not only quantitative but also qualitative differences. In order for an electron fluid approximation to be applicable, the model must at least have the capability of resolving the strong anisotropic thermodynamics for electrons. Finally, a fully kinetic PIC simulation is carried out to study the emission of a 2D collisionless plasma plume. All the important processes, including the ion-electron coupling, beam neutralization, and plume expansion, relevant to the beam plasma emission from an Abstract viii ion thruster are taken into account in this model. The macroscopic plume structure is shown to have several distinctive regions, including an undisturbed core region, an electron cooling expansion region, and an electron isothermal expansion region. The properties of each region are determined by the microscopic electron kinetic characteristics. TABLE OF CONTENTS ix Table of Contents Epigraph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xvii List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xviii Chapter 1:Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Significance of Collisionless, Mesothermal Plasma Flows . . . . . . . . 2 1.1.2 Necessity of Kinetic Investigations . . . . . . . . . . . . . . . . . . . . 5 1.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.1 Expansion of a Collisionless Semi-infinite Plasma . . . . . . . . . . . 7 1.2.2 Expansion of a Collisionless Finite Size Plasma . . . . . . . . . . . . 12 1.2.3 Collisionless, Mesothermal Plasma Flow over an Object . . . . . . . . 15 1.2.4 Electric Propulsion Plasma Plume . . . . . . . . . . . . . . . . . . . . 17 1.3 Motivation and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.4 Dissertation Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Chapter 2:Kinetic Simulation Method . . . . . . . . . . . . . . . . . . . . . . 24 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2 Particle-in-Cell Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.2.1 Fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 TABLE OF CONTENTS x 2.2.2 Full PIC vs. Hybrid PIC . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2.3 Normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.3 Remarks on Full Particle Simulation . . . . . . . . . . . . . . . . . . . . . . 35 2.3.1 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.3.2 Structured Stretching Mesh . . . . . . . . . . . . . . . . . . . . . . . 36 2.3.3 Poisson Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Chapter 3:Collisionless, Mesothermal Plasma Flow over a Large Unbiased Conducting Plate: Breakdown of the Electron Fluid Approxi- mation in the Wake . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 Formulation and Analytical Solution . . . . . . . . . . . . . . . . . . . . . . 43 3.3 Simulation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.4 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4.1 Plasma field: fluid electron vs kinetic electron . . . . . . . . . . . . . 47 3.4.2 Kinetic characteristics and fluid approximation breakdown of electron 50 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Chapter 4:Collisionless, Mesothermal Plasma Flow over a Large Biased Conducting Plate: Effects of Surface Charging on the Electron Fluid Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.2 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.3.1 Plasma field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.3.2 Ion impact and current collection . . . . . . . . . . . . . . . . . . . . 64 4.3.3 Kinetic characteristics and fluid approximation breakdown of electron 66 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Chapter 5:Expansion of a Collisionless Hypersonic Plume into Vacuum: Fully Kinetic and Hybrid Particle-in-Cell Simulations . . . . . . 71 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.2 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.3.1 Plasma plume structure . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.3.2 Ion dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.3.3 Effects of electron characteristics . . . . . . . . . . . . . . . . . . . . 86 5.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 TABLE OF CONTENTS xi Chapter 6:ExpansionofaCollisionlessHypersonicPlumeintoVacuum: 3D Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.2 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.3.1 Evolution of plasma plume . . . . . . . . . . . . . . . . . . . . . . . . 93 6.3.2 Steady state plasma plume . . . . . . . . . . . . . . . . . . . . . . . . 95 6.3.3 Electron thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . 98 6.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Chapter 7:Fully Kinetic Simulation of Collisionless, Mesothermal Plasma Emission: Macroscopic Plume Structure and Microscopic Elec- tron Characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 7.2 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.3.1 Macroscopic Plume Structure . . . . . . . . . . . . . . . . . . . . . . 106 7.3.2 Microscopic Electron Characteristics . . . . . . . . . . . . . . . . . . 112 7.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 Chapter 8:Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 117 8.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 8.1.1 Collisionless, Mesothermal Plasma Flow over a Large Conducting Plate117 8.1.2 Expansion of a collisionless plasma plume into a vacuum . . . . . . . 119 8.2 Recommendations for Further Study . . . . . . . . . . . . . . . . . . . . . . 122 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 LIST OF FIGURES xii List of Figures 1.1 Artist’s rendition of the plasma iteracting with a space station. . . . . . . . . 2 1.2 Artist’s concept of a solar electric propulsion system. . . . . . . . . . . . . . 4 2.1 Schematic of the distribution function loosely represented by the particles in a cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2 Flow schematic for the ES-PIC scheme. . . . . . . . . . . . . . . . . . . . . . 28 2.3 Schematic of force weighting and particle weighting based on the linear weight- ing scheme in 2D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.4 Illustration of a stretching mesh being used for the plasma plume simulation in this dissertation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.1 Schematics for (a) plasma wake and (b) simulation setup. . . . . . . . . . . . 43 3.2 Instantaneous ion density profiles along thex-axis from the Boltzmann electron fluid based hybrid PIC simulation at ~ y =y= D0 = (a) 0, (b) 35, (c) 70 and (d) 105. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3 Instantaneous ion density profiles along the x-axis from the full particle PIC simulation at ~ y =y= D0 = (a) 0, (b) 35, (c) 70 and (d) 105. . . . . . . . . . 46 3.4 Plasma wake. Left panel: (a) ion density and (b) electron density. The color lines are the iso-density lines from the Boltzmann electron fluid based hybrid PIC simulation. The color contours are the density errors of the hybrid PIC results against the full PIC results defined by Eq. (3.4). The gray dashed lines the analytical solution of Eq. (3.3). Right panel: potential contours (color), electric field vectors (yellow arrows) and ion streamlines (black solid lines) from the (c) Boltzmann electron fluid based hybrid PIC and (d) full particle PIC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.5 (a) Ion density ~ n i contours from the hybrid PIC simulations with various values and the full particle PIC simulation. (b) Weighted root mean square error of the hybrid PIC results against the full PIC result defined by Eq. (3.5) as a function of . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.6 Electron temperature (color contours) and streamlines (black solid lines) from full PIC simulation: (a) ~ T e;x and (b) ~ T e;y . The magenta solid line shows the boundary between electron forward-flow and back-flow. . . . . . . . . . . . 51 LIST OF FIGURES xiii 3.7 Local electron VDFs ( ^ f S ) from the full PIC simulation (color surface plots) and fitted Maxwellian VDFs ( ^ f M ) based on local ~ v d and ~ v t (black solid lines). (a) ^ f(~ v x ; ~ x) along ~ y = 0, (b) ^ f(~ v x ; ~ y) along ~ x = 0, (c) ^ f(~ v y ; ~ x) along ~ y = 0, and (d) ^ f(~ v y ; ~ y) along ~ x = 0. The insets highlight the comparison between the actual VDF ( ^ f S ) and the fitted Maxwellian ( ^ f M ) at several representative positions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.8 (a) Contour of non-equilibrium parameter overlaid with the iso-lines of ~ n i;err = 10% (pink dotted) and ~ n e;err = 10% (black solid). (b) Dependence of the conditional averages of hybrid PIC prediction errors on . . . . . . . . . 53 4.1 Schematic of the plasma field around a large biased plateform in the (a) “quasi-neutral”, (b) “thin sheath”, and (c)“thick sheath” regimes [1]. . . . . . 57 4.2 Schematics of the simulation setup for a collisionless, mesothermal flow over a negatively biased large thin conducting plate. . . . . . . . . . . . . . . . . . 59 4.3 Ion and electron densities (color isolines) from the hybrid PIC simulations, and density errors of the hybrid PIC results against the full PIC results (color contours) defined by Eq. (3.4). . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.4 Contours of normalized potential ~ = ~ p (color), electric field vectors (yellow arrows) and ion streamlines (black solid lines) from both hybrid and full PIC simulations. The region enclosed by the white line is where ~ < ~ p , indicating that there is a potential well in the wake. The potential well structure in the wake appears only for the cases of ~ p =2 and5. . . . . . . . . . . . . . 63 4.5 Ion current density profiles on the front and back faces of the plate. . . . . . 65 4.6 Effects of surface potential on: (a) ion current collection, (b) Wake side current ratio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.7 Electron temperatures (color contours) and streamlines (black solid lines) from full PIC simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.8 Contours of non-equilibrium parameter overlaid with the isolines of ~ n i;err = 10% (pink dotted) and ~ n e;err = 10% (black solid). . . . . . . . . . . . . . . . 69 5.1 Schematics of the simulation setup for a collisionless hypersonic plume expan- sion into a vacuum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.2 One-dimensional electric potential profiles along the x-axis (plume direction) inside the core region fort! pi = 30 (blue), 40 (red) and 50 (black). The profiles are extracted along ~ y = (a) 0, (b) 10, and (c) 18. . . . . . . . . . . . . . . . 76 5.3 (Color) Iso-contours of the electric potential att! pi = 50 for the full PIC results with (a) H + , (b) Ar + , (c) Xe + , and (d) the hybrid PIC result. Yellow arrows show the electric field vectors. Black dashdot lines (“ ”) are the calculated first Mach lines. The inner white and outer green solid lines correspond to the ion density iso-lines with the values of ~ n i (~ x = 50[~ v 0 =C s0 1] + ~ x s ; ~ y = 0) and ' 0, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 LIST OF FIGURES xiv 5.4 (Color) Iso-contour lines of the ion (“ ”) and electron (“ ”) number density at t! pi = 50 for the full PIC results with (a) H + , (b) Ar + , (c) Xe + , and (d) the hybrid PIC result. Black dashdot lines (“ ”) are the calculated first Mach lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.5 One-dimensional electron temperature profiles inside the core plume region for t! pi = 50 extracted along ~ y = (a) 0, (b) 10, and (c) 18. . . . . . . . . . . . . 79 5.6 Fitting value of the normalized slope of the electric potential profile in the self-similar expansion region as a function of time. . . . . . . . . . . . . . . . 80 5.7 Ion velocity contours at t! pi = 50 for the full PIC results with (a) H + , (b) Ar + , (c) Xe + , and (d) the hybrid PIC result. The x-component and the y-component ion velocities normalized by C s0 are shown on the left panel (1) and the right panel (2), respectively. . . . . . . . . . . . . . . . . . . . . . . 82 5.8 Velocity increase of fastest ions as a function of time: (a) v max i;x =C s0 M 0 = ~ v max i;x , (b) v max i;y =C s0 = ~ v max i;y . The present simulation results from different models are represented by symbols. The red solid line is the prediction of Eq. (5.3) from Mora’s work [2]. The dotted lines are the curves best fitted to the full PIC simulation results with different ion species according to Eq. (5.4). 83 5.9 Dependence of final velocity increase of the fastest ions on mass ratio. . . . . 85 5.10 One-dimentional electron temperature profiles of H + beams for beam sizes R 0 = 20 D0 and 40 D0 att! pi = 50. The profiles are extracted along the beam direction at the position 20 D0 away from the beam edge, namely along ~ y = 0 (beam center) for R 0 = 20 D0 and ~ y = 20 for R 0 = 40 D0 , respectively. . . . 87 6.1 Schematics of the simulation setup for a collisionless hypersonic plume expan- sion into a vacuum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.2 One-dimensional electric potential profiles along the x-axis (plume direction) inside the core region for t! pi = 30 (blue), 40 (red) and 50 (black). The upper (a, b, c) and lower (d, e, f) panles show the results from the hybrid PIC and full PIC simulations, respectively. The profiles are extracted along ~ r = (a, d) 0, (b, e) 10, and (c, f) 18. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.3 (Color) Iso-contour lines of the ion (“ ”) and electron (“ ”) number density at the steady state: (a) two-dimensional hybrid PIC, (b) axi-symmetric hybrid PIC, (c) two-dimensional full PIC and (d) axi-symmetric full PIC. . 95 6.4 One-dimensional plasma density profiles along the x-axis (plume direction) inside the core region for the steady state results, a comparison between the two-dimensional and axi-symmetric results from both the hybrid and full PIC simulations. The upper (a, b, c) and lower (d, e, f) panles show the ion and electron densities, respectively. The profiles are extracted along ~ r = (a, d) 0, (b, e) 10, and (c, f) 18. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.5 Ion streamlines originated from the injection plane ~ x = 6 for ~ r = 2, 10, and 18, a comparison between the two-dimensional and axi-symmetric results from both the hybrid and full PIC simulations. . . . . . . . . . . . . . . . . . . . . 96 LIST OF FIGURES xv 6.6 (Color) Iso-contour lines of the ion (“ ”) and electron (“ ”) number density at the steady state: (a) two-dimensional hybrid PIC, (b) axi-symmetric hybrid PIC, (c) two-dimensional full PIC and (d) axi-symmetric full PIC. . 97 6.7 One-dimensional electron temperature profiles along the x-axis (plume direc- tion) inside the core region for the steady state full PIC results, a comparison between the two-dimensional and axi-symmetric results. The profiles are extracted along ~ r = (a) 0, (b) 10, and (c) 18. . . . . . . . . . . . . . . . . . . 97 6.8 log 10 ( ~ T e ) log 10 (~ n e ) along the ion streamlines originated from the injection plane ~ x = 6 at various ~ r positions: (a) 2, (b) 4, (c) 6, (d) 8, (e) 10, (f) 12, (g) 14, (h) 16, and (i) 18. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 7.1 SimulationsetupforthefullykineticPICmodelingofcollisionless, mesothermal plasma emission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 7.2 Electric potential contours (a) and total charge density contours (b) at t! pi0 = 40 (t! pe0 = 1713:6). The values of ~ plotted are with respect to the ambient reference potential. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 7.3 Number density contours (normalized by ~ n i0 ) for ions (colored “ ”) and electrons (colored “ ”) at t! pi0 = 40 (t! pe0 = 1713:6) for m i =m e = 1836 and M 0 = 4:29. The black “ ” represents the Mach line based on M 0 originating from the edge of the emission source, (~ x 0 ; ~ z 0 ) = (20; 6). . . . . . . 106 7.4 Contours of electron temperature ~ T e;x (a) and ~ T e;z (b) and four selected ion streamlines at t! pi0 = 40 (t! pe0 = 1713:6). The selected ion streamlines originated at (~ x 0 ; ~ z 0 ) = (2; 6); (6; 6); (12; 6); (16; 6). The boundaries for different plume regions are also shown. The iso-potential line with a value equal to the average at the plume emission plane ( 9:125) is also shown. . . . . . . . . . 108 7.5 1-D profiles of potential (a), electron number density (b), and electron tem- perature (c) along the z direction for ~ x = 0, 1 2 R b , and R b at t! pi0 = 40 (t! pe0 = 1713:6). The division between the different regions is based on the ~ x = 0 profile. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 7.6 log 10 (T e ) log 10 (n e ) relation along selected ion streamlines at t! pi0 = 40 (t! pe0 = 1713:6). The T e vs. n e points in Region B are clustered together at the end of the curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.7 log 10 (T e )log 10 (n e ) relation (a) - (h) and variations of polytropic coefficients (i) with different ion streamlines at t! pi0 = 40 (t! pe0 = 1713:6). The T e vs. n e points in Region B are clustered together at the end of the curve. . . . . . . 111 7.8 EVDFs at selected locations along the central axis of the emitting source at t! pi = 40. ^ f(~ v x ) (“ ”) and ^ f(~ v z ) (“ ”) are the instant electron velocity distributions sampled from the macro-particles at that time step. The EVDFs are calculated from the macro-particles contained in one cell. The Maxwellian velocity distribution fits using the temperatures and velocities from the macro- particles at each location, ^ f M (~ v x ) (“4”) and ^ f M (~ v z ) (“”), are also plotted for comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 LIST OF FIGURES xvi 7.9 EVDFs at selected locations along two off-center axis, ~ x = 20 and ~ x = 35, at t! pi = 40. ^ f(~ v x ) (“ ”) and ^ f(~ v z ) (“ ”) are the instant electron velocity distributions sampled from the macro-particles at that time step. The EVDFs are calculated from the macro-particles contained in one cell. The Maxwellian velocity distribution fits using the temperatures and velocities from the macro- particles at each location, ^ f M (~ v x ) (“4”) and ^ f M (~ v z ) (“”), are also plotted for comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 LIST OF TABLES xvii List of Tables 1.1 Typical plasma environmental parameters at LEO, GEO and in SW at 1 AU. 3 1.2 Typical parameters of ion thruster’s plasma plume. . . . . . . . . . . . . . . 5 2.1 Comparison of the performances among different sparse linear solvers. . . . . 38 5.1 Ion-to-electron mass ratio, plasma plume speed, and electron-to-ion plasma frequency ratio at emission surface . . . . . . . . . . . . . . . . . . . . . . . 74 5.2 Best fitting parameters to ion velocity increase in the x-direction Eq. (5.4a). 84 5.3 Best fitting parameters to ion velocity increase in the y-direction Eq. (5.4b). 84 6.1 Ion-to-electron mass ratio, plasma plume speed, and electron-to-ion plasma frequency ratio at emission surface . . . . . . . . . . . . . . . . . . . . . . . 92 LIST OF SYMBOLS xviii List of Symbols mfp Mean free path . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 L c Characteristic length . . . . . . . . . . . . . . . . . . . . . . . . . 1 v te Electron thermal velocity . . . . . . . . . . . . . . . . . . . . . . 1 v ti Ion thermal velocity . . . . . . . . . . . . . . . . . . . . . . . . . 1 v bulk Bulk velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 v d Drifting velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 v ave Average velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 k b Boltzmann constant . . . . . . . . . . . . . . . . . . . . . . . . . 1 T e Electron temperature . . . . . . . . . . . . . . . . . . . . . . . . 1 T i Ion temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 m e Electron mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 m i Ion mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 n e Electron number density . . . . . . . . . . . . . . . . . . . . . . . 2 n i Ion number density . . . . . . . . . . . . . . . . . . . . . . . . . 2 D Debye length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 0 Vacuum permittivity . . . . . . . . . . . . . . . . . . . . . . . . . 2 e Elementary charge . . . . . . . . . . . . . . . . . . . . . . . . . . 2 c;ei Mean collision frequency of electron-ion Coulomb collisions . . . 2 LIST OF SYMBOLS xix mfp;ei Mean free path of electron-ion Coulomb collisions . . . . . . . . . 3 ! p Plasma frequency . . . . . . . . . . . . . . . . . . . . . . . . . . 3 C s Ion acoustic velocity . . . . . . . . . . . . . . . . . . . . . . . . . 3 f(r;v) Particle’s phase-space distribution function . . . . . . . . . . . . 5 r = (x;y;z) position vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 v = (v x ;v y ;v z ) velocity vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 f M (r 0 ;v) Maxwellian velocity distribution function . . . . . . . . . . . . . 5 F Force vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Kn Knudsen number . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Electric potential . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Self-similar parameter in 1D collisionless plasma expansion . . . . 8 Polytropic coefficient . . . . . . . . . . . . . . . . . . . . . . . . . 10 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 E(r;t) Electric field vector . . . . . . . . . . . . . . . . . . . . . . . . . 24 B(r;t) Magnetic field vector . . . . . . . . . . . . . . . . . . . . . . . . . 24 0 Vacuum permeability . . . . . . . . . . . . . . . . . . . . . . . . 25 J Current density vector . . . . . . . . . . . . . . . . . . . . . . . . 25 M Mach number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Non-equilibrium parameter . . . . . . . . . . . . . . . . . . . . . 53 d sh Sheath thickness . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 1 Chapter 1: Introduction 1.1 Overview A plasma is an ionized gas consisting of approximately equal amount of positively charged ions and negatively charged electrons. Due to the long range Coulomb interactions among charged particles, the characteristics of plasmas are significantly different from those of solids, liquids, and neutral gases so that plasmas are considered as the “fourth fundamental state of matter” in the universe. On the one hand, according to the magnitude of the mean free path ( mfp ) relative to the characteristic length (L c ), a plasma can be classified into collisional or collisionless. A collisionless plasma is a type of plasma with the property of mfp L c such that the collisions between individual particles can be neglected in the entire flow region of interest. On the other hand, the thermal energy, characterized by v ti and v te , and the bulk kinetic energy, characterized by v bulk , are also very important parameters to describe a plasma flow. Here, v te = p k b T e =m e and v ti = p k b T i =m i represent the electron and ion thermal velocity, respectively, and v bulk indicates the bulk or the average velocity of a plasma as a fluid. Sometimes people use v d , the drifting velocity, or v ave , the average velocity, to express the same physical meaning as v bulk . In this dissertation, these terms are used interchangeably. In the definition of v te and v ti , k b denotes the Boltzmann constant, T e and T i denote the electron and ion temperature, respectively; m e and m i represent the electron and ion mass, respectively. If the relative magnitudes of the thermal and bulk kinetic energies are considered, a plasma may be categorized as a cold, thermal or mesothermal plasma. In particular, a mesothermal plasma is identified as the flow having a characteristics of v ti v bulk v te . The studies of collisionless, mesothermal plasmas are mainly inspired by their wide applications in the space environment and space engineering. The spacecraft(S/C)-plasma interaction and the electric propulsion are two typical examples in which the behavoirs of the collisionless, mesothermal plasma flows play a critical role. The focus of this dissertation 1.1 Overview 2 Figure 1.1: Artist’s rendition of the plasma iteracting with a space station. is to address the most significant but not well discussed scientific problems in the field of astronautical engineering involving the collisionless, mesothermal plasma flows. 1.1.1 Significance of Collisionless, Mesothermal Plasma Flows The space is not empty. S/C operates in and interacts with the space environment. The properties of space environment where S/C operates are completely different from those near the ground. One of the key features seen in the space environment is the ubiquitous composition of ionized gases. As a result, much of the emphasis has been put on the interactions between the plasma environment and S/C, seeing Fig. 1.1 for an artist’s rendition of the plasma iteracting with a space station. Table 1.1 lists several typical plasma environments at the low Earth orbit (LEO), the geostationary Earth orbit (GEO) and in the solar wind (SW) at1 AU. The data for LEO and GEO are from [3, 4], and for SW at 1 AU are from [5]. The relevant plasma parameters shown in the table are given as follows: n is the plasma number density, and thereby n e and n i represent the electron and ion number density, respectively. D = p 0 k b T e =n 2 e is the Debye length, where 0 is the vacuum permittivity and e represents the elementary charge; c;ei = 1.1 Overview 3 4 p 2 3 e 2 4 0 2 n i ln e m 1=2 e (k b Te) 3=2 is the mean collision frequency of electron-ion Coulomb collisions, where ln e is the Coulomb logarithm; mfp;ei = p 2v te = c;ei is the mean free path of electron-ion Coulomb collisions; ! p = p ne 2 = 0 m is the plasma frequency; C s = p k b (ZT e +T i )=m i is the ion acoustic velocity, where Z is the ion charge number. It needs to be pointed out that for the calculations involving the ion, the properties of the species whose concentration is generally higher (listed in the row of ! pi ) are used. The drifting velocity v d for LEO and GEO at each altitude is estimated based on the orbital velocity of a circular orbit. The dimensions of S/C on orbit vary from 1 to 100 m. This indicates that the plasma interacting with S/C on orbit can be readily treated as collisionless. The mesothermal condition (v ti v d v te ) is valid for almost all LEO flights and the SW plasmas, but not for the GEO, as seen in Table 1.1. However, we must note that the LEO flights are dominant in space missions,and the celestial bodies people are mostly interested in, i.e. the Moon and asteroids, interact with SW. Therefore, the collisionless, mesothermal flow condition dominates the plasma environment which the current space missions are playing with. 100 400 km 500 600 km 1000 km 40000 km Solar Wind n e (cm 3 ) 10 5 10 6 10 6 10 4 10 10 2 10 T e (eV) 0.10.3 0.3 0.20.3 18 12 T i (eV) 0.1 0.1 0.2 0.44 3.5 n(O + )=n e 0.951 0.31 10 2 10 1 10 2 0 n(H + )=n e 0.020.1 0.10.6 0.51 0.91 0.96 n(He + )=n e 10 3 10 1 10 3 10 1 0.1 0.01 0.04 D (m) ' 1 10 2 4 10 3 ' 4 10 2 ' 2:10 8.14 mfp;ei (m) 1:15 10 4 1:24 10 3 1:08 10 5 5:38 10 9 10 11 c;ei (Hz) 28.5 264.0 3.05 3:2 10 4 10 5 ! pe (Hz) ' 1:78 10 7 5:64 10 7 5:64 10 6 ' 5:64 10 5 10 5 ! pi (Hz) 1:04 10 5 (O + ) 3:29 10 5 (O + ) 1:32 10 5 (H + ) 1:32 10 4 (H + ) 4000 (H + ) C s (km/s) 1.55 1.55 6.92 33.90 38.53 v te (km/s) 229.71 229.71 229.71 1186.19 1452.78 v ti (km/s) 0.77 0.77 4.38 19.57 18.31 v d (km/s) 7:67 7:84 7:56 7:61 7.35 2.93 468 collisionless? Yes Yes Yes Yes Yes mesothermal? Yes Yes Yes No Yes Table 1.1: Typical plasma environmental parameters at LEO, GEO and in SW at 1 AU. 1.1 Overview 4 Figure 1.2: Artist’s concept of a solar electric propulsion system. Electric propulsion (EP) has played a key role in many space applications [6], seeing Fig. 1.2 for an artist’s concept of a electric propulsion system on a deep space explorer. The emission and the subsequent expansion of a plume from an EP thruster is important in the interaction of exhaust particles with S/C. Great effort has been made into this fundamental EP problem [7]. Table 1.2 shows the plasma plume parameters from two ion thrusters, one is developed at the Laboratory of Advanced Plasma Dynamics (LAPD) of the University of Southern California (USC) [8], and the other is the NSTAR ion thruster on Deep Space 1 (DS1, the first NASA deep space mission) [9–11]. The plume paramters shown are able to represent the typical plasma plume condition from an ion engine. As exhibited in Table 1.2, the plasma plume from an EP thruster definitely belongs to the collisionless, mesothermal category. Clearly, the collisionless, mesothermal plasma can be found in broad applications of space engineering. It is the significance of the physics and the potential applications involving the collisionless, mesothermal plasma flow that stimulates people to make continous effort to improve the understanding of this subject with more and more advanced methodologies. 1.1 Overview 5 1.1.2 Necessity of Kinetic Investigations In general, two models or theories, the fluid and the kinetic theories, are practically used to describe the dynamics of a plasma flow. In the fluid model, the plasma is regarded as continuummedium, andtherebythemotionofindividualparticleisneglected. Thefluidmodel is based on the approach in fluid mechanics, which basically applies the conservation laws, i.e. of mass, momentum and energy, to the fluid element (see Chapter 3 of Chen’s textbook [12]). The kinetic model, however, considers the evolution of the particle’s distribution function f(r;v) in the phase space, where r = (x;y;z) is the position vector and v = (v x ;v y ;v z ) is the velocity vector. The most important distribution function is the Maxwellian velocity distribution, f M (r 0 ;v) = n (2k b T=m) 3=2 exp m(vv d ) 2 2k b T (1.1) where n, T and v d represent the number density, temperature and drifting velocity in at the position r =r 0 , respectively, and (vv d ) 2 = [(v x v d;x ) 2 + (v y v d;y ) 2 + (v z v d;z ) 2 ] in three-dimensional velocity space. The Maxwellian represents the state of the particles in USC Micro Ion Thruster NSTAR Ion Thruster on Deep Space 1 propellant argon xenon diameter(cm) 4 30 n e (cm 3 ) 2:97 26:7 10 8 1:54 3:22 10 9 (n e 'n i ) n i (cm 3 ) 1:23 11:1 10 8 1:54 3:22 10 9 T e (eV) 2.5 2 T i (eV) 2.5 (T i 'T e ) 2 (T i 'T e ) D (m) 2:3 6:8 10 4 1:9 2:7 10 4 mfp;ei (m) 82:1 687:2 18:7 38:1 c;ei (Hz) 1:38 11:6 10 3 2:23 4:54 10 4 ! pe (Hz) 9:72 29:2 10 8 2:21 3:20 10 9 ! pi (Hz) 2:31 6:94 10 6 4:51 6:52 10 6 C s (km/s) 3.46 1.71 v te (km/s) 663.10 593.10 v ti (km/s) 2.45 1.21 v d (km/s) 72.68 29:8 38:7 collisionless? Yes Yes mesothermal? Yes Yes Table 1.2: Typical parameters of ion thruster’s plasma plume. 1.1 Overview 6 phase space when they are in thermal equilibrium. The governing equation for the kinetic model of plasma is the Boltzmann equation, @f @t +v @f @r + F m @f @v = @f @t c (1.2) where F is the external force acting on the particles, m is the particles’ mass, and (@f=@t) c is the change rate of f due to inter-particle collisions. It can be seen from Eq. (1.2) that the kinetic theory indeed describes the conservation of f in phase space (see Chapter 7 of Chen’s textbook [12]). In principle, the fluid model is applicable only when there are sufficient collisions between particles to keep them in thermal equilibrium or near-equilibrium. The kinetic theory, however, does not rely on the satisfaction of thermal equilibrium and thus is a more general method. In fact, the governing equations used in the fluid model can be derived from the Boltzmann equation, Eq. (1.2), in the continuum limit, seeing Chapter 7 of Chen’s textbook [12] and the monograph by Chapman and Cowling [13]. The dimensionless parameter Knudsen number, Kn= mfp =L c , is used to determine whether the fluid theory or the continuum assumption is valid for modelling a gas flow. When Kn 1, collisions between particles are frequent and thus the fluid model is capable. It is surperising to see that the fluid model works well for plasmas in some places where collisions between particles are not so frequent. There are various reasons for this. For example, Irving Langmuir first reported the observation of Maxwellian distribution for electrons under the condition that the collisions were too few to allow the existence of such a distribution [14]. This is the famous “Langmuir’s Paradox” named after Gabor et al. [15] Another situation that helps the validation of fluid model occurs in the magnetized plasmas, where the magnetic field forces the particles to gyrate, thereby limiting particles’ free-streaming motion. The gyromotion perpendicular to the magnetic field plays a role somewhat similar to the collisions between particles, resulting in a good approxmation of the fluid theory for particles’ motion perpendicular to the magnetic field in an essentially collisionless plasma. 1.2 Literature Review 7 Despite the unusual equilibration observed in the plasma with infrequent collisions, the kinetic model is still necessary or more favorable for the study of collisionless, mesothermal plasma flows in this dissertation. First, as exhibited in Tables 1.1 and 1.2, mfp is at least several orders of magnitude larger than L c in the cases we are interested in. This implies that the plasmas are probably highly nonequilibrium such that the unusual equilibration mechanisms mentioned above are far from enough to maintain the validation of fluid model, and thus the kinetic treatment is required. Second, the fluid model, especially for electrons, is widelyadoptedinstudiesoftheexpansionofEPplumeandS/C-plasmainteraction, sometimes without any justifications. It would be of great significance to carry out investigations using both the fluid and kinetic models under the same flow conditions, and then to answer whether or to what extent the fluid model is applicable by direct comparisons of the results gained from the two different models. Finally, in terms of plasma physics, there exist interesting kinetic phenomena which are not able to be uncovered or explained by the fluid theory. 1.2 Literature Review The expansion is a signature phenomenon in a collisionless, mesothermal plasma flow. The expansionofacollisionlessplasmaintoavacuumoralowerdensityplasmahasbeenintensively studied, especiallyfortheone-dimensional(1D)caseswheretheanalyticalsolutionsareusually able to be derived. This section will first review the milestone investigations on the general physics of a collisionless plasma expansion. The applications of collisionless, mesothermal plasma flows in space engineering will be mostly found in plasma-S/C interactions and plasma plumes from EP thrusters. 1.2.1 Expansion of a Collisionless Semi-infinite Plasma The expansion of a collisionless plasma into a vacuum was first studied theoretically by Gurevich et al. [16] The plasma was considered semi-infinite, initially occupying half-space x< 0, and stationary. In the pioneering work of Gurevich et al. [16], the ion was treated 1.2 Literature Review 8 kinetically with the 1D Vlasov equation (Eq. (1.3)) and the electron was considered as fluid via the well-known Boltzmann relation (Eq. (1.4)). @f i @t +v @f i @x e m i @ @x @f i @v = 0 (1.3) where f i denotes the ion velocity distribution function, is the electric potential. n e =n 0 exp e ( 0 ) k b T e0 (1.4) where the subscript 0 appeared inn e ,T e and indicates the the values at the reference point. The Boltzmann relation in Eq. (1.4) is valid when the electron is fluid (in thermal equilibrium or near-equilibrium) and isothermal. By further assuming quasi-neutrality of the plasma (n e =Zn i ) and introducing the self-similar parameter =x=C s t (Ref. [16] used = p 2=) whereC s = p Zk b T e0 =m i is the ion acoustic velocity withT i = 0, a set of self-similar solutions was obtained. The calculation in [16] exhibits that the ion temperature drops exponentially with , T i ()T i0 exp(2) (1.5) where T i0 is the ion temperature of the unperturbed plasma. Eq. (1.5) shows that the ion temperature descreases very rapidly (exponentially with ), implying that the ion thermal motion can be neglected in most part of the expanding region. If T i 0, the ion can be described by the fluid model without loss of too much accuracy, yielding @n i @t + @ @x (n i v i ) = 0 (1.6) @v i @t +v i @v i @x = C 2 s n i @n i @x (1.7) Eqs. (1.6) and Eq. (1.7) are very similar to the governing equations for a neutral gas in 1D. Eq. (1.6) is nothing but the continuity equation, and Eq. (1.7) is the momentum equation 1.2 Literature Review 9 with the right-hand-side accounting for the electrostatic force. With the help of the paramter , an elegant self-similar solution to Eqs. (1.6) and (1.7) can be obtained, v i =C s ( + 1) (1.8) n i =n e =Z =n 0 exp [( + 1)] (1.9) Plugging Eq. (1.9) into Eq. (1.4) gives e ( 0 ) =k b T e0 ( + 1) (1.10) The analytical self-similar solution given by Eqs. (1.9)-(1.10) for the cold ion model was also obtained by Allen and Andrew [17]. Most recently, Huang et al. obtained a more general self-similar solution by allowing the unperturbed electrons to have arbitrary velocity distribution functions but still keeping the isothermal and quasineutral assumptions [18]. The self-similar solution based on the assumptions of quasi-neutrality and isothermal electron in [16–18] predicts no ion front and unlimited ion acceleration in the expansion region. Such phenomena contradict the observation, and are mainly attributed to the assumptions used to get the self-similar solution. The quasi-neutrality condition is invalid near the ion front, where the local Debye length D is comparable to the density scale length C s t [19]. Dropping the assumption of quasi- neutrality, one needs to take the Poisson’s equation into account, 0 @ 2 @x 2 =e (Zn i n e ) (1.11) which has to be solved numerically. Crow et al. did one of the earliest such calculations using a Lagrangian method with satisfactory accuracy[20]. The numerical results of Crow et al. [20] demonstrates the validity of the self-similar solution in the quasi-neutral region. The formation of an ion front was also revealed. However, the acceleration of the ion front was still indefinite, implying the breakdown of the Boltzmann relation (Eq. (1.4)) at some 1.2 Literature Review 10 point. The identical problem was revisited by Mora recently with a more accurate Lagrangian scheme [2]. A time-dependent expression, matching the numerical result precisely, for the electric field at the ion front was obtained. E front ' 2E 0 = 2e N +! 2 pi t 2 1=2 (1.12) where E 0 = (n e0 k b T e = 0 ) 1=2 and e N in the denominator denotes the constant 2.71828 ..., the base of natural logarithm. Based on Eq. (1.12), other parameters describing the ion front structure, such as v i;front ,x i;front ,n i;front , andn e;front , were derived analytically. Espeically, the ion front speed is given by v i;front ' 2C s ln + p 2 + 1 (1.13) where =! pi t= p 2e N . It is evident that the acceleration of ion is indefinite. A more realistic model is to drop the isothermal assumption and to consider the electron kinetics. In general, the collisionless nature requires a kinetic treatment on the electron behaviors. Denavit performed one of the earliest fully kinetic simulations on this problem [21]. In [21], both ions and electrons were treated as particles with a periodic smoothing method described in [22]. Note that, in [21], the ion still had temperature much less than the electron, and thus may be able to be described by the fluid model properly. The most important finding by Denavit’s fully kinetic simulation is that, rather than an isothermal relation, the electron exhibits a polytropic relation T e n 1 e = const or p e n e = const (1.14) where is the polytropic coefficient. Denavit’s results indicates that = 1 should be proportional to (Zm e =m i ) 1=2 in the expansion region. Futhermore, when the cooling effect of the electron during the expansion is considered, the maximum ion velocity is also 1.2 Literature Review 11 bounded. Mora and Pellat[23] partly considered the electron kinetics by working with the Vlasov equation for electrons, @f e @t +v @f e @x + e m e @ @x @f e @v = 0 (1.15) where f e denotes the electron velocity distribution function. The results in [23] shows that if @f e =@t is negligible as the electrons can be considered to be in equilibrium with the slowly varying electric field and the unperturbed electrons follow a Maxwellian distribution, then the Boltzmann relation relation Eq. (1.4) and thereby the self-similar solutions Eqs. (1.9) - (1.10) can be reproduced. They only represent the zeroth-order solution. If f e allows to vary with time, correction terms need to be added to the zeroth-order self similar solutions. By doing the first-order expansion and assuming the electric potential differs from the self-similar one (1.10) by a correction of order 1=R;R = p m i =Zm e , Mora and Pellat get v i =C s + 1 1 2 d d (1.16) n i =n e =Z =n 0 exp ( + 1) 1 2 (1.17) e ( 0 ) =k b T e0 + 1 1 2 (1.18) where ' (2 p 2=3R) 3=2 when 1R 2 . The most profound outcome from Mora and Pellat’s results is that the first order solution restores the correct electron mean velocity and the energy transfer from the electrons to the ions. v e 'C s ( + 1) (1.19) q e 'n i m i C 3 s (1.20) wherev e is the electron mean velocity andq e is the heat flow. These results show the evidence of electron cooling in the expansion of a collisionless semi-infinite plasma. Sack and Schamel carried out extensive numerical simulations on the same problem [24, 25]. Both ions and electrons were described by the fluid model, but the electron cooling was considered by 1.2 Literature Review 12 incorporating the polytropic relation with the polytropic coefficient varying between 1 and 2. Their simulations show that at the ion front, a spikelike structure is developed and quickly grown up until the numerical solution eventually collapses. These results actually imply that the hydrodynamic model may not be able to take into account the full physics of the collisionless electrons and thus fails to work. 1.2.2 Expansion of a Collisionless Finite Size Plasma Though the electron cooling exists in the expansion of a collisionless plasma, this cooling effect is relatively weak in the semi-infinite situation. As shown by Denavit [21], / (Zm e =m i ) 1=2 is a small amount in general. Therefore, the self-similar solutions derived based on the assumption of electrons’ isothermality reviewed above work fairly well for most of the semi- infinite plasma situations. A physical explanation is that the available amount of energy from the electrons in a semi-infinite plasma can be viewed as infinite. Hence, it is able to drive the expansion and, in the mean time, to keep the electron temperature constant in the expansion region. However, for a finite size plasma, the situation will be completely different because the available amount of electron energy is finite if no external energy source exists. Applying the so-called rescaling methods to the Vlasov-Poisson system, Manfredi et al. studied the expansion of a finite size collisionless plasma into vacuum analytically and numerically [26]. The ion to electron mass ratio was chosen to be relatively small, varying from 2 to 10, to accelerate the process of reaching the asymptotic solution. Their results show that the self-similar solution characterized by v = x=t can be achieved at long time scales. More importantly, the local electron temperature is found to decrease with time as T e /t 2 and, instead of the more commonly seen form pn = const, the following form of polytropic relation should exist for both the ion and electron species. D Dt pn = @ @t +v @ @x pn = 0 (1.21) 1.2 Literature Review 13 where p and n are the partial pressure and number density of each single specie, respectively. For the 1D planar geometry, Manfredi et al. give the polytropic coefficient = 3, indicating an adiabatic cooling process [26]. These results evidently invalidates the assumption of electrons’ isothermality which is almost universally used in the study of semi-infinite plasma expansions. By assuming the quasi-neutrality, cold ions and adiabatic electrons, Baitin and Kuzanyan obtained a self-similar solution to the expansion of a collisionless plasma bunch [27]. Their results shows that the electron temperature changes with time, but the initial Maxwellian distribution function is conserved during the expansion even without collisions. Moreover, the self-similar solution predicts the quadratic dependence of the electric potential on the coordinates /x 2 , in comparison with the linear relation between andx in the isothermal cases. Inspired by the work of Baitin and Kuzanyan [27], Dorozhkina and Semenov [28] derived the exact self-similar solution to the Vlasov equations for the expansion of collisionless plasma bunch only with the quasi-neutral assumption. The quadratic dependence of the electric potential on the positions revealed in [27] is also given in this study. In addition, the electron temperature is shown to inversely proportional to t 2 . Furthermore, in the limit of m e m i , the exact self-similar solution is able to reproduce the adiabatic approximation used in [27]. However, the solution found in [28] is limited to the cases where the ions and electrons in the plasmas need to have similar distribution functions. The self-similarity of the solution also restricts the spatial disctribution of the initial plasma bunch to be the Gaussian one. Kovalev and Bychenkov applied the renormalization-group method to get analytical solutions of the 3D adiabatic expansion of a plasma bunch described by the Vlasov equations [29]. The solutions obtained in [29] is able to deal with arbitary initial velocity distributions for electrons and ions. Two particular examples, the expansion with two temperature electrons and the expansion with two ion species, were studied by the newly derived exact solutions. 1.2 Literature Review 14 The self-similar solutions presented in [27–29] are valid only when the quasineutral condition holds. Mora carried out a hydrodynamic simulation on the expansion of a thin-foil plasma with initially a sharp boundary [30]. Both the electron cooling and the charge separation effects were taken into account, but the electron temperature was assumed to be homogeneous in space. The results show that the electron temperature, in the long time regime, follows T e / t 2 for the nonrelativistic (low energy) limit, while T e / t 1 for the ultrarelativistic (high energy) limit. Note that the T e /t 2 relation for the nonrelativistic electron is also given in [26, 28]. The ion front can be observed because the charge separation effect is considered by this model. Different from the indefinite acceleration of the ion front in the isothermal case, only a finite value of increase in the ion velocity can happen. For the nonrelativistic electrons, a fitting formula of the final ion velocity is given as follows: v i;nal ' 2C s0 ln (0:32L= D0 + 4:2) (1.22) where L is the initial length of the thin-foil, and C s0 and D0 represent the ion acoustic velocity and the Debye length calculated using the initial plasma conditions, respectively. Mora further applied this hydrodynamic model to investigate the expansion of a plasma with a Gaussian distribution of density in space [31], because not only the velocity distribution functions but also the spatial distribution of the plasma can affect the expansion process of a finite size plasma (seeing [28]). Both the isothermal and adiabatic electron expansion were studied. One obvious difference of the expansion of Gaussian plasma from that of steplike plasma is that the ion wave-breaking occurs in the Gaussian case, no matter the electron is isothermal or adiabatic. The self-similarity of the plasma properties is preserved behind the ion front. The finial ion velocity of the Gaussian plasma shows a much weaker dependence on the initial size of the plasma, and is always lower than the steplike plasma. A fitting formula for the final velocity of the Gaussian plasma is given by Eq. (1.23), in comparison with Eq. (1.22). v i;nal ' 0:77C s0 ln (6:3L= D0 + 60) (1.23) 1.2 Literature Review 15 The effects of a finite gradient for the initial ion density profile on the expansions of both semi-infinite and finite size plasmas were also investigated [32]. The study reveals the finite gradient of ion density causes an ion wave-breaking, no matter the plasma is initially semi-infinite or finite size. This configuration has less influence on the subsequent expansion of the initially semi-infinite plasma in comparison with the sharp boundary situation. For the finite size plasma, the ion is accelerated less efficiently when its initial density scale is slightly larger than the total initial size of plasma. Though it provides great insights to the expansion of finite size plasma, the main drawback of the hydrodynamic model used in [30–32] is the assumptions of homogeneous temperature in space as well as the equilibrium state of the electron. To overcome this drawback, one has to consider the electron kinetics. Mora and Grismayer [33] treated the electron kinetically based on an adiabatic invariant [34], while still handled the ion with the hydrodynamic model. It is found that the electron temperature is strongly inhomogeneous and the expansion drives the electron velocity distribution function towards a top-hat shape. Though the electron thermal velocity v te decreases with time as expected, the ion acoustic velocity increases with time in the inner part of plasma at the early stage of expansion. Self-similar solutions to the expansions of finite size non-quasineutral plasmas were derived by Murakami and Basko [35]. However, the solutions are limited to the special polytropic condition of T e (t)/ [n e (t)] 12= , where = 1, 2, and 3 corresponds to the planar, cylindrical, and spherical symmetrical geometries, respectively. 1.2.3 Collisionless, Mesothermal Plasma Flow over an Object The problem of a collisionless, mesothermal plasma flow over an object can be seen in many space-relevant problems. For instance, the solar wind plasma interactions with the Moon or an asteroid at about 1 A.U. can be characterized by the collisionless, mesothermal plasma flow model. Another example is the plasma environment seen by a spacecraft in LEO, which is indeed a collisionless, mesothermal plasma flow problem as well. This problem has been receiving continous attention from the whole community focused on space explorations. 1.2 Literature Review 16 Simulations have been playing a more and more important role in the study of plasma–S/C interactions. Since late 1980’s, the particle-in-cell (PIC) method has been widely applied to investigate the plasma–S/C interactions. Early two- or three-dimensional PIC simulations on the collisionless, mesothermal plasma flow over an object were usually based on the hybrid electron-fluid/ion-particle model. The plasma dynamics in the near wake of a highly negatively biased body immersed in the collisionless, mesothermal plasma environment was investigated by Morgan et al. with both a laboratory experiment and a 2D hybrid PIC simulation [36]. The simulation results were directly compared to the experimental measurement. However, the simulation predicted ion and electron densities in the wake were higher than those obtained from the experiment. A 2D hybrid PIC simulation was carried out to investigate the wake potential and density of a collisionless, mesothermal plasma perturbed by a negatively charged cylinder by Coggiola and Soubeyran [37]. The cylinder was biased to a moderate negative voltage to represent the floating potential. The simulations on both large and small cylinder to Debye length ratios were performed. The ionic emptying phenomenon was observed when the cylinder to Debye length ratio was small and the Mach number was low. Theoretical and 2D hybird PIC modeling study of a collisionless, mesothermal plasma flow over large high-voltage idealized space platforms was reported by Wang and Hastings [1, 38]. The most profount finding in their simulations is abundant ion-rich sheath structures embedded in the quasineutral background wake under different flow conditions. A 2D hybrid PIC simulation was conducted to study a collsionless, mesothermal plasma flow over a cylinder by Guio and Pécseli [39]. Unstable vortices in ion phase-space were observed in the wake of the cylinder. Many distinct macroscopic flow field properties, such as irregular filamented density depletions, and strong velocity shear, are attributed these phase-space structures. Though the hybrid PIC modeling is powerful, the assumption under which this method is applicable does not always hold for the collisionless, mesothermal plasma flow problems. For instance, in situ measurements [40, 41] as well as laboratory experiments [42–45]have reported the abnormal electron temperature enhancement in the wake of an object in the collisionless, mesothermal plasma environment. The fully kinetic approaches are required to 1.2 Literature Review 17 explain such abnormal electron behaviors [46, 47]. However, these early fully kinetic studies are limited idealized 1D model. It is only very recently that the fully kinetic simulations on the collisionless, mesothermal plasma flow of over a more realistic configuration beyond 1D are possible [48–50]. 1.2.4 Electric Propulsion Plasma Plume An electric propulsion (EP) system is a type of device that converts the electric energy to the kinetic energy of the propellant and thereby propels an S/C [51–54]. The main advantage of the electric rockets over the conventional chemical rockets is their highly efficient utilization of the propellant mass. Consequently, EP is almost the only choice for deep space missions [9, 55], as well as the workhorse for many other space applications such as station keeping, orbit raising and orbit insertion [6]. Among the many important scientific problems, e.g. plasma creation, ion acceleration, beam neutralization, etc., involved with EP, we are particularly interested in the neutralization of ion beam and the subsequent expansion of plasma plume from an EP thruster in this dissertation. TheplasmaplumefromanEPthrustermayinteractwiththeS/Candcausecontamination or damage to the facilities on the S/C. Most of these modeling studies of EP thruster plume have focused on the physics driven by the ion dynamics. The electron, due to its small mass, is usually considered to already reach a steady state at the ion time scales and can be modeled as a fluid. Accordingly, the hybrid (electron-fluid/ion-particle) PIC [56] has become dominant in the modeling study of plasma plume from an EP thruster. The simplest and most commonly used model is to treat the electron as massless, isothermal fluid. This model almost ignores the electron dynamics completely, and gives the well-known Boltzmann relation (Eq. (1.4)) which links the electron density to the electric potential. Incorporating Eq. (1.4) into the Poisson’s equation yields a nonlinear Poisson’s equation, which is computationally 1.2 Literature Review 18 expensive to solve. If one further assume the quasi-neutrality, namely n e ' n i , then the electric potential is directly avialable = 0 + k b T e0 e ln n e n e0 (1.24) in Eq. (1.24) is then used to compute the the electrostatic force for pushing the ion particles. A great portion of the previous hybrid PIC studies restored to the simple Boltzmann electron fluid model. Wang and his collaborators have carried out extensive hybrid Boltzmann-electron-fluid/ ion-particle PIC modeling studies of the plasma plume flows [57–61]. 3D hybrid PIC with Monte-Carlo collisions (MCC) simulations were carried out by Wang et al. to study the near S/C plasma environement with an ion thruster plume [57, 58]. The electrons were modeled with the Boltzmann relation, and the nonlinear Poisson’s equation was solved self- consistantly to get the electric potential in the hybrid PIC modelings. The charge-exchange ion density, the charge-exchange ion backflow current density and the spatial electric potential distribution are obtained from the simulations. The simulation results concluded that, with low charging S/C, the thruster’s operation has little effects on the ram side plasma environment. The plasma environment of Deep Space 1 (DS1) was studied by the hybrid PIC-MCC simulation [59]. It is found that the dominant plasma environment around the S/C is the charge-exchange species, which can backflow to the S/C with a moderate charging potential through an expansion process similar to that of a mesothermal plasma expansion into a vacuum. The charge-exchange plasma shields out the S/C potential, and thus the sputterd Mo ions from the thruster are weakly influenced by the S/C charging potential. An immersed-finite-element (IFE)-PIC code was developed to perform large scale simulations of ion thruster plume–spacecraft interations [60, 61]. The incorporation of the IFE method into hybrid PIC model mainly aims at accurately handling complex geometries associated with the realistic spacecrafts. 1.2 Literature Review 19 Roy et al. developed a 3D hybrid PIC model to investigate S/C contamination due to ion thrusters [62, 63]. Two ion thruster plumes and the corresponding the interactions with S/C were studied in [62]. It is found that the plume structure is highly asymmetric, and has a potential profile not just a simple superposition of the backflow from two thrusters. In [63], 3D plume backflow is examined by including the realistic S/C configuration. The “corner-effect”, that is the charge-exchange plasma flowed over the front face of the S/C, is able to be identified by the 3D hybrid PIC simulation. The comparsion between 3D and axisymmetric results concludes that the axisymmetric model may be used to give an upper bound the charge-exchange ion densities for the backflow. The hybrid Boltzmann-electron-fluid/ion-particle PIC model is also the workhorse for many other softwares dedicated to thruster plume-spacecraft interaction modeling [64–66]. In particular, a hybrid PIC code with adaptive mesh refinement (AMR) was reported by Korkut et al. recently [67, 68]. The electron model adopted was the simplest Boltzmann electron fluid plus the quasi-neutral model. This model has been applied to investigate the backflows from single and multiple ion thruster plumes [67] and the ion thruster plumes in ground facilities [68]. The limitations of the hybrid PIC model taking advantage of the Boltzmann relation (Eq.(1.4)) are obvious. For instance, Boyd pointed out that, rather than constant, the electron temperature varies from 10 eV at the thruster exit to 1–2 eV in the far field in Hall thruster plumes [7]. Consequently, the fluid electron models that include more detailed physics are needed in order to accurately simulate plasma plumes from EP thrusters. Roy et al. proposed a modified Boltzmann relation that allows the electron temperature T e in Eq. (1.4) to vary spatially [69, 70]. The electron energy equation was solved to obtain T e . This new hybrid PIC model was applied to study the axisymmetric plume and backflow contamination from an ion thruster. VanGilder et al. [71] employed a local Boltzmann relation similar to the one proposed by Roy et al. [69, 70], whereas the spatial profile of T e was obtained from the experimental measurements [72]. Comparing to the simple isothermal electron model, the varying electron 1.2 Literature Review 20 temperature model is found to greatly affect the properties of the near-field plume from an Hall thruster. The simulation results with the improved model agree better with the measured near-field current density and far-field electron number density. Boyd proposed a detailed electron fluid model that closely coupled the plasma potential and electron energy equation[73]. A hybrid PIC simulation with the detailed electron fluid model was carried out to investigate the near-field plume of a Hall thruster. The detailed model overwhelm the simple Boltzmann electron model in terms of accurately predicting the near-field plasma potential and electron temperature. Recently, Professor Ahedo’s group has been developing the two-fluid model that closes the system of equations by using the polytropic relation Eq. (1.14) to study the collisionless plasma expanions in the magnetic nozzle [74, 75] and far downstream of the thruster plume [76]. The two-fluid modeling of collisionless plasma plumes from ion thrusters is capable of reproducing the fully kinetic PIC results by properly tuning the the polytropic coefficients [77]. Martinez-Sanchez et al. self-consistently solved the ambipolar electric field and the distribution functions of both ions and electrons by considering the energy and magnetic moment conservation, and the quasi-neutrality condition [78]. Their solutions shows a finite potential drop in a magnetized plasma expansion due to the cooling effects, a much more reasonable result than the unbounded potential drop predicted by the isothermal model. While the fluid model is acceptable in studying the collisionless plasma plume where the ion-electron interaction is not very strong (far-field region for instance), it is not applicable to the situations where the electron dynamics is important, such as the ion-electron coupling in the beam neutralization process. The electron dynamics during the ion thruster beam neutralization process was studied by Othmer et al. through 3D full-particle PIC simulations [79, 80]. Their simulations have shown that the electron velocity distribution function might be far from the Maxwellian, indicating the highly non-equilibrium feature, in the plasma plume from an ion thruster. 1.3 Motivation and Objectives 21 Wheelock et al. performed 2D and 3D full PIC simulations to investigate the ion-electron mixing process for an ion thruster beam neutralization [81]. The full PIC simulations have demonstrated the achievement of charge neutralization, but not current neutralization. They argued that either propagating electrostatic waves in the beam or Coulomb collisions needs to be included in the PIC method in order to get current neutralization. Full PIC simulations have been conducted to investigate a plasma beam from an ion thruster with the focus on the neutralization process by Wang and his co-workers [82, 83]. The mechanism of ion beam neutralization revealed from their simulations is mainly attributed to the interaction between the trapped electrons and the potential well established between the emission source and the propagating beam front. The same beam neutralization problem was examined by Usui et al. with the emphasis on the transient behavior of electrons [84]. Their full PIC simulations show that the electrons are trapped by the potential well similar to the structure discussed in [83] and move back and forth in the directions both parallel and perpendicular to the beam propagation direction. Two different streams of electrons are observed near the beam front, but they are too weak to induce two-stream instabilities. In addition, the heating of electrons is observed. The heating mechanism is identified as the electrons’ interaction with the potential well, not due to the local beam instabilities. 1.3 Motivation and Objectives The extensive theoretical analysis of the 1D plasma expansion problems seems to provide very limited guidance on the community of space engineering. However, the 2D/3D plasma wake expansion and plasma plume expansion often seen in real space applications are indeed closely related to the 1D semi-infinite and finite size plasma expansion. A plasma wake expansion can be analogous to a 1D semi-infinite plasma expansion into a vacuum in a way that the plasma entering into the expansion region is from an unperturbed source with unlimited amount of particles and energy, whereas a plasma plume expansion is analogous to a 1D 1.3 Motivation and Objectives 22 finite-size plasma expansion into a vacuum in a way that the plasma’s mass and energy are limited. It has been shown from the literature review that the behavior of a finite-size plasma expansion into a vacuum can be quite different from that of a semi-infinite plasma expansion. The 1D theoretical conclusions will be beneficial to the analysis of 2D/3D realistic problems. More importantly, it can be seen from the previous section that the electron fluid modeling dominates the study of collisionless, mesothermal plasma flows related to space-engineering (plasma-S/C interaction and EP plume), mainly due to the computational tractability. However, in view of the collisionless nature, the kinetic approach is necessary for these problems. Verification of when or to what extent the fluid modeling is valid for collisionless, mesothermal plasma flows will be important to the community of space engineering as well as an interesting scientific problem. This disseration aims at advancing the current understanding of the applicability of the fluid modeling for space-engineering relevant collisionless, mesothermal plasma flows. Fully kinetic simulations will be carried out and the results will be compared with those obtained from the fluid modeling. The connection between the macroscopic flow structure and microscopic ion/electron kinetic characteristics will be analyzed in detail. Criteria for the failure of the fluid modeling will be provided by systematic comparison between the fully kinetic and the fluid modeling results. The theoretical models in the literature will be extended to describe the more complex realistic space applications in 2D and 3D. 1.4 Dissertation Overview 23 1.4 Dissertation Overview The rest of this dissertation is outlined as follows: 1. Chapter 2 presents a brief overview of the commonly used numerical methods for studying the collisionless plasma problems. 2. Chapter 3 studies the collisionless, mesothermal flow over a large unbiased conducting plate. 3. Chapter 4 studies the effects of surface charging on the collisionless, mesothermal plasma flow over a large conducting plate. 4. Chapter 5 studies the expansion of a collsionless, pre-neutralized plume from an EP thruster into vacuum. 5. Chapter 6 studies the 3D effects on the expansion of a collsionless, pre-neutralized beam from an EP thruster into vacuum. 6. Chapter 7 carries out a comprehensive study of the beam neutralization and the subsequent plume expansion. 7. Chapter 8 briefly summarizes the work done in this dissertation and highlights the main conclusions. Several potetial works deserving further investigations are also outlined. 24 Chapter 2: Kinetic Simulation Method The collisionless, mesothermal plasma flows essentially happen in two (2D) or higher dimen- sions, and therey are qualitatively different from one-dimensional (1D) problems. The kinetic studies of collisionless plasma problems in dimensions higher than 1D mostly rely on numer- ical simulations. This chapter will first give a brief overview of commonly used kinetic simulation methods for studying the collisionless plasma problems. Next, the particle-in-cell (PIC) method will be discussed because it is the conventional tool for kinetic simulations of the collisionless mesothermal plasma flows. Finally, several remarks on full particle PIC simulations of the collisionless, mesothermal plasma flows are addressed. 2.1 Introduction As discussed in the last chapter, a kinetic discription is necessary for collisionless plasma dynamics. In a collisionless plasma system, the interactions between charged particles are through the self-consistent collective fields created by themselves and the pair collisions are negligible. Then, the Boltzmann equation [Eq. (1.2)] is specialized as @f j @t +v @f j @r + q j m j (E +vB) @f j @v = 0 (2.1) where f j (r;v;t), q j and m j represent the distribution function, charge and mass of the j th species of the plasma, respectively, and the electric field E(r;t) and the magnetic field B(r;t) are solved from the Maxwell’s equations, rE = 0 (2.2) rB = 0 (2.3) 2.1 Introduction 25 rE = @B @t (2.4) rB = 0 J + 0 @E @t (2.5) Here, 0 and 0 denote the permeability the permittivity of vacuum, respectively. is the charge density and J is the current density, determined as follows, (r;t) = X j Z q j f j dv; J(r;t) = X j Z q j vf j dv Eq. (2.1) is called the Vlasov equation. Eqs. (2.1)-(2.5) are the Vlasov-Maxwell system of equations, which are the complete set of the governing equations for the kinetic description of collisionless plasmas. When self-induced magnetic field B of the plasma can be neglected, the Vlasov-Maxwell system is reduced to the Vlasov-Poisson system, @f j @t +v @f j @r + q j m j E @f j @v = 0 (2.6) E =r (2.7) r 2 == 0 (2.8) where (r;t) represents the electric potential. There are mainly two categories of numerical simulation methods, the grid-based method or the particle-based method, to deal with the kinetic equation [Eq. (2.1) or Eq. (2.6)] for the evolution of plasmas. A grid-based or Eulerian kinetic method solves the Vlasov equation by directly discretizing the mesh in phase space. The grid-based kinetic method takes great advantage of the broad numerical techniques for partitial differential equations (PDEs) because Eq. (2.1) or Eq. (2.6) is just a first order hyperbolic PDE in multi-dimensions. A large portion of the grid-based kinetic methods belongs to the semi-Lagrangian (SL) type of method, which tracks the information propagating along the characteristics in each simulation step. The SL type of method allows one to use relatively larger steps than the conventional Eulerian numerical 2.1 Introduction 26 methods without sacrificing accuracy due to the Lagrangian evolution mechanism. The SL type of method for the Vlasov equation has been under active developments for finite difference [85–91], finite volume [92–95], finite element [96], and hybrid finite element-finite difference [97] methods. The abovementioned SL methods are caste into the so-called the backward SL (BSL) method. In a BSL method, to find a quantity f for a given mesh point in phase space at time t + t, one starts from the point and follow the characteristic curve backward in time until t. Then, the information at the arrival point at time t is used to assessf at timet + t. The arrival point at timet is most likely not located at a mesh point, and thereby interpolation is required. The BSL type of method is the most commonly used Eulerian method for solving hyperbolic PDEs. Another type of kinetic SL method for solving the Vlasov equation is the forward SL (FSL) method first introduced by Denavit [22] and recently extented to the modern fashion [98, 99]. The FSL method is very similar to the particle method, except that the “particle” in the FSL method is an element in phase space with the “mass” equal tof(r;v)drdv. All the BSL and FSL methods discussed above followed Cheng and Knorr [85] and utilized the Strang operator splitting to solve f in ther-space and the v-space separately. However, the Strang splitting is 2nd order accurate in time and can cause problems in long time simulations. High order splitting methods have been proposed to resolve these issues [100, 101]. Another remedy is to employ the unsplitting schemes with the cost of higher program complexity and computational efforts. The unsplitting schemes have been incorporated into the convetional non-SL type of finite-volume [102–104] and finite-element [105] methods, as well as the SL type of finite-element method [106]. A particle-based method uses a few number of particles to statistically represent the distribution function f(r;v), and solves the motions of particles in phase space. The particle-based method can be divide into the particle-mesh (PM), particle-particle (PP), and particle-particle-particle-mesh (P 3 M) [107]. The PM, or most frequently called the particle-in-cell (PIC), method is the least expensive method, but able to retain most of the physics. Consequently, the PIC method is the most widely used particle-based method in studying the plasma physics. In the PIC method, the velocity distribution function (VDF) 2.1 Introduction 27 Figure 2.1: Schematic of the distribution function loosely represented by the particles in a cell. f(v) in a cell is loosely represented by the particles in a statistical way, as shown in Fig. 2.1. The drawback of the particle representation of VDF is that the statistical noise proportional to 1= p N p , where N p is the number of particles, might cause numerical challenges in some applications. Comparing to the particle-based method like PIC, the grid-based methods of solving the Vlasov equation are noise free, and less difficult to incorporate high-order accurate algorithms. The cost of these advantages is the extremely fast increase in computational demands as the problem dimension increases. For instance, a grid method usually requires O(10 2 ) mesh points in one direction in the v-space to properly caputure the distribution function f. If we consider a fully 3D (3D in both the r- and v-space) problem with only 100 mesh points in each direction of the x-space, the total number of mesh points needed to represent f(r;v) in the six-dimensional phase space is then O(10 12 ). Though solving a problem with such a tremendous size is possible on the most powerful supercomputers, the grid-based method is definitely not desired in terms of the economic efficiency, let alone the problem with only 100 mesh points in one dimension is indeed very small in the r-space. Nevertheless, 100 particles per cell are usually enough for a typical PIC simulation. Then, the size of the fully 3D problem that the PIC simulation needs to handle is only O(10 8 ), which is feasible even with a desktop. As a result, the PIC method is adopted as the kinetic simulation tool in this 2.2 Particle-in-Cell Method 28 Figure 2.2: Flow schematic for the ES-PIC scheme. dissertation. Moreover, in all the applications of this dissertation, the effects of magnetic field are negligible so we only focus on the Vlasov-Poission system of Eqs. (2.6)-(2.8). 2.2 Particle-in-Cell Method 2.2.1 Fundamentals The particle-in-cell (PIC) method solving Eqs. (2.6)-(2.8) is referred as the electrostatic particle-in-cell (ES-PIC) method. The flow of the ES-PIC scheme is illustrated in Fig. 2.2. The PIC method defines the particles continuously in both the r-space and v-space, and the fields discretely in ther-space. The particle in a PIC simulation is the so-called macroparticle which represents many real charged particles. The dynamics of particles are solved according to the Newton’s law Eqs.(2.9) - (2.10). d (mv) dt =qE (2.9) dr dt =v (2.10) where the electric field E is obtained according to Eqs. (2.7) and (2.8), and q and m are the charge and mass that one macroparticle carries, respectively. 2.2 Particle-in-Cell Method 29 Figure 2.3: Schematic of force weighting and particle weighting based on the linear weighting scheme in 2D. As shown in Fig. 2.2, a typical ES-PIC loop consists of four steps: (i) force weighting, (ii) particle pushing, (iii) particle weighting, and (iv) solving the electric field. One can find the detailed discussion on the PIC algorithms from the excellent monographs [107, 108] as well as some other review articles [109–111]. We will only briefly describe the fundamentals for each of the four steps relevant to the implementation in this dissertation. (i) Force weighting. The electric field is available only at discrete mesh points. Hence, to get the electrostatic force exerted on a particle, an interpolation or a weighting of the force from the mesh points to the particles’ positions is needed. The linear weighting scheme is a good balance between the accuracy and efficiency, and thus adopted as a standard weighting scheme. Fig. 2.3 shows how to get the force F exerted on a particle within a rectangle cell using a 2D example. The cell size is x =x i+1;j x i;j =x i+1;j+1 x i;j+1 and y =y i;j+1 y i;j =y i+1;j+1 y i+1;j . Supposing the particle’s coordinate is (x p ;y p ), we can compute w x = x p x i;j x ; w y = y p y i;j y (2.11) 2.2 Particle-in-Cell Method 30 Then, the force on the particle is obtained according to F =q [(1w x )(1w y )E i;j +w x (1w y )E i+1;j + (1w x )w y E i;j+1 +w x w y E i+1;j+1 ] (2.12) It is clear that the contribution of theE field at a mesh node to the force is proportional to the fractional area/volume of the subcell away from the node on the diagonal. For instance, in Fig. 2.3, the contribution of E at node (i;j), shown by the green circle, to the particle is proportional to the ratio of the green subcell area to the total cell area. The other node-subcell pairs illustrated by various colors in Fig. 2.3 follow the same rules. (ii) Particle pushing. The motion of a particle in phase space is integrated according to Eqs.(2.9) - (2.10). Numerically, this is done with the leap-frog scheme: v n+ 1 2 =v n 1 2 + q m E n t (2.13) r n+1 =r n +v n+ 1 2 t (2.14) where t is the time step used in the PIC simulation, and the superscript n denotes the simulation performed at time t =n t. Note that, in the leap-frog scheme, the times at which the velocity v and the position r of a particle are available are offset from each other by a half time step 1 2 t. The leap-frog scheme is the 2nd order accurate in time, and numerically stable for the oscillatory motion with a frequency ! as long as t 2=!. (iii) Particle weighting. A particle in a cell should distribute its charge to the cell nodes, so the charge density on the nodes will be available and used to solve the electric field. This process is called particle weighting. The particle weighting should employ the same weighting scheme as the force weighting to avoid self-force on the particles [108]. 2.2 Particle-in-Cell Method 31 As shown in Fig. 2.3, the particle in the cell distributes its charge to the four nodes according to Q i;j =q(1w x )(1w y ) Q i+1;j =qw x (1w y ) Q i;j+1 =q(1w x )w y Q i+1;j+1 =qw x w y (2.15) whereq is the charge carried by the macroparticle (charge of one real particle multiplied by the number of real particles that this macroparticle represents), and Q i;j is the charge deposited on node (i;j). The charge on the mesh nodes is accumulated according to Eq. (2.15) for all the particles. Afterwards, the charge density at node (i;j) is computed as follows i;j = Q i;j xy (2.16) (iv) Solving E field. After the charge density is available, we can numerically solve Eqs. (2.7) and (2.8) for the electric fieldE. Specifically, this dissertation first solves the Poisson’s equation (2.8) by using the finite difference scheme with central differencing, i1;j 2 i;j + i+1;j x 2 + i;j1 2 i;j + i;j+1 y 2 = i;j 0 (2.17) Here, the cell size x or y is assumed to be the same among all the cells. Once is solved, the electric field at a given node (i;j) is calculated based on the central difference scheme as well, E x i;j = i+1;j i1;j 2x ; E y i;j = i;j+1 i;j1 2y (2.18) where E x i;j and E y i;j indicate the x- and y-components of the electric field vector E at node (i;j), respectively. Note that schemes shown by Eqs. (2.17) and (2.18) are both 2nd order accurate. 2.2 Particle-in-Cell Method 32 2.2.2 Full PIC vs. Hybrid PIC According to how the electrons are handled in a numerical model, one can distinguish between the fully kinetic or full particle PIC model, in which both the electrons and ions are treated as particles, and the hybrid PIC model, in which the ions are treated as particles while the electrons are considered to be a fluid. In particular, the charge density needed to solve the Poission’s equation [Eq. (2.17)] can be written as =e X k Z k n i (k)n e ! (2.19) where e is the elementary charge, n e is the electron number density and n i (k) is the ion number density for the k th species whose charge number equals Z k . Both n i and n e are obtained from particles’ information in full PIC simulations while n i is from particles but n e is based on the fluid models in hybrid PIC simulations. Implementing the full PIC model is straightforward by following the four steps presented in the last subsection. However, in order to resolve the electron dynamics, the time steps used in the full PIC simulations are limited to be smaller than the electron time scales. Then, the full PIC simulations usually need to be run for a long time because the dynamics relevant to the collisionless mesothermal plasma of interest happen at the ion time scales. Another issue with performing accurate full PIC simulations on the collisionless mesothermal plasma flows is that a very large simulation domain is normally required to reduce the boundary condition effects, which will be discussed in detail in the next section. Consequently, the full PIC models are much more computationally demanding than the hybrid PIC models. Many previous studies considered that the collisionless mesothermal plasma flows were governed by the ion dynamics, and the electron dynamics could be approximated by simplified fluid models. Hence, the hybrid PIC approach has been widely used. The electron fluid models used in past hybrid PIC simulations may be categorized as globally isothermal fluid models [10, 59, 61, 62, 67, 112], simple local temperature dependent models [69, 71, 113, 114], and detailed thermal energy transfer models [73, 115, 116]. The isothermal electron model 2.2 Particle-in-Cell Method 33 plus the assumption of massless electrons yields the Boltzmann relation Eq. (1.4), which is rewritten here, n e =n 0 exp e ( 0 ) k b T e0 (2.20) where “0” in the subscript indicates the parameters at the reference point, and n e , , T e , k b and e denote the electron density, electric potential, electron temperature, Boltzmann constant and elementary charge, respectively. Eq. (2.20) is incorporated into Eqs. (2.19) and (2.8) to formulate the most widely used hybrid PIC model. In this model, one needs to the non-linear Poisson’s equation (2.21), r 2 =e Zn i n 0 exp e ( 0 ) k b T e0 (2.21) where only one ion species is supposed to exist for simplicity. The non-linear equation (2.21) is solved by using the Newton’s method [108]. It is worth mentioning that solving Eq. (2.21) for is relatively more difficult to be programmed compared with solving the linear Poisson’s equation in full PIC simulations. Moreover, it is generally more time consuming to solve Eq. (2.21) in one PIC iteration. The major benifit of employing a hybrid PIC model is to use a much larger time step, and thus to save the overall computational time compared with performing full PIC simulations. The assumption of isothermal electron used to derive the Boltzmann relation Eq. (2.20) cannot always be true. Some authors worked with the more general polytropic relation of state Eq. (1.14), in which the effects of electron cooling/heating can be taken into account. Incorporating the polytropic relation into the electron fluid equations yields the dependence of electron density on the potential and the polytropic coefficient , n e =n 0 1 + 1 e ( 0 ) k b T e0 1=( 1) ; 6= 1 (2.22) 2.2 Particle-in-Cell Method 34 The corresponding non-linear Poisson’s equation is given by r 2 =e Zn i n 0 1 + 1 e ( 0 ) k b T e0 1=( 1) ! ; 6= 1 (2.23) Inthisdissertation, thehybridPICmodelsbasedonboththeBoltzmannrelationEq. (2.20) and the polytropic relation Eq. (2.22) are implemented and used to compare with the fully kinetic ES-PIC model. The purpose is to assess the applicability of the widely used electron fluid models benchmarked against the high-fidelity full PIC model. 2.2.3 Normalization In numerical practice, normalization of the originial governing equations is usually necessary in order to reduce the round-off error due to the finite precision representation of floating numbers in computers. The normalization scheme used for the full particle PIC is as follows, ~ r = r D0 ; ~ v = v v te0 ; ~ t =t! pe0 ; ~ m i;e = m i;e m e ; ~ n i;e = n i;e n 0 ; ~ = e k b T e0 ; ~ E = E k b T e0 =e D0 (2.24) and the normalization scheme for the hybrid PIC is given by, ~ r = r D0 ; ~ v = v C s0 ; ~ t =t! pi0 ; ~ m i = m i m i = 1; ~ n i;e = n i;e n 0 ; ~ = e k b T e0 ; ~ E = E k b T e0 =e D0 (2.25) In eqs (2.24) and (2.25), the subscript “0” denotes the reference condition, D0 = p 0 k b T e0 =n 0 e 2 is the Debye length, ! pe0 = p n 0 e 2 = 0 m e and ! pi0 = p n 0 e 2 = 0 m i are the electron and the ion plasma frequency, respectively, v te0 = p k b T e0 =m e is the electron thermal velocity, andC s0 = p k b T e0 =m i is the ion acoustic velocity. It is noted that the normalization for the fully kinetic [Eq. (2.24)] and that for the hybrid PIC modeling [Eq. (2.25)] differ in the mass ( ~ m), time ( ~ t) and velocity (~ v). As ! pi0 ! pe0 , the hybrid PIC simulation is performed with a much larger time step than that for the fully kinetic PIC, thereby accelerating the simulation. 2.3 Remarks on Full Particle Simulation 35 The normalized governing equations are the same for both the full PIC and hybrid PIC models, which are written as d~ r d ~ t = ~ v (2.26) d ( ~ m~ v) d ~ t = ~ F (2.27) ~ F = ~ q ~ E (2.28) ~ r 2 ~ = (Z~ n i ~ n e ) (2.29) ~ E = ~ r ~ (2.30) Note that ~ q = Z and1 for ions and electrons, respectively. In this dissertation, all the simulations performed and the results presented are based on the normalized physical parameters Eqs. (2.24)-(2.25) and the dimensionless governing equations (2.26)-(2.30). 2.3 Remarks on Full Particle Simulation 2.3.1 Boundary Conditions One of the key issues affecting the accuracy of simulations is the implementation of boundary conditions. In a PIC simulation, both the particle and the field boundary conditions must be specified. The standard field conditions at a boundary with the Poisson’s solver can be chosen as the Dirichlet condition, the Neumann condition, or the mix of the two. Special attention must be paid to the particle boundary conditions with the collisionless, mesothermal plasma flows. The property ofv ti v bulk v te for a mesothermal flow indicates that the ion is cold and hypersonic, while the electron is thermal and subsonic. In a particle simulation of hypersonic gas flow, one can apply the so-called “stream” condition at the inlet and the “vacuum” condition at the outlet. At the inlet, a Maxwellian distribution can be assumed. This is true for our cases of collisionless mesothermal plasma flows around an object and the plasma beam emitted from an EP thruster. In addition, the macroscopic flow 2.3 Remarks on Full Particle Simulation 36 conditions, such as the density, bulk velocity, and temperature, etc., of the gas are usually known at the inlet. Then, one can use the known inflow macroscopic parameters to easily construct the local Maxwellian, from which the velocities of injected particles are sampled. When a particle arrive at the hypersonic outlet, it is simply removed from the simulation domain. As a result, it is easy to deal with the particle boundary condition for the ions. For the electrons in our applications, we can still consider the local Maxwellian is kept and the macroscopic flow conditions are known at the inlet. Then, the same “stream” boundary condition as that for injecting ions is applied to the electron inflow condition. The diffculty arises from accurately implementing the outlet boundary condition for the electrons. The common strategy is to use a relatively small simulation domain so the downstream outflow boundary is inside the plasma region. Then, a certain amount of particles need to be re- injected into the simulation domain according to the outlet condition because the electrons are thermal and subsonic. This is the standard strategy in particle simulations of neutral gases. Several methods have been proposed to deal with the subsonic outflow boundary conditions for the particle-based methods [117–123]. Regardless of the difficulties in the accurate implementation itself, these subsonic flow conditions for the particle-based methods implicitly assume the gas to be in equilibrium at the boundaries such that appropriate local VDFs can be constructed from a few macroscopic flow parameters. However, as we will see in the following chapters, the local VDFs are very likely to be highly non-equilibrium in the plasma region due to the collisionless nature of the flows we are dealing with. Therefore, we cannot use these methods for the electron outflow condition if we put the outlet boundary inside the plasma region. A practical remedy is to use a very large simulation domain so outlet boundary is still distant from the plasma propagation front when the simulation ends to minimize the boundary condition effects. 2.3.2 Structured Stretching Mesh To correctly capture the physics of plasma, the cell size for a PIC modeling is required to be not bigger than the local Debye length D . D is proportional to 1= p n where n is the 2.3 Remarks on Full Particle Simulation 37 Figure 2.4: Illustration of a stretching mesh being used for the plasma plume simulation in this dissertation. plasma density. Most PIC simulations use the uniform mesh. Such a treatment is acceptable in those problems as the plasma density and thus the local Debye length are relatively uniform. Unfortunately, the plasma density can vary across many orders of magnitude in the problems of this dissertation. Then, it is very wasteful to use uniform meshes because we have to use small cell sizes to resolve the high density region. As discussed in the last subsection, a very large simulation domain is necessary for a full particle PIC simulation to alleviate the effects of boundary conditions. This will make the waste of computational resources much worse if uniform meshes are used. In this dissertation, we adopt the structured stretching mesh to resolve this issue. The advantage of structured stretching mesh is the computational convenience and the cheap usage of memory, compared to other strategies such as the unstructured adaptive mesh refinement (AMR) [67, 68]. For instance, the simulation only needs to storeN x +N y grid points for a 2D simulation if a structured stretching mesh is used. Fig. 2.4 shows an example of the computational mesh used for simulating a 2D plasma plume expansion. 2.3.3 Poisson Solver The discretization of the Poisson’s equation leads to a large sparse linear system to be solved. Though numerous methods and tools exist, solving a large sparse linear system could be 2.3 Remarks on Full Particle Simulation 38 Case 1 Computing time (sec.) Iterations to converge CG 3.8165 3788 PCG(Jacobian) 0.4753 360 PCG(ILU(0)) 0.2333 111 Direct LU 0.0139 (3.3926) N.A. Case 2 Computing time (sec.) Iterations to converge CG 39.4473 3788 PCG(Jacobian) 5.1155 713 PCG(ILU(0)) 2.4852 213 Direct LU 0.2894 (124.96) N.A. Table 2.1: Comparison of the performances among different sparse linear solvers. very computationally expensive. It is interesting to note that people carrying out ES-PIC simulations have been trying to avoid using the direct method to solve the Poisson’s equation. All kinds of iterative solvers, e.g. successive over relaxtion (SOR), diagonalized alternating direction implicit (DADI), (preconditioned-) conjugate gradient (PCG), multi-grid (MG), etc., have been reported to be incorporated into ES-PIC solvers in the literatures. Among these, the modern PCG and MG methods are considered to have very good performance against the others. For a linear system with n equations, a direct method, such as the LU method, scales asO(n 3 ), which is indeed expensive. However, in an ES-PIC simulation, the coefficient matrix of the linear system, once set up, never changes during the entire simulation. Consequently, the expensive LU factorization is done only once at the beginning of the simulation and the factorized lower and upper triangular matrices are stored. Afterwards, only the cheap forward and backward substituitions are needed to solve the system of equations at each PIC loop. For more details about the direct methods for sparse linear systems, one can refer to the great monograph by Davis[124]. A comparison between the PCG solver (with different preconditioners) and the direct method of LU is carried out. Two test cases produced by the IFE-PIC code developed in our group [61] are used. The first case (Case 1) has 10302 IFE nodes and the second one (Case 2) has 66330 IFE nodes. The performances with different methods are listed in Table 2.1. The pure conjugate gradient (CG), PCG with Jacobian preconditioner, PCG with ILU(0) preconditioner, and the direct method of LU are used 2.3 Remarks on Full Particle Simulation 39 for this comparison. The time used for performing the LU factorization is presented in the parentheses. Table 2.1 shows a clear evidence that using the direct LU solver for the ES-PIC modeling could substantially increase the performance against an iterative solver. 40 Chapter 3: Collisionless, Mesothermal Plasma Flow over a Large Unbiased Conducting Plate: Breakdown of the Electron Fluid Approximation in the Wake This chapter presents a comparative study using fluid based analytic solution, hybrid particle- ion fluid-electron particle-in-cell (PIC), and fully kinetic PIC for a collisionless, mesothermal plasma flow over a large conducting plate. The plate is considered to have the same potential as the ambient plasma. The main focus is on the validity of the fluid approximation for electrons in the wake behind the plate. It is found that the plasma wake may be characterized into two regions based on electron characteristics: a fluid electron expansion region and a kinetic electron expansion region. In the fluid electron expansion region, the electrons may be considered as an equilibrium fluid and all three approaches leads to a similar result. In the kinetic electron expansion region, the local electron velocity distributions are strongly non-Maxwellian and the commonly used fluid approximation for electrons breaks down. It is also shown that the error of the hybrid PIC simulation strongly correlates with a sophisticated non-equilibrium parameter taking into account both the electron temperature anisotropy and the deviation of the electron velocity distribution from the equilibrium. 3.1 Introduction The formation of plasma wake behind an object is relevant to many application problems in space physics and spacecraft-plasma interactions, and has been a subject of extensive 3.1 Introduction 41 modeling and simulation studies. Previous modeling and simulation studies range from magnetohydrodyamics (MHD) and electromagnetic hybrid particle-in-cell (PIC) simulations of solar wind plasma flow around planets and large asteroids [125–128], electrostatic (ES) hybrid PIC simulations of plasma flow around spacecraft and spacecraft charging [1, 38, 129, 130], to ES grid-free particle simulations of plasma flow around small asteroids [50]. In these problems, the plasma flow is mesothermal (v ti v 0 v te ) with a Mach numberM=v 0 =C s 1, where v ti , v 0 , v te , and C s are the ion thermal velocity, plasma drifting velocity, electron thermal velocity, and ion acoustic velocity, respectively. For spacecraft and many small asteroids, their size L is smaller than the collision mean-free-path of the plasma as well as in the range of D <L<r ci , where D is the plasma Debye length and r ci is the ion gyro-radius. Hence, in the near wake region, the plasma flow can be considered collisionless and the effects of the magnetic field are not important. This paper considers the plasma wake in a collisionless, unmagnetized, mesothermal plasma. We note that this subject is also closely related to the plasma expansion into vacuum, a topic studied extensively within the context of plasma dynamics [2, 16, 17, 20, 21, 131]. As previous studies of plasma wake/plasma expansion mostly focus on the physics driven by ion dynamics, a common approach has been to consider the electrons as an equilibrium massless fluid. Moreover, in almost all studies of unmagnetized plasma wake or expansion, the electrons were further assumed to be an isothermal, ideal gas. Thus, from the electron momentum equation and state equation 0'm e @v e @t +v e rv e =errp e ; p e =n e k b T e ; (3.1) where m e , n e , u e , p e , and T e denote the electron mass, density, average velocity, pressure, and temperature, respectively, k b the Boltzmann constant, e the elementary charge, and the electric potential, one obtains the commonly used Boltzmann relation between n e and , 3.1 Introduction 42 Eq. (2.20). If one allows the electron cooling in the process, the dynamics of electrons can also be explicitly solved from Eq. (3.1) by introducing the polytropic law T e n 1 e = const; or p e n e = const (3.2) where is the polytropic coefficient. The resultant relation between n e and is given by Eq. (2.22). Recently, the electron fluid model taking into account the cooling effect through the simple polytropic law, Eq. (3.2), has been revisited and applied to study the expansion of collisionless plasma plume [74–76]. It is worth noting that = 1 in Eq. (3.2) corresponds to the isothermal case. A fluid approximation for electrons simplifies the solution of electron dynamics and leads to significant savings in computational time in simulations. However, for a fluid electron model to be valid, the electrons need to be in a near thermodynamic equilibrium state. The primary mechanism for electrons to reach the equilibrium state is through random collisions between particles. The long-range Coulomb interaction between charged particles in a plasma is generally not considered as an effective mechanism for equilibrium [132]. Hence, while a fluid approximation can be justified for electrons in plasma wakes behind planets or large asteroids, its validity has not been resolved for electrons in plasma wakes behind spacecraft or small asteroids where the characteristic length scale is often smaller than the collision mean free path. Indeed, several studies have indicated that a fluid electron model is not adequate for certain plasma wake or plasma expansion problems. For instance, various observations have found anomalous electron temperature enhancement in the wake of ionospheric satellite [133–136] which were speculated due to a non-Maxwellian electron velocity distribution [46]. Mora and Pellat [23] examined the use of Eq. (2.20) for electrons in the derivation of the self-similar solution of collisionless plasma expansion into vacuum and found that the exact electron density differs substantially from that predicted based on Eq. (2.20). 3.2 Formulation and Analytical Solution 43 Figure 3.1: Schematics for (a) plasma wake and (b) simulation setup. This chapter investigates the validity of the fluid approximation for electrons in a plasma wake. We consider a large, unbiased conducting plate of size L p D in a collisionless, unmagnetized, mesothermal plasma flow. The plasma wake is studied using ES hybrid PIC and full particle PIC simulations, as well as a fluid based analytical solution. In the hybrid PIC simulation, the ions are treated as macro-particles and the electron dynamics is simplified by either the Boltzmann relation, Eq. (2.20), or the more general polytropic relation, Eq. (2.20). In the full particle PIC simulation, both ions and electrons are treated as macro-particles and the actual proton to electron mass ratio m i =m e = 1836:15 is used. A direct comparison of these three approaches is carried out to study the effects of electron kinetic characteristics on the plasma wake, and to determine the validity of the fluid approximation for electrons. 3.2 Formulation and Analytical Solution If one considers that the plasma is quasi-neutral in the entire flow field, an analytical solution for the near wake region can be readily obtained by considering both the ions and electrons as fluid [1]. To include the plasma sheath around the plate and ion kinetic characteristics while retaining the assumption that the electrons can be modeled as a fluid, this problem can be simulated using the hybrid particle-ion fluid-electron PIC approach. To further include electron kinetic characteristics, the full particle PIC approach is required. For a collisionless, hypersonic (M 0 1) quasi-neutral plasma flow, the equation of motion for the ions can be written equivalently in the form of the steady state ion momentum equation except in the far wake region where the ions cross. Hence, in the near wake, the ion trajectories are the same as the streamlines. If one further assumes that the electrons can be 3.3 Simulation Setup 44 approximated by an isothermal, perfect gas, it can be shown that the plasma flow follows the same governing equations mathematically as that for an isothermal compressible gas flow. Thus, as shown in Fig. 3.1 (a), the perturbation generated by the plate is analogous to that of the Prandtl-Meyer expansion of a supersonic gas flow around a corner [1]. The edge of the plate generates an expansion fan. The first characteristic line of the expansion fan is the Mach line which separates the undisturbed plasma flow from the expansion region. The Mach line extends from the plate edge outward into ambient plasma at an angle 0 ' sin 1 (1=M 0 ) with respect to the incoming flow direction. The expansion fan is a presheath. Behind the Mach line, both the plasma density and electric potential decrease as the plasma expands into the wake. As the ions go through the expansion, their trajectories are turned into the wake in such a way that the velocity component normal to each expansion characteristic line is always at the ion acoustic speed C s . The plasma density and electric potential behind the Mach line as a function of expansion angle from one edge of the plate is given by [1]. ln n i n 0 ' ln n e n 0 = e k b T e = q M 2 0 1 1 2 () 2 (3.3) where is measured from the Mach line in the direction toward the wake. The detailed derivation of the analytical quasi-neutral solution and its comparison with hybrid PIC simulations are available in [1, 38]. 3.3 Simulation Setup We consider a thin conducting plate with a size L p = 140 D in a collisionless, unmagnetized, mesothermal plasma flow. The inflow plasma considered are: n i =n e =n 0 , v 0 =C s =M 0 = 8, and ion temperature T i0 = 0:01T e0 , where the subscript “0” denotes the ambient plasma condition. We consider that the plate is unbiased and has a negligible surface potential difference with respect to the ambient plasma potential. In PIC simulations, only half of the plate is included in the simulation domain due to symmetry. The incoming plasma flow is taken to be along the x-direction. The plate is 3.3 Simulation Setup 45 0 50 100 150 200 250 300 0 0.2 0.4 0.6 0.8 1 1.2 0 50 100 150 200 250 300 0 0.2 0.4 0.6 0.8 1 1.2 0 50 100 150 200 250 300 0 0.2 0.4 0.6 0.8 1 1.2 0 50 100 150 200 250 300 0 0.2 0.4 0.6 0.8 1 1.2 Figure 3.2: Instantaneous ion density profiles along the x-axis from the Boltzmann electron fluid based hybrid PIC simulation at ~ y =y= D0 = (a) 0, (b) 35, (c) 70 and (d) 105. placed along the left boundary of the simulation domain (x = 0). The plate potential is set to be the same as the ambient plasma potential ~ p = p =(k b T e0 =e) = 0. In the hybrid PIC simulation, the only particle species is ion. Macro-particles representing the ions are injected from the upstream domain boundary as a drifting Maxwellian while the electrons are simplified by the Boltzmann relation Eq. (2.20). As the ions are cold and supersonic, all the ion particles reaching the outer simulation boundary are removed. The simulation domain used in the hybrid PIC simulation is L x L y = 400 D0 300 D0 . In the full PIC simulation, macro-particles representing both the ions and electrons are injected from the upstream domain boundary as drifting Maxwellian distributions with the same bulk drifting velocity v 0 . Since the electrons are thermal, a re-injection of electrons based on local flow condition at the simulation domain boundary is necessary. However, as the results in Sec. 3.4 will show, the local electron VDF in the wake region deviates from the Maxwellian and cannot be pre-determined. Hence, it is not possible to construct an accurate scheme to re-inject 3.3 Simulation Setup 46 0 50 100 150 200 250 300 0 0.2 0.4 0.6 0.8 1 1.2 0 50 100 150 200 250 300 0 0.2 0.4 0.6 0.8 1 1.2 0 50 100 150 200 250 300 0 0.2 0.4 0.6 0.8 1 1.2 0 50 100 150 200 250 300 0 0.2 0.4 0.6 0.8 1 1.2 Figure 3.3: Instantaneous ion density profiles along the x-axis from the full particle PIC simulation at ~ y =y= D0 = (a) 0, (b) 35, (c) 70 and (d) 105. electrons at simulation boundaries for this problem a priori. Therefore, we choose to use a very long simulation domain along the x-direction, L x L y = 3000 D0 300 D0 , so that the electron front does not reach the simulation boundary during the simulation. Themeshresolutionusedis x = y = D0 . Thetimestep tusedis 0:02! 1 pi0 and 0:1! 1 pe0 in the hybrid and full PIC simulations, respectively. The results presented in the next section will focus on the steady state results in the near wake region, L x L y = 300 D0 150 D0 . The PIC methods is an explicit time-advancing scheme. The steady state are obtained by running the simulations long enough so the results in the region of interested no longer change with time. Figs. 3.2 and 3.3 show the instantaneuous ion density profiles along the x-axis extracted at several locations obtained from the Boltzmann electron fluid based hybrid PIC and the full particle PIC simulations, respectively. Clearly, the steady state is reached when t 60! 1 pi0 and t 2400! 1 pe0 in the hybrid and full PIC simulations, respectively. 3.4 Results and Discussions 47 The PIC simulations performed in this chapter typically use 400 macro-particles per cell for each species in the ambient plasma region, and adopt a uniform weight so the number of macro-particles per cell scales with the local density. To further reduce the statistical noise in both the hybrid and full PIC simulations, over 1200 samples of the flow field are processed after the steady state is achieved to obtain the average flow field properties. The sampling is done once every 5 simulation steps (0:1! 1 pi0 ) and every 40 simulation steps (4! 1 pe0 ) in the hybrid and full PIC simulations, respectively, such that the particles within a cell are refreshed between two successive samplings. 3.4 Results and Discussions 3.4.1 Plasma field: fluid electron vs kinetic electron Fig. 3.4 shows the plasma flow field from the Boltzmann electron fluid based hybrid PIC modeling and the fully kinetic PIC simulation. The left panel shows the ion and electron densities (~ n i = n i =n 0 and ~ n e = n e =n 0 ), showing the hybrid PIC simulation results (color lines), the difference between the hybrid and full PIC (color contours), and the analytical solution of Eq. (3.3) (grey dashed lines). The right panel shows the potential ~ = =(k b T e0 =e) (color contour), electric field vector (yellow arrows) and ion streamlines (black solid lines) obtained from the hybrid and full PIC simulations. As predicted by the analytical solution, an expansion fan structure is generated by the edge of the plate. The plasma expansion starts at the Mach line (marked as the iso-density lines of ~ n i;e = 1). The Mach line angle from the PIC simulation results matches exactly the analytic solution of 0 ' sin 1 (1=M 0 )' 1=M 0 . Both the plasma density and electric potential decrease from the Mach line towards the wake. A potential well is established behind the plate in the wake region. From Fig. 3.4, we find that the analytical solution agrees very well with both the hybrid and full PIC simulations within the region where the plasma density ~ n i;e 0:3. In the region ~ n i;e < 0:3, while the analytical solution still mostly agrees with the PIC simulation results, the plasma iso-density lines from the PIC simulations start to exhibit a curvature. The curvature 3.4 Results and Discussions 48 Figure 3.4: Plasma wake. Left panel: (a) ion density and (b) electron density. The color lines are the iso-density lines from the Boltzmann electron fluid based hybrid PIC simulation. The color contours are the density errors of the hybrid PIC results against the full PIC results defined by Eq. (3.4). The gray dashed lines the analytical solution of Eq. (3.3). Right panel: potential contours (color), electric field vectors (yellow arrows) and ion streamlines (black solid lines) from the (c) Boltzmann electron fluid based hybrid PIC and (d) full particle PIC. in the far wake region is due to the counter streaming of ions, which are clearly shown by the ion streamlines in the PIC results but unaccounted for in the fluid based analytic solution. To examine the difference between the hybrid and full PIC simulations, we define the hybrid PIC error as the local density difference between the hybrid PIC and the full PIC ~ n i;err = j~ n i;hybrid ~ n i;full j ~ n i;hybrid ; ~ n e;err = j~ n e;hybrid ~ n e;full j ~ n e;hybrid (3.4) where the subscripts hybrid and full denote the results from the Boltzmann electron fluid based hybrid and full PIC simulations, respectively. The ~ n i;err contours for ~ n i;err 10% are shown in Fig. 3.4 (a), and the ~ n e;err contours for ~ n e;err 10% are shown in Fig. 3.4 (b). We find that the discrepancy between the hybrid and full PIC results is very small until the iso-density ines start to show significant curvature. It is interesting to note that the ~ n i;err = ~ n e;err = 10% line approximately coincides with the ~ n i ' ~ n e ' 0:15 iso-density line, which also matches approximately the expansion characteristics line at = 2 0 . Consistent with the plasma density results, the potential from the hybrid PIC simulation also agrees 3.4 Results and Discussions 49 1.0 0.6 0.2 0.1 0.05 0 50 100 150 200 250 300 0 50 100 150 0.9 0.95 1 1.05 1.1 1.15 1.2 0.022 0.0225 0.023 0.0235 0.024 0.0245 0.025 Figure 3.5: (a) Ion density ~ n i contours from the hybrid PIC simulations with various values and the full particle PIC simulation. (b) Weighted root mean square error of the hybrid PIC results against the full PIC result defined by Eq. (3.5) as a function of . with that from the full PIC simulation within the same region, The initial expansion fan region within the ' 2 0 characteristics line is also marked in Fig. 3.4. The effects of electron cooling on the plasma wake have been examined by further carrying out the hybrid PIC simulations based on the polytropic relation, Eq. (2.22), for the electron dynamics. Fig. 3.5 (a) shows the comparison of the ion density contour in the near wake region between the hybrid PIC simulations with various values and the full PIC simulation. Note that the case of = 1 correponds to the Boltzmann electron fluid model. Clearly, even through the cooling/heating effects are taken into account, the electron fluid model based on the polytropic law with a constant is not able to fully correct the errors. When ~ n i 0:2, the difference between the hyrbid PIC results with various values and the full PIC results is not obvious. A close examination shows that the hybrid PIC with 1 agrees better with the full PIC comparing to other values. When ~ n i 0:1, the ~ n i isolines from all the hybrid PIC simulations start to exhibit notable discrepancies in the deep wake region. The increase of value slightly improve the agreement between the hybrid and full PIC results when ~ n i 0:05. This implies an effective overall electron cooling in that region. 3.4 Results and Discussions 50 To quantitatively measure the difference between the hybrid and full PIC results, a weighted root mean square relative error at each mesh point (ix,iy) over the entire near wake region, L x L y = 300 D0 150 D0 , is employed err(~ n i ) = v u u t X ix;iy w ix;iy ~ n i;hybrid (ix; iy) ~ n i;full (ix; iy) ~ n i;hybrid (ix; iy) 2 (3.5) where w ix;iy = ~ n i;hybrid (ix; iy)= P ~ n i;hybrid (ix; iy). Eq. (3.5) was used for analyzing the errors of a two-fluid modeling and full PIC simulation in the plume expansion [77]. Fig. 3.5 (b) shows the accuracy of the hybrid PIC measured by Eq. (3.5) as a function of . It is found that, though none of the hybrid PIC simulations reproduces the full PIC result, = 1:02 yields the least error. 3.4.2 Kinetic characteristics and fluid approximation breakdown of electron The difference between the hybrid PIC and full PIC results is due to the electron model. The full PIC simulation resolves electron kinetic characteristics. Fig. 3.6 shows both the electron temperature (color contours) and electron “streamlines” (blacklines) obtained from the full PIC simulation. Here, the local electron temperature T e at each mesh point (ix,iy) refers to the second electron velocity moment averaged over the macro-particles within each cell, and is calculated from ~ T e;j (ix; iy) = ~ v 2 te;j (ix; iy) =h~ v 2 j ih~ v j i 2 ; h~ v n j i = Np X p=1 ~ v n j;p =N p (3.6) where j =x or y, n = 1 or 2, ~ v j;p is individual particle velocity, and N p the total number of macro-particles in a given cell. ~ v j =v j = p k b T e0 =m e and ~ T e;j =T e;j =T e0 are the normalized velocity and electron temperature, respectively. The electron temperature contours in Fig. 3.6 show that ~ T e;x and ~ T e;y behave very differently in the wake region. In the initial expansion fan region, ~ T e;x is approximately 3.4 Results and Discussions 51 Figure 3.6: Electron temperature (color contours) and streamlines (black solid lines) from full PIC simulation: (a) ~ T e;x and (b) ~ T e;y . The magenta solid line shows the boundary between electron forward-flow and back-flow. constant with T e;x ' T e0 while ~ T e;y slowly decreases from T e0 as the plasma expands into the wake region. As the plasma is collisionless, the anistropic electron temperature behavior reflects the different effects of the electric field on electron motion in different directions. The electric field in the wake region is mostly in the y-direction. In the y-direction, since the electrons need to overcome the electrostatic force as they expand into the wake, this leads to a decrease in electron kinetic energy and thus ~ T e;y . In the x-direction, since there is neither electrostatic force nor collisions to affect electron thermal velocities, ~ T e;x remains approximately constant. In the deep wake region, both ~ T e;x and ~ T e;y exhibit a complex pattern. ~ T e;x shows a significant decrease immediately behind the plate and a moderate increase further downstream along the x-direction. ~ T e;y shows a moderate decrease near the edge behind the plate and a significant increase near the center along the y-direction. To understand this pattern, Fig. 3.7 further shows the local electron VDFs ( ^ f S ) by the color surface plots along the x = 0 and y = 0 axis. Additionally, with the bulk velocity and temperature computed from the macro-particles, we also construct the Maxwellian VDFs ( ^ f M ), shown by the black lines at selected locations, assuming a local thermal equilibrium state for electrons at these locations. In the x-direction, as the electrons enter into the wake, those with insufficient energy to overcome the potential well in the wake are attracted to the plate while those with higher energy flow downstream (see electron stream lines). The division between the back-flow and forward-flow is marked approximately by the magenta color line in Fig. 3.6. As the electrons 3.4 Results and Discussions 52 Figure 3.7: Local electron VDFs ( ^ f S ) from the full PIC simulation (color surface plots) and fitted Maxwellian VDFs ( ^ f M ) based on local ~ v d and ~ v t (black solid lines). (a) ^ f(~ v x ; ~ x) along ~ y = 0, (b) ^ f(~ v x ; ~ y) along ~ x = 0, (c) ^ f(~ v y ; ~ x) along ~ y = 0, and (d) ^ f(~ v y ; ~ y) along ~ x = 0. The insets highlight the comparison between the actual VDF ( ^ f S ) and the fitted Maxwellian ( ^ f M ) at several representative positions. are collisionless, the ~ T e;x pattern observed simply reflects the fact that the region immediately behind the plate is populated by the electrons with lower initial kinetic energies in x and the region further downstream by the electrons with higher initial energies in x. In the y-direction, the decrease in ~ T e;y near the edge behind the plate is a continuation of electron cooling through the expansion fan. Near the center behind the wake (~ y 0), the significant increase of ~ T e;y is not a result of actual electron heating. Rather, it is due to the counter streaming of electrons. As Figs. 3.7 (c) and (d) clearly show, ^ f S (~ v y ) is a single population near the edge of plate but a double-peak shaped distribution from two counter-streaming electron populations near ~ y 0. A Maxwellian VDF fit to the double-peaked distribution would lead to an artificially high electron temperature. This offers a possible explanation 3.4 Results and Discussions 53 0 50 100 150 200 250 300 0 50 100 150 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 0.02 0.04 0.06 0.08 0 0.2 0.4 0.6 0.8 1 Figure 3.8: (a) Contour of non-equilibrium parameter overlaid with the iso-lines of ~ n i;err = 10% (pink dotted) and ~ n e;err = 10% (black solid). (b) Dependence of the conditional averages of hybrid PIC prediction errors on . for the observed anomalous electron temperature enhancement in the wake of ionospheric satellites [133–136]. Figs. 3.7 shows that the electron VDF in the wake region is anisotropic and strongly non-Maxwellian. Hence, a fluid approximation for electrons in general is not valid in that region. In order to evaluate the effects of non-equilibrium electrons on macroscopic plasma flow field, we further define a non-equilibrium parameter as = 1 ~ T e X j=x;y ~ T e;j Z 1 1 j ^ f M;j ^ f S;j j ^ f M;j d~ v j (3.7) where ~ T e;j is defined in Eq. (3.6), ~ T e = ( ~ T e;x + ~ T e;y )=2, and ^ f M;j = ^ f M (~ v j ) and ^ f S;j = ^ f S (~ v j ) are the assumed Maxwellian and actual VDF correspoding to ~ v j ,respectively. is the sum of the temperature weighted mean deviation of the Maxwellian VDF ^ f M from the actual VDF ^ f S in the x- and y-directions. We note that some other non-equilibrium or fluid breakdown parameters were used in the studies of rarefied gas flows [e.g. 137–140]. The shortcoming of those parameters is that only the macroscopic flow field information was considered. In contrast, the non-equilibrium parameter defined in Eq. (3.7) takes into account the temperature anisotropy and the local VDF’s deviation from the Maxwellian. Fig. 3.8 (a) shows the contour of for > c = 0:02, where c is taken to be a threshold that the effects of non-equilibrium electrons on macroscopic flow field becomes significant. The regions where the hybrid PIC starts to deviate from the full PIC, ~ n i;err = 10% (pink 3.4 Results and Discussions 54 dotted) and ~ n e;err = 10% (black solid), are also shown in Fig. 3.8 (a). Fig. 3.8 shows that the region where the Boltzmann electron fluid based hybrid PIC becomes different from the full particle PIC simulation falls into the non-equilibrium region ( >= 0:02). More importantly, the shape of the boundary as defined by ~ n i;err and ~ n e;err = 10% coincides very well with that defined by = 0:02. Fig. 3.8 (b) further shows the calculated average of the hybrid PIC errors (both ~ n i;err and ~ n e;err ) based on the binned value of over the entire plasma region (both ~ n i and ~ n e > 0). The conditional averages of the hybrid PIC errors,h~ n i;err j i andh~ n e;err j i, are also found to scale well with . Based on the strong correlation between the hybrid PIC error and the non-equilibrium parameter , we conclude that, in a collisionless, unmagnetized, mesothermal plasma flow, the plasma wake may be characterized into two regions based on electron behavior. Within the initial expansion fan region from the Mach line to ' 2 0 ' 2=M 0 , the effects from electron kinetic characteristics on macroscopic flow are not significant. While the electrons are not exactly isothermal, an isothermal fluid electron based approach, either analytical or hybrid PIC, leads to a plasma flow field similar to that from the full PIC simulation. Hence, this region may be considered as the “fluid” electron expansion region where a fluid approximation for electrons can be applied. Beyond ' 2 0 ' 2=M 0 , the local electron VDFs are strongly non-Maxwellian, and the electron “temperature” is strongly anisotropic and exhibits both “cooling” and “heating” in the wake. Since the energy transfer between the electrons is not through inter-particle collisions, a more sophisticated fluid based approach than that considered here (such as replacing the isothermal condition with the polytropic law, T e n 1 e = const, and incorporating the electron energy equation) still cannot be applied to accurately describe the macroscopic electron properties. Therefore, this region is the kinetic electron expansion region where a fluid approximation for electrons breaks down and a fully kinetic approach must be adopted. 3.5 Summary 55 3.5 Summary We studied the plasma wake formation behind a large plate in a collisionless, unmagnetized, mesothermal plasma through the fluid based analytic solution, hybrid particle-ion fluid- electron PIC, and fully kinetic PIC simulations. The focus is to investigate the validity of the commonly used assumption that the electrons may be considered as a fluid in such a problem. The electron kinetic effects on macroscopic plasma flow are quantified using both the hybrid PIC error and a non-equilibrium parameter based on the deviation of the electron VDF from the Maxwellian VDF. We find that the plasma wake may be characterized into two regions based on electron behavior, a fluid expansion region between the initial Mach line and the expansion fan angle at ' 2 0 ' 2=M 0 , and a kinetic electron expansion region beyond' 2 0 ' 2=M 0 . In the fluid electron expansion region, the effects of electron kinetic characteristics are not significant, all three approaches lead to a similar result, and thus a fluid approximation for electrons may be applied. The perturbation generated by the plate in this region is analogous to that of the Prandtl-Meyer expansion of a supersonic gas flow. In the kinetic electron expansion region, the fluid approximation for electrons breaks down and a fully kinetic approach must be adopted. Full particle PIC simulations show that the local electron VDFs in this region are strongly non-Maxwellian and the electron “temperature” is strongly anisotropic and exhibits a complex pattern of “cooling” and “heating”. The effects from a biased object and the magnetic field on the conclusions of this paper will be addressed in future studies. 56 Chapter 4: Collisionless, Mesothermal Plasma Flow over a Large Biased Conducting Plate: Effects of Surface Charging on the Electron Fluid Approximation The last chapter investigates the plasma physics in the wake of a large plate of zero potential with respect to the ambient plasma. In this chapter, the surface charging of the plate and its interactions with the collisionless, mesothermal plasma flow are studied with by using both the hybrid and full PIC simulations. The formation and the steady-state structure of the space-charge wake as well as the ion impact and current collection are presented. In particular, the validity of the fluid approximation for electrons is assessed by the direct comparsion between the hybrid and full PIC simulations. 4.1 Introduction Many space platforms, such as a solar array or a biased radiator, attached to space systems in orbit are charged and can interact with the environmental plasmas. The structures of the plasma field and thus the interaction between the plasma flow and the space platforms or the spacecraft mainly depends on the dimension of platform L p , the surface potential p and the sheath thickness d sh . According to the three pamaters, the collisionless plasma flow over a space plateform may be divided into the “quasi-neutral” regime, the “thin-sheath” regime, and the “thick-sheath” regime [1], as shown in Fig. 4.1. The flow characteristics of the plasma is dinstinct in each of the regime. The quasi-neutral limit corresponds to the case 4.1 Introduction 57 Figure 4.1: Schematic of the plasma field around a large biased plateform in the (a) “quasi- neutral”, (b) “thin sheath”, and (c)“thick sheath” regimes [1]. we deal with in Chapter 3. In the initial expansion region, the wake is simply composed of a expansion fan similar to the Prandtl-Meyer expansion wave and can be accurately described by the analytical solution, Eq. (3.3), as well as the fluid model. In the thin sheath regime of a negatively biased plate, the front face forms a thin sheath, d sh L p , which can be approximated by the Child-Langmuir law [12]. On the wake side, two embedded sheaths generated at the edges of the plate will modify the expansion struture. As the plate becomes more negatively charged, the embedded sheath can be bendded to form the “hook”-shaped structure and fill in the ion void region. It is speculated that the “hook”-shaped sheath will eventually disappear when the collisionless, mesothermal plasma flow over a large high-voltage platform falls into the thick sheath regime where d sh L p . In this regime, only a single sheath exists in the near wake region, as illustrated in Fig. 4.1 (c). In the “thin-sheath” and “thick-sheath” regimes, analytical solutions are difficult to be obtained and numerical simulations are the primary tool for investigations. The dimension of the space platforms and spacecraft in low Earth orbit (LEO) is mostly in the range of 1 m L p 20 m. The surface potential is typically in the range of 1 V. p . 10 2 V. In a typical LEO environment, the Debye length D0 is 0:01 m as the 4.2 Simulation Model 58 plasma densityn 0 is 210 11 m 3 and the electron temperatureT e0 is 0:1 eV. Accordingly, the sheath ratio =d sh =L p varies from 1 to 1 if the sheath thicknessd sh is estimated according to the Child-Langmuir law [12], d sh D0 ' 0:8 p M 0 j ~ p j 3=4 (4.1) where ~ p = p =(k b T e0 =e) is the normalized surface potential, and the ambient Mach number M 0 is typically 8 in LEO. Based on these information, it can be seen that the interaction of the space plasmas with the spacecraft in LEO can involve all the three regimes. This chapter carries out both the hybrid and fully kinetic PIC simulations on the collision- less, mesothermal plasma flow over a negatively biased large thin plate. The surface potential varies in a wide range so all the three regimes, namely the quasi-neutral, thin sheat, and thick regimes, can be included. The main purpose is to investigate the effects of surface potential or the sheath ratio on the interactions between the space plasma and the space platform. The validity of applying the commonly used electron fluid approximation and the hybrid model to the collisionless, mesothermal plasma flow over a negatively biased conducting object in a wide range of surface potential will be assessed by the high fidelity fully kinetic method. 4.2 Simulation Model This chapter considers a large thin conducting plate with the length L p = 150 D0 and the thickness d p = 2 D0 as a model space platform. The surface potential, ~ p = p =(k b T e0 =e), is considered to change from -1 to -2000 so that all the plasma space-platform interaction regimes, from the “quasi-neutral” regime, “thin sheath” regime to “thick sheath” regime, are included. The ambient plasma flow condition is the same as that in Chapter 3, namely,n i =n e =n 0 , v 0 =C s = M 0 = 8, and T i0 = 0:01T e0 , where the subscript “0” denotes the ambient plasma condition. 4.2 Simulation Model 59 Figure 4.2: Schematics of the simulation setup for a collisionless, mesothermal flow over a negatively biased large thin conducting plate. The geometric setup is slightly different from that in the last chapter. The simulation domain L x L y is set to be 800 D0 400 D0 and 4000 D0 400 D0 for the hybrid PIC and fully kinetic PIC simulations, respectively. As shown in Fig. 4.2, the ambient plasma is injected from the inflow boundary at ~ x =100 and the plate is put at ~ x = 100, 200 D0 downstream of the inlet boundary, to avoid the sheath of plate interacting with the inlet boundary for all the values of ~ p considered. Only half of the plate is included in the simulation domain due to symmetry, same as the last chapter. The boundary conditions are also set to be the same as those in the last chapter, illustrated in Fig. 4.2. The mesh resolution used is x = y = D0 . The time step t used is 0:02! 1 pi0 and 0:1! 1 pe0 in the hybrid and full PIC simulations, respectively. We are interested in the steady state results in the sheath region upstream of the plate and the near wake region. Therefore, the results presented in the next section will focus on the region of (~ x; ~ y) = (0 600; 0 200). Based on our monitoring, a steady state in the region of interest is guaranteed for all the cases whent 120! 1 pi0 andt 6000! 1 pe0 for the hybrid and full PIC simulations, respectively. Both the hybrid and full PIC simulations use 400 macro-particles per cell for each species in the ambient plasma region. A uniform weight is adopted so the number of macro-particles per cell scales with the local density. Over 800 flow fields are sampled and averaged after the steady state is achieved to further reduce the statistical noise. 4.3 Results and Discussion 60 4.3 Results and Discussion 4.3.1 Plasma field Fig. 4.3 shows the hybrid PIC results of ion (~ n i ) and electron (~ n e ) density by the color isolines and its relative errors against the full PIC results for different surface potentials. Fig. 4.4 shows the electric potential, electric field vectors and ion streamlines. It is clear that the plasma flow field around the plate is greatly influenced by the change of the plate surface potential. When the plate is slightly negatively biased ( ~ p =2), the plasma wake structure is very similar to that associated with an unbiased plate presented in the last chapter despite that a very thin sheath forms on the front side of the plate. As the plate gets more negatively biased, the sheath becomes thicker at the front side. On the wake side, the remarkable phenomenon on the ion profiles is the formation of the “embedded” sheath generated at the plate edge first recognized by Wang and Hastings [1]. The “embedded” sheath, as clearly seen in Fig. 4.3, is an ion density enhanced stream penetrating into the ion void region and separating it from the expansion fan region. As the plate surface potential further decreases, the “embedded” sheath starts to be bent, forming a “hook” shaped structure. Meanwhile, an electron void region forms around the plate because ~ p ~ T e . When ~ p 80, the core (maximum ion density) of the “embedded” sheath is attached to the plate. The bending of the “embedded” sheath and its attachment to the plate are also reflected by the ion streamlines in Fig. 4.4. The core of the “embedded” sheath on the back face of the plate moves towards the center of the plate as the potential keeps decreasing. When the potential is negative enough, the two ion streams generated at the upper and lower edges of the plate will cross each others and create “secondary” streams hitting the plate surface again, as seen from the case of ~ p =2000. The crossing of ion streams from the upper and lower sides of the plate is clearly shown by the ion streamlines. In addition, the electron void region around the plate gets larger and the sheath gets thicker as the surface potential becomes more negative. In the thick sheath regime, an expansion fan region still exists. However, it is created at the edge of 4.3 Results and Discussion 61 thick sheath, instead of the plate edge. This is different from the speculation illustrated in Fig. 4.1, which is probably true for the situation of thermal, low speed plasma flows. The color contours in Fig. 4.3 show the hybrid prediction errors of the ion and electron densities against the full PIC results calculated according to Eq. (3.4). Same as the last chapter, only the regions of ~ n i;err %10 and ~ n e;err %10 are shown. For all the plate surface potential considered here, the hybrid PIC results agree well with the full PIC result on the ram side of the plate. In the wake, the hybrid PIC predictions of both ~ n i and ~ n e agree well with the full PIC results in the expansion fan region, which is identified as the fluid expansion region in the last chapter. It is surprising, at the first appearance, to see that the hybrid PIC simulation is able to accurately predict the ion profiles, especially the “embedded sheath”, near the plate. The success of the hybrid PIC near the plate is mainly due to the absence of electrons in the wake next to the largely negatively biased plate. The failure of hybrid PIC simulation mainly happens in the deep expansion region. It is interesting to note that the accuracy of the hybrid PIC simulation in the deep expansion region changes non-monotonically with ~ p . The smallest discrepancy occurs at ~ p '10. The potential contour shown in Fig. 4.4 may be able to provide a clue for the reason of the non-monotonic dependence of the hybrid prediction errors on ~ p . For ~ p 5, there is a potential well in the wake of the plate. As ~ p becomes more negative, the depth of the potential decreases. For ~ p 10, no potential well exists in the wake of the plate. As discussed in the last chapter, the potential well plays an important role in the unusual observations on the electron macroscopic flow quantities. Therefore, the discrepancy of the hybrid PIC simulation from the full PIC simulation in the wake is initially alleviated as ~ p becomes more negative, and achieves the minimum around ~ p =10. For ~ p > 10, the electron void region is becoming larger as ~ p keeps decreasing. The border of where the electron kinetic effects are important becomes larger. This is verified by the severe hybrid errors for the electron density near the ~ n e = 0 isoline. Therefore, for ~ p > 10, the area where the hybrid PIC simulation is inapplicable is getting larger as ~ p becomes more negative. 4.3 Results and Discussion 62 (a) ~ p =2 (b) ~ p =5 (c) ~ p =10 (d) ~ p =20 (e) ~ p =40 (f) ~ p =80 (g) ~ p =200 (h) ~ p =1000 (i) ~ p =2000 Figure 4.3: Ion and electron densities (color isolines) from the hybrid PIC simulations, and density errors of the hybrid PIC results against the full PIC results (color contours) defined by Eq. (3.4). 4.3 Results and Discussion 63 (a) ~ p =2 (b) ~ p =5 (c) ~ p =10 (d) ~ p =20 (e) ~ p =40 (f) ~ p =80 (g) ~ p =200 (h) ~ p =1000 (i) ~ p =2000 Figure 4.4: Contours of normalized potential ~ = ~ p (color), electric field vectors (yellow arrows) and ion streamlines (black solid lines) from both hybrid and full PIC simulations. The region enclosed by the white line is where ~ < ~ p , indicating that there is a potential well in the wake. The potential well structure in the wake appears only for the cases of ~ p =2 and5. 4.3 Results and Discussion 64 4.3.2 Ion impact and current collection The ion impact and current collection on the plate surface is of crucial importance for the spacecraft design and operation. Fig. 4.5 shows the ion current density profiles on both the front and back faces of the plate from both the hybrid and full PIC simulations. In general, the hybrid PIC results agree well with the full PIC results. When the plate surface potential is slightly lower than the plasma potential ( ~ p 5), the ion collection only happens on the front face and the current density equals en i v 0 , the prediction from the thin sheath theory. When the surface potential has a moderate negative value (10 ~ p 20), the back face of the plate is still ion void so the ions are only collected on the front face. However, near the corner on the front face, the plate collects slightly more ions than the results of the thin sheath theory due to the edge effect. When ~ p decreases to -40, the potential difference is big enough to attract some ions back to the plate. The back face is observed to collect ions, with a current density plateau near the center of the plate. For80 ~ p 200, a current density peak is seen somewhere between the edge and the center of the back face due to the “embedded” sheath attachment. The position of the current density peak moves towards the center of the plate as ~ p decreases. On the front face, the edge effect, which increases the ion current density compared to the theoretical value, becomes more significant as the surface potential gets more negative. The flow enters the “thick sheath” regime when ~ p is reduced to1000. On the front side, the ion current density is higher than the value predicted by the “thin sheath” theory for the entire plate. From the edge to the center, the ion current density slightly rises first, getting the maximum a few Debye length away, and then slowly decreases to the minimum. On the back side, an extremely high current density is observed near the center of the plate, to which the very high density “embedded” sheath core is attached. As ~ p further reduces to2000, the pattern of ion current density distribution on the front side is similar to that for the ~ p =1000 case, but with a slightly higher value. On the back side, it is observed that the location of the peaked ion current density moves towards the edge as a result of the “secondary” ion 4.3 Results and Discussion 65 0 0.2 0.4 0.6 0.8 1 1.2 0 15 30 45 60 75 front - hybrid PIC front - full PIC back - hybrid PIC back - full PIC (a) ~ p =2 0 0.2 0.4 0.6 0.8 1 1.2 0 15 30 45 60 75 front - hybrid PIC front - full PIC back - hybrid PIC back - full PIC (b) ~ p =5 0 0.2 0.4 0.6 0.8 1 1.2 0 15 30 45 60 75 front - hybrid PIC front - full PIC back - hybrid PIC back - full PIC (c) ~ p =10 0 0.2 0.4 0.6 0.8 1 1.2 0 15 30 45 60 75 front - hybrid PIC front - full PIC back - hybrid PIC back - full PIC (d) ~ p =20 0 0.2 0.4 0.6 0.8 1 1.2 0 15 30 45 60 75 front - hybrid PIC front - full PIC back - hybrid PIC back - full PIC (e) ~ p =40 0 0.2 0.4 0.6 0.8 1 1.2 0 15 30 45 60 75 front - hybrid PIC front - full PIC back - hybrid PIC back - full PIC (f) ~ p =80 0 0.2 0.4 0.6 0.8 1 1.2 0 15 30 45 60 75 front - hybrid PIC front - full PIC back - hybrid PIC back - full PIC (g) ~ p =200 0 1 2 3 4 5 6 0 15 30 45 60 75 front - hybrid PIC front - full PIC back - hybrid PIC back - full PIC (h) ~ p =1000 0 1 2 3 4 5 6 0 15 30 45 60 75 front - hybrid PIC front - full PIC back - hybrid PIC back - full PIC (i) ~ p =2000 Figure 4.5: Ion current density profiles on the front and back faces of the plate. stream hitting the plate mentioned in the iso-density plot. Note that the scales of J i =en 0 v 0 for ~ p =1000 and2000 in Fig. 4.5 are different from the others in the same figure. The ion current collected by the plate can be obtained through the integration of current density profile shown in Fig. 4.5. Fig. 4.6 (a) shows the dependence of ion current collection per unit length on the surface potential for both the hybrid and full PIC simulations. It is worth mentioning that the “thin sheath” theory gives ~ I= ~ L p = 1 at the front side. The front side ion current collection initially increases very slowly as the potential decreases. The contribution to this increase mainly comes from the edge effects. When entering the “thick sheath” regime, the front side current collection increases more rapidly. For ~ p =2000, 4.3 Results and Discussion 66 -2000-1000 -200 -80 -40 -20 -10 -5 -2 0 0.5 1 1.5 2 2.5 -2000-1000 -200 -80 -40 -20 -10 -5 -2 0 0.2 0.4 0.6 0.8 1 Figure 4.6: Effects of surface potential on: (a) ion current collection, (b) Wake side current ratio. about 50% increase of the ion current collection is achieved, compared with the “thin sheath” theory prediction. The back side collects ion when ~ p 40. The amount of back side ion collection is shown to increase faster than its front side counterpart as the surface potential decreases. This can be seen from the ratio of the back side to front side current collection C w = I back =I front shown in Fig. 4.6 (b). It is clear that C w increases monotonically as ~ p decreases. In particular, C w scales logarithmically with ~ p in the “thick sheath” regime. The scaling factor is found to be 0:16 The present simulation results show that the back side collects as much as about 60% of the current collected by front side for ~ p =2000. 4.3.3 Kinetic characteristics and fluid approximation breakdown of electron It is shown that the hybrid PIC simulations give very good predictions for the ion dynamics in the expansion fan region and the near plate region, including the density profiles and the current collection by the plate. However, in the deep expansion region, the hybrid PIC results are inadequate, especially for the cases with low and high negative surface potentials. The difference between the hybrid PIC and full PIC results is due to the electron model. The full PIC simulation resolves electron kinetic characteristics. Fig. 4.7 shows both the electron temperatures (color contours) and electron streamlines (blacklines) obtained from the full PIC simulations for different surface potentials. The 4.3 Results and Discussion 67 pattern of electron temperature contours in the wake of the plate for the low voltage cases ( ~ p =2 and5) is quite similar to that for the unbiased plate in the last chapter. The cause of the such a pattern is the potential well which serves to filter the electrons with low and high energies and kicks them into different regions. A detailed explanation is presented in Sec. 3.4.2. ~ T e;x in the deep expansion region is initially high when ~ p =2, then decreases and keeps a relatively low value when ~ p =1040. After ~ p 80, it increases again and gets a very high value when ~ p =2000. The change in ~ T e;y in the deep expansion region shows the opposite trend. The interesting phenomenon observed here is that the hybrid PIC prediction errors, in terms of ~ n i;err and ~ n e;err , are significant when ~ T e;x is high or ~ T e;y is low in the deep expansion region. For instance, Fig. 4.3 shows that the area with ~ n i;err or ~ n e;err 10% is large for ~ p =2 and2000. In the mean time, it is observed from Fig. 4.7 that ~ T e;x increase is significant in the deep expansion region for the same values of ~ p . In contrast, for ~ p =1020, the hybrid PIC simulation is quite satisfactory, and ~ T e;x enhancement is observed to be mild in the deep expansion region. Sec. 3.4.2 claims that the failure of the hybrid PIC modeling is due to the breakdown of the electron fluid approximation, and defines a microscopic based parameter to quantify the non-equilibrium degree, Eq. (3.7). Fig. 4.8 examines the reason for the hybrid PIC failure with various surface potential by showing the contour of the non-equilibrium parameter defined in Eq. (3.7) and the border (~ n i;err and ~ n e;err = 10%) of the hybrid PIC breakdown region. The comparison shows that the regions where the hybrid PIC starts to deviate from the full PIC, ~ n i;err = 10% (pink dotted) and ~ n e;err = 10% (black solid), falls into the non-equilibrium region ( crit ), where the threshold crit is set to be 0.024 here. We note that the shape of the hybrid PIC breakdown border (10% isolines of ~ n i;err and ~ n e;err ) in general follows the crit boundary very well in the wake region, despite that slightly overestimate the area of hybrid PIC failure region for the cases of moderate voltage ( ~ p '=40). These facts confirms that the failure of the hybrid PIC modeling is primarily due to the breakdown of electron fluid approximation and that is in general an appropriate parameter to quantify the non-equilibrium degree of the flow. 4.3 Results and Discussion 68 (a) ~ p =2 (b) ~ p =5 (c) ~ p =10 (d) ~ p =20 (e) ~ p =40 (f) ~ p =80 (g) ~ p =200 (h) ~ p =1000 (i) ~ p =2000 Figure 4.7: Electron temperatures (color contours) and streamlines (black solid lines) from full PIC simulations. 4.3 Results and Discussion 69 0 100 200 300 400 500 600 0 50 100 150 200 0.024 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (a) ~ p =2 0 100 200 300 400 500 600 0 50 100 150 200 0.024 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (b) ~ p =5 0 100 200 300 400 500 600 0 50 100 150 200 0.024 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (c) ~ p =10 0 100 200 300 400 500 600 0 50 100 150 200 0.024 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (d) ~ p =20 0 100 200 300 400 500 600 0 50 100 150 200 0.024 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (e) ~ p =40 0 100 200 300 400 500 600 0 50 100 150 200 0.024 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (f) ~ p =80 0 100 200 300 400 500 600 0 50 100 150 200 0.024 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (g) ~ p =200 0 100 200 300 400 500 600 0 50 100 150 200 0.024 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (h) ~ p =1000 0 100 200 300 400 500 600 0 50 100 150 200 0.024 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 (i) ~ p =2000 Figure 4.8: Contours of non-equilibrium parameter overlaid with the isolines of ~ n i;err = 10% (pink dotted) and ~ n e;err = 10% (black solid). 4.4 Summary 70 4.4 Summary The collisionless, mesothermal plasma flow over a large negatively biased plate with a wide range of surface potential is studied by using both the hybrid and fully kinetic PIC simulations. The plasma wake structures, ion impact and current collection, and the electron kinetics in both the “thin sheath” and “thick sheath” regimes are investigated. In particular, the present simulations of collisionless mesothermal plasma flows have shown that the quasi-neutral expansion fan still exists even in the “thick sheath” regime, in contrast to the speculation in the literature. The main focus is to investigate the effects of surface potential on the validity of the isothermal electron fluid approximation. It is found that the hybrid PIC modeling is able to describe the ion dynamics in the vicinity of the plate very well. The difference bewteen the hybrid and full PIC results of the ion current collection is very small for all the surface potential values considered in this study. The electron kinetic effects on the macroscopic plasma flow in the expansion region are quantified using both the hybrid PIC error and the non-equilibrium parameter based on the deviation of the electron VDF from the Maxwellian VDF. We find that the area where the hybrid PIC errors are significant (> 10%) first decreases and then increases as the surface potential ~ p changes from2 to2000. Together with this trend, the enhancement of ~ T e;x or the dropping of ~ T e;y , compared with the ambient ~ T e , in the deep expansion region is observed to first decrease and then increase. The failure of the hybrid PIC modeling is primarily due to the breakdown of electron fluid approximation. This is confirmed by the local non-equilibrium degree of electron flow measured by the non-equilibrium parameter , which is shown to have a very good generalizability. 71 Chapter 5: Expansion of a Collisionless Hypersonic Plume into Vacuum: Fully Kinetic and Hybrid Particle-in-Cell Simulations In this chapter, we conduct both the fully kinetic and hybrid particle-in-cell (PIC) simulations to investigate the two-dimensional (2D) expansion of a collisionless, collisionless plasma plume into a vacuum. The fully kinetic PIC simulations are carried out using the real ion-to-electron mass ratios of H + , Ar + and Xe + while the hybrid PIC model assumes the electrons to be a massless, isothermal fluid. The electron fluid model used in the hybrid PIC is the Boltzmann electron fluid model [Eq.(2.20)]. The simulations show that the hypersonic plasma plume exhibits four distinct regions, that is, the unperturbed, quasi-steady expansion, self-similar expansion and electron front regions. The kinetic behavior of electrons is strongly anisotropic, causing considerably different expansion characteristics in the plume direction and the transverse direction. Along the plume direction, the expansion dynamics is demonstrated to be similar to that of the classical one-dimensional (1D) semi-infinite plasma expansion in which the electrons may be considered almost isothermal. In the transverse direction, the expansion process can be analogous to the 1D expansion of a finite plasma in which the effect of electron cooling is important. This anisotropic characteristic is attributed to the amount of electron thermal energy available from the source in different directions. The direct comparison between the hybrid and full PIC simulations shows that the widely used equilibrium isothermal electron fluid model is not valid in studying the expansion of a collisionless plasma plume. 5.1 Introduction 72 5.1 Introduction This chapter considers the expansion of a collisionless mesothermal plasma plume into a vacuum. The plasma plume at the emission source is quasi-neutral and has a plume velocity v 0 for both the ions and electrons with a flow characteristic of v ti0 v 0 v te0 , where v ti0 and v te0 denote the ion and electron thermal velocity of the emission source, respectively. In addition, the plume Mach number M 0 =v 0 =C s0 is usually much larger than 1, and thus can be considered to be hypersonic. The associated expansion process not only is a fundamental problem in plasma dynamics but also appears in many plasma engineering applications. An important example is the plasma plume emitted by an electric propulsion (EP) thruster. An EP thruster typically emits a low density plasma plume of cold beam ions and thermal electrons with a mesothermal flow characteristic [10, 59]. The plume becomes quasi-neutral almost immediately downstream of the thruster exit. The collision mean free path is also much larger than the flow characteristic length of the plume. The expansion process not only determines the plume structure but also plays an important role in EP plume-spacecraft interactions [59, 83, 141]. As a classical topic in plasma physics, plasma expansion has been studied extensively using one-dimensional (1D) models. According to the initial size of plasma, one may distinguish between the semi-infinite plasma [2, 16–18, 20, 21, 23, 142] and the finite plasma [26–28, 30– 35] in a 1D plasma expansion. These 1D studies showed that there are significant differences between the semi-infinite plasma expansion and the finite plasma expansion. For instance, the self-similar solution to a semi-inifinite plasma expansion in x showed that the potential drop is proportional tox[16, 17], while that to a finite plasma expansion predicted a quadratic dependence of potential decrease along x[27]. Moreover, in a semi-infinite plasma expansion, the ion front velocity exhibits a logarithmic dependence on time and the final ion velocity is indefinite [2, 21]. In contrast, in a finite-size plasma expansion, the ion acceleration is limited [30–32, 34]. The expansion process of a plasma plume into a vacuum happens in two (2D) or three (3D) dimensions, and is qualitatively different from the 1D process. The studies of plasma 5.1 Introduction 73 expansion in dimensions higher than 1D mostly rely on numerical simulations. The PIC method is a conventional tool for kinetic simulations of the collisionless plasma expansion. Many previous studies considered that the plume expansion was governed by the ion dynamics, and the electron dynamics can be approximated by simplified fluid models. Hence, the hybrid PIC approach has been widely used. The fluid treatment of electrons is valid only when the electrons are in an equilibrium or near-equilibrium state. This requires frequent collisions between particles, a condition not satisfied in a collisionless plasma plume. None of the fluid electron models for the hybrid PIC in the literature has considered the effects of non-equilibrium electron characteristics or assessed the validity of using a fluid approach to model electrons. Recently, fully kinetic PIC simulations have been carried out to study the expansion of a plasma plume from an EP thruster [79, 81–84]. To address the validity or the effects of using a fluid treatment of electrons in the expansion of a finite-size collisionless plasma plume, a direct comparison between the hybrid and full PIC simulations is needed. Few studies have carried out such direct comparisons which are very computationally demanding for the fully kinetic simulations. Very recently, Pfeiffer et al. [143] made one of the first such attempts. However, their results are meaningful only within a very short time scale due to the small simulation domain used. This chapter carries out both hybrid PIC and large scale fully kinetic PIC simulations of collisionless plasma plume expansion and provides a direct comparison between the two approaches. Results are presented on the plume dynamics, and the electron and ion characteristics involving the expansion. A quantitative assessment of the effects of a electron fluid model on the collisionless hypersonic plasma plume expansion is provided. As Ar + and Xe + are two commonly used propellent species in EP, the effects of ion-to-electron mass ratio on the expansion are also examined by considering H + , Ar + , and Xe + in the full PIC simulations of this study. 5.2 Simulation Model 74 Table 5.1: Ion-to-electron mass ratio, plasma plume speed, and electron-to-ion plasma frequency ratio at emission surface Model m i =m e v 0 =C s0 v 0 =v te0 ! pe0 =! pi0 Hybrid PIC 1 20 - - Full PIC (H + ) 1836 20 0.46674 42.85 Full PIC (Ar + ) 72820 20 0.074115 269.85 Full PIC (Xe + ) 241073 20 0.040734 491 5.2 Simulation Model The hybrid and full PIC methods are described in Sec. 2.2. The most commonly used Boltzmann electron fluid model given by Eq. (2.20) is used for the hybrid simulation. All the quantities are expressed in the normalized form, and denoted by the notation “ ~”. The reference condition, denoted by the subscript “0” in Eqs. (2.24) and (2.25), used for normalization is chosen to be the emission source condition. At the plume exit surface, the plume speed is taken to be v 0 =C s0 = (20; 0; 0) and the ion temperature is T i0 =T e0 = 0:01. We note that an ion thruster typically operates with the plume Mach number ofM 0 =v 0 =C s0 = 20 30. In the fully kinetic PIC simulations, the real ion-to-electron mass ratios for H + , Ar + and Xe + are used. The real ion-to-electron mass ratios are used in the fully kinetic PIC simulations. Table 5.1 compares some of the mass ratio related parameters at the emission surface used by different simulation runs in this paper. It is noted that m i =m e in the present hybrid PIC which uses the Boltzmann electron fluid model can be considered to be1. Fig. 5.1 illustrates the simulation setup in this study. The simulation domain is a rectangular box with L x L y = 2000 D0 1000 D0 . The emission source is positioned in the region of (~ x; ~ y) = (0 6; 0 20). Particles are injected into the simulation domain through the injection plane located at (~ x; ~ y) = (6; 0 20). The electric potential of the source is ~ s = 0, and a zero-electric field condition is used at the simulation domain boundary. Hence, ~ s is floating with respect to the ambient. Any particles (mainly the electrons in the fully kinetic PIC simulations) attracted back to the emission source are removed from the simulation domain. At the boundaries of ~ x = 0 and 2000, and ~ y = 1000, an absorbing 5.3 Results and Discussions 75 Figure 5.1: Schematics of the simulation setup for a collisionless hypersonic plume expansion into a vacuum. boundary condition for particles is also applied. The simulation is considered to be symmetric with respect to ~ y = 0, so the specular reflection condition for the particles and the symmetric condition for the Poisson’s solver are imposed at ~ y = 0. Because of the symmetric boundary with respect to ~ y = 0, the emission source is considered to have a radius of R 0 = 20 D0 . The simulation domain is discretized with a uniform mesh. The cell size is ~ x = ~ y = 1. The time step adopted to push the particles is ~ t = t! pe0 = 0:05 in all the full PIC simulations, and ~ t = t! pi0 = 0:01 in the hybrid PIC simulation, respectively. The simulations are run sufficiently long (t! pi0 50) such that the expansion of plasma plume is well developed, but are terminated when the plasma front is still far away from the outflow boundary to minimize the boundary effects on the expansion. A uniform particle weight is used so the number of macro-particles in a cell scales with the local plasma number density (~ n). Typically, more than 160 million macro-particles are tracked when the simulations are terminated. 5.3 Results and Discussions 5.3.1 Plasma plume structure The expansion of a hypersonic plasma plume is an ion time scale process. Here, the simulation results up to t! pi = 50 are presented. Longer simulations (up to t! pi0 = 60) were run, but not shown, to ensure the conclusions are still valid at later time. 5.3 Results and Discussions 76 0 200 400 600 800 1000 1200 1400 1600 ˜ x −14 −12 −10 −8 −6 −4 −2 0 2 ˜ Φ (a) ˜ k= −1 30 tω pi =30 ˜ x=30(v0/Cs0−1)+ ˜ xs ˜ k= −1 40 tω pi =40 ˜ x=40(v0/Cs0−1)+ ˜ xs ˜ k= −1 50 tω pi =50 ˜ x=50(v0/Cs0−1)+ ˜ xs H + Ar + Xe + hybrid 0 200 400 600 800 1000 1200 1400 1600 ˜ x −14 −12 −10 −8 −6 −4 −2 0 2 ˜ Φ (b) ˜ k= −1 30 tω pi =30 ˜ x=30(v0/Cs0−1)+ ˜ xs ˜ k= −1 40 tω pi =40 ˜ x=40(v0/Cs0−1)+ ˜ xs ˜ k= −1 50 tω pi =50 ˜ x=50(v0/Cs0−1)+ ˜ xs H + Ar + Xe + hybrid 0 200 400 600 800 1000 1200 1400 1600 ˜ x −14 −12 −10 −8 −6 −4 −2 0 2 ˜ Φ (c) ˜ k= −1 30 tω pi =30 ˜ x=30(v0/Cs0−1)+ ˜ xs ˜ k= −1 40 tω pi =40 ˜ x=40(v0/Cs0−1)+ ˜ xs ˜ k= −1 50 tω pi =50 ˜ x=50(v0/Cs0−1)+ ˜ xs H + Ar + Xe + hybrid Figure 5.2: One-dimensional electric potential profiles along the x-axis (plume direction) inside the core region for t! pi = 30 (blue), 40 (red) and 50 (black). The profiles are extracted along ~ y = (a) 0, (b) 10, and (c) 18. Fig. 5.2 shows the time evolution of the 1D electric potential profiles along the x-axis (plume direction) inside the core plume region after the expansion is already well developed. Four distinct plasma regions from the plasma emission source downstream to the vacuum can be observed: (1) the unperturbed plasma region, (2) the quasi-steady propagation region, (3) the self-similar expansion region, and (4) the pure electron gas region. In the unperturbed region (1), the electric potential and plasma density is almost constant. The process is similar to the expansion process of a supersonic neutral gas flow, as pointed out in Refs.[1, 38]. The boundary separating the unperturbed region from the plasma expansion region should be the first Mach line originating at the edge of plasma emission source exit (~ x; ~ y) = (6; 20). This is evident in the 2D plots shown in Figs. 5.3 and 5.4. Fig. 5.3 shows the potential iso-contours and the electric field vectors, and Fig. 5.4 shows the iso-contour lines of the ion and electron densities. The black dashdot lines overlaid in Figs. 5.3 and 5.4 is the first Mach line. The slope of the line k is computed as k = tan(); where' sin 1 1 M 0 : (5.1) It is clear that, within the region bounded by the plasma emission plane and the first Mach line, there are negligible variations in the potential and the plasma density. In Fig. 5.2(c), this unperturbed region is nearly invisible because the two boundaries are too close to each other. 5.3 Results and Discussions 77 Figure 5.3: (Color) Iso-contours of the electric potential at t! pi = 50 for the full PIC results with (a) H + , (b) Ar + , (c) Xe + , and (d) the hybrid PIC result. Yellow arrows show the electric field vectors. Black dashdot lines (“ ”) are the calculated first Mach lines. The inner white and outer green solid lines correspond to the ion density iso-lines with the values of ~ n i (~ x = 50[~ v 0 =C s0 1] + ~ x s ; ~ y = 0) and' 0, respectively. Next to the unperturbed region behind the Mach line is the quasi-steady expansion region (2), where the plasma plume structure is similar to that of the Prandtl-Meyer expansion fan in a supersonic gas flow. It can be seen from Fig. 5.2 that in this region the potential profiles at the earlier time overlaid with the ones at the later time for all the models. The quasi-steady expansion region ends at a deflection point where the potential starts to decrease rapidly. The lines obtained from the full PIC runs for different ion-to-electron mass ratios 5.3 Results and Discussions 78 (d) 0 200 400 600 800 1000 1200 1400 0 100 200 300 400 Figure 5.4: (Color) Iso-contour lines of the ion (“ ”) and electron (“ ”) number density at t! pi = 50 for the full PIC results with (a) H + , (b) Ar + , (c) Xe + , and (d) the hybrid PIC result. Black dashdot lines (“ ”) are the calculated first Mach lines. are almost indistinguishable in Fig. 5.2. The hybrid PIC model predicts a faster decline in the potential profiles than the full PIC simulations in this region. The quasi-steady expansion region is followed by a region with a linear decrease in potential. It is apparent in Fig. 5.2 that the slope of the potential profile almost precisely equals1=(t! pi ). It is notable that the self-similar solution by Gurevich et al. [16] predicts e=(k b T e ) =1x=(tC s ) or ~ =1 ~ x=(t! pi ). Therefore, this linearly decreasing potential region is recognized as the self-similar expansion region (3). One of the most remarkable features observed is that the deflection point between regions (2) and (3) propagates downstream along the x-axis with a constant speed for all the models considered in this study. In Fig. 5.2, we show, for t! pi = 30; 40 and 50, the vertical dotted lines calculated from ~ x ~ x s =t! pi (v 0 C s0 )=C s0 =t! pi (M 0 1); (5.2) 5.3 Results and Discussions 79 0 200 400 600 800 1000 1200 1400 1600 ˜ x 0.0 0.5 1.0 1.5 2.0 ˜ T ˜ T e,x ˜ T e,y (a) H + Ar + Xe + 0 200 400 600 800 1000 1200 1400 1600 ˜ x 0.0 0.5 1.0 1.5 2.0 ˜ T ˜ T e,x ˜ T e,y (b) H + Ar + Xe + 0 200 400 600 800 1000 1200 1400 1600 ˜ x 0.0 0.5 1.0 1.5 2.0 ˜ T ˜ T e,x ˜ T e,y (c) H + Ar + Xe + Figure 5.5: One-dimensional electron temperature profiles inside the core plume region for t! pi = 50 extracted along ~ y = (a) 0, (b) 10, and (c) 18. where ~ x s = 6 denotes the position of the source injection plane. It is evident that these vertical dotted lines match the boundaries, as indicated by the deflection points on the potential profiles, between the quasi-steady and the self-similar expansion regions. The self-similar solution to the expansion of a 1D semi-infinite collisionless plasma into a vacuum by Gurevich et al. [16] predicts a rarefaction wave propagating inward the unperturbed plasma with the ion acoustic velocity. The dynamics of the moving boundary of the quasi-steady expansion region is determined by the combination of the plasma plume’s downstream propagation and the rarefaction wave’s upstream propagation. Consequently, the propagation speed of the boundary (deflection point) should equal v 0 C s0 instead of v 0 . Gurevich’s self-similar solution is based on the isothermal electron assumption so the ion acoustic velocity does not change in the process. Accordingly, the observed constant propagation speed along the plume (x-) direction implies that the electrons in the quasi-steady expansion region must be isothermal in the plume direction. Fig. 5.5 shows the electron temperature profiles at t! pi0 = 50 extracted at the same locations as those in Fig. 5.2. ~ T e;x and ~ T e;y denote the x- and the y-component of electron temperature, respectively. The electron temperature exhibits a strong anisotropic feature in the plasma plume. We find that ~ T e;x is almost constant (' 1) in the entire unperturbed and quasi-steady expansion regions while ~ T e;y experiences a considerable decrease as soon as the electrons are emitted from the source. The nearly uniform ~ T e;x is consistent with the observed constant propagation speed along the plume streaming (x-) direction of the quasi-steady expansion region’s boundary. 5.3 Results and Discussions 80 10 20 30 40 50 60 tω pi −2.0 −1.5 −1.0 −0.5 0.0 ˜ k ·tω pi H + Ar + Xe + hybrid Figure 5.6: Fitting value of the normalized slope of the electric potential profile in the self-similar expansion region as a function of time. Such a behavior is closely related to the electric field characteristics. The boundary of the quasi-steady expansion region, shown by the inner white solid line in Fig. 5.3, is considered to be the ion density iso-line with the value equal to ~ n i at (~ x; ~ y) = (t! pi [M 0 1]+ ~ x s ; 0). A careful examination reveals that, except for some random noise inherent in particle simulations, the x-component of the electric field, E x , is nearly negligible. On the other hand, there is a significant positive y-component, E y , near the boundary of the unperturbed region. E y decreases with ~ x inside the quasi-steady propagation region. As the plasma expands into the vacuum, the electrons experience no force along the plume (x-) direction and hence no energy loss. In the transverse (y-) direction, the electrons must overcome the electrostatic force , which leads to the temperature decrease in that direction. The transverse electric field component in the quasi-steady expansion region is also responsible for the slow decrease in the potential profiles observed in Fig. 5.2. Fig. 5.3 clearly shows that E x jumps from nearly zero in the quasi-steady expansion region abruptly to a finite value in the self-similar region immediately downstream of the boundary between the two regions. It also shows that E y is nearly negligible for ~ y< 20 and then slowly increases with ~ y for ~ y > 20 in this self-similar expansion region. In response to the electric field, ~ T e;x decreases quickly as the electrons overcome the electrostatic force in the x-direction while ~ T e;y keeps almost constant in the self-similar expansion region, as clearly shown in Fig. 5.5. The self-similar expansion feature of this region, once established, 5.3 Results and Discussions 81 is preserved over time. This is demonstrated by Fig. 5.6, which shows the product of t! pi and the fitted dimensionless slope ~ k of the potential profile in the self-similar expansion region as a function oft! pi . We find that the values of ~ kt! pi are approximately1, with the variation less than 5%, for all the models considered in this study when t! pi 15. The self-similar expansion region ends near the ion front. It is worth mentioning that the self-similar solution to the 1D semi-infinite collisionless plasma expansion problem fails near the ion front due to effects of charge separation [20]. The 2D ion front is displayed with the green solid line in Fig. 5.3, and marked with ~ n i = 0 in Fig. 5.4. The region beyond the ion front contains only electrons of very low density, as shown in Fig. 5.4. Previous 1D simulations [2, 20, 21] have shown that the strength of electric field achieves a maximum near the ion front. The present simulations show that, in general, such an observation still holds in the 2D expansion, despite of the much more complex 2D electric field structures. The electrostatic force at the ion front drives the propagation of the ion front. 5.3.2 Ion dynamics In a plasma plume expansion, the overall plume dynamics and plume structure are mainly determined by the ion dynamics. Fig. 5.7 shows the ion velocity contours at t! pi = 50 for all the models considered in this study. In contrast to the 1D plasma expansion where only ion acceleration is possible, the 2D plasma expansion exhibits a more complex ion velocity field. Fig. 5.7 (left panel) shows that both acceleration and deceleration regions exist for the x-component of ion velocity, v i;x . The acceleration region (v i;x =C s0 >M 0 + 1%~ v max i;x ) and deceleration region (v i;x =C s0 <M 0 0:2%~ v max i;x ) are illustrated with the white solid lines as the borders on top of the v i;x contour. Here, ~ v max i;x is defined as v max i;x =C s0 M 0 , the velocity increase in the x-direction of the fastest ions. In the plume core region, the ions are neither accelerated nor decelerated along the x-direction. Substantial accelerations for the ions are seen beyond the plume core region. It is interesting to note that the ions are decelerated in region outside of the core region in the transverse direction. The acceleration and deceleration 5.3 Results and Discussions 82 Figure 5.7: Ion velocity contours at t! pi = 50 for the full PIC results with (a) H + , (b) Ar + , (c)Xe + , and (d) the hybrid PIC result. Thex-component and they-component ion velocities normalized by C s0 are shown on the left panel (1) and the right panel (2), respectively. regions are consistent with the electric field shown in Fig. 5.3. The electric field has a large positive E x ^ x in the acceleration region, and a negative E x ^ x in the deceleration region. Fig. 5.7 (right panel) shows that the y-component ion velocity (v i;y ) is accelerated immediately outside the unperturbed region. v i;y experiences no deceleration. The boundary of the acceleration region, shown by the white solid line in Fig. 5.7 (right panel), is set to be v i;y =C s0 = %1~ v max i;y , where ~ v max i;y =v max i;y =C s0 is the velocity increase in the y-direction of the fastest ions. This line approximately agrees with the boundary of the unperturbed region. One can see from Fig. 5.7 that both v max i;x and v max i;y in the hybrid simulation result are apparently larger than those in the full PIC results. The difference among the full PIC simulations with various ion masses is almost invisible in the 2D contours. Therefore, a careful examination of the model effects on the ion acceleration process is necessary. Fig. 5.8 shows the velocity increase of the fastest ions varying with time for t! pi = 5 50. In Fig. 5.8, the symbols denote the simulation results from different models used in this 5.3 Results and Discussions 83 0 10 20 30 40 50 tω pi 0 1 2 3 4 5 6 7 8 v max i,x /C s0 −M 0 8.215− 90.465 tω pi +11.88 8.338− 88.694 tω pi +11.21 8.418− 94.407 tω pi +12.27 (a) H + Ar + Xe + hybrid Mora [10] 0 10 20 30 40 50 tω pi 0 1 2 3 4 5 6 7 8 v max i,y /C s0 7.365− 106.775 tω pi +17.11 7.616− 107.698 tω pi +16.23 7.752− 112.906 tω pi +16.76 (b) H + Ar + Xe + hybrid Mora [10] Figure 5.8: Velocity increase of fastest ions as a function of time: (a)v max i;x =C s0 M 0 = ~ v max i;x , (b) v max i;y =C s0 = ~ v max i;y . The present simulation results from different models are represented by symbols. The red solid line is the prediction of Eq. (5.3) from Mora’s work [2]. The dotted lines are the curves best fitted to the full PIC simulation results with different ion species according to Eq. (5.4). study. For both v max i;x (Fig. 5.8(a)) and v max i;y (Fig. 5.8(b)), the hybrid model predicts higher velocity increases than all the three full PIC models. The values of ~ v max i;x and ~ v max i;y are similar to each other for all the time in the hybrid PIC simulation. However, the full PIC simulations predict a considerably larger ion velocity increase in the x-direction than that in the y-direction. Furthermore, different ion masses also affect the ion accelerations in the full PIC simulations. A higher ion mass yields a larger amount of acceleration in both v max i;x and v max i;y . Also presented by the red solid line in Fig. 5.8 is the time dependence of ion front velocity from Mora’s semi-analytical formula [2], v i =C s0 = 2 ln + p 2 + 1 ; (5.3) where = t! pi = p 2e N and the constant e N = 2:718 28::: is the Euler’s number. We find that the prediction of Eq. (5.3) agrees well with the present hybrid PIC results, even though Eq. (5.3) is derived from the expansion of a 1D semi-infinite plasma while the present simulation is for a 2D finite-size plume. Nevertheless, such an agreement is not surprising since the derivation of Eq. (5.3) is based on the same assumption that the electrons can be 5.3 Results and Discussions 84 Table 5.2: Best fitting parameters to ion velocity increase in the x-direction Eq. (5.4a). Model (ion species) p 0 p 1 p 2 RMSE(~ v max i;x ) H + 8.215 90.465 11.88 0.95% Ar + 8.338 88.694 11.21 0.42% Xe + 8.418 94.407 12.27 0.97% Table 5.3: Best fitting parameters to ion velocity increase in the y-direction Eq. (5.4b). Model (ion species) q 0 q 1 q 2 RMSE(~ v max i;y ) H + 7.365 106.775 17.11 2.73% Ar + 7.616 107.698 16.23 2.73% Xe + 7.752 112.906 16.76 2.34% considered as a massless, isothermal fluid and described by the Boltzmann relation Eq. (2.20). This assumption gives the electrons an infinite amount of energy to drive the plasma expansion. Therefore, the hybrid model tends to overestimate the ion acceleration, as shown in Fig. 5.8. When the electron kinetics is fully taken into account, no analytical expression for the dependence of the ion acceleration is available in the literature. We fit the present full PIC data of ion velocity increases in the x- and y-directions to the expression given as follows, v max i;x (t)=C s0 M 0 = ~ v max i;x =p 0 p 1 t! pi +p 2 ; (5.4a) v max i;y (t)=C s0 = ~ v max i;y =q 0 q 1 t! pi +q 2 ; (5.4b) wherep 0;1;2 andq 0;1;2 are the fitting parameters to be determined. The best fitting parameters to Eqs. (5.4a) and (5.4b) are listed in Tables 5.2 and 5.3, respectively. RMSE in the tables denotes the root mean square error defined as, RMSE(f) = s P n t=1 ( ^ f t f t ) 2 n : (5.5) where f represents either ~ v max i;x or ~ v max i;y , and ^ f t denotes the predicted value by the fitting formula with respect to the t th data point f t from the simulation results. The error of the fitting expression given in Eq. (5.4) is less than 1% and 3% for the acceleration of the fastest 5.3 Results and Discussions 85 42.85 269.85 491 7 7.5 8 8.5 9 Maximum Ion Velocity Increase Figure 5.9: Dependence of final velocity increase of the fastest ions on mass ratio. ions in the plume and transverse directions, respectively. The best fitting curves of ~ v max i;x and ~ v max i;y are also presented in Fig. 5.8 by the dotted lines. The final velocity gains of the fastest ions can be extrapolated from Eq. (5.4), that is, v max i;x (t!1)=C s0 M 0 =p 0 andv max i;y (t!1) =q 0 . Fig. 5.9 shows the ion mass dependence of the final velocity gain of the fastest ions. It is found that both ~ v max i;x and ~ v max i;y scale with p m i =m e , but ~ v max i;y shows a stronger dependence on p m i =m e . Denavit [21] showed that the electron cooling was inversely proportional to p m i =m e . Therefore, a larger value of m i =m e corresponds to less effective electron cooling, and thus a higher electron temperature. This can be clearly seen from the ~ T e;y profiles in Fig. 5.5. The energy to drive the ion front motion comes from the electron thermal energy. A plume with heavier ions tends to preserve a higher electron temperature. This allows the ions to be accelerated to a higher final velocity. Finally, the observation that ~ v max i;x shows less dependence on the ion mass than ~ v max i;y can be attributed to the fact that the electron cooling is more effective in the y-direction. The reason is that, as shown in Fig. 5.5, ~ T e;x is almost constant until the self-similar expansion region while ~ T e;y decreases throughout most of the plume region. 5.3 Results and Discussions 86 5.3.3 Effects of electron characteristics The full PIC results are directly compared to the results obtained from the hybrid simulation usingtheBoltzmannelectronfluidmodel. Thecomparisonshowsthat, besidesthequantitative difference between the hybrid and fully kinetic results, the oversimplified electron fluid model is not able to capture the different plasma expansion dynamics in both the plume and the transverse directions. This qualitative difference is attributed to the strong anisotropy in electron temperature that is omitted by the hybrid PIC modeling. Past studies have shown that the 1D plasma expansion may be classified as the semi- infinite plasma expansion or the finite plasma expansion. For a freely expanding plasma in 1D, there exists a rarefaction wave propagating inward the plasma with the ion acoustic speed of C s . Therefore, for a time scale t, a semi-infinite plasma refers to a plasma with a size L > C s t such that the rarefaction wave front cannot reach the plasma center and an unperturbed plasma region always exists, while a finite plasma has a size L<C s t and thus the entire plasma is affacted by the expansion. Hence, for the 2D hypersonic plasma plume considered in this study, the expansion along the plume direction can be considered similar to a 1D semi-infinite plasma expansion, while that along the transverse expansion is similar to a 1D finite plasma expansion. The main difference between the semi-infinite plasma and finite plasma expansion is the completely different electron thermodynamics involved. The semi-infinite plasma serves as a reservoir with an unlimited amount of energy such that the change of electron temperature may be ignored. Consequently, the isothermal electron fluid model is a good approximation for the expansion of a semi-infinite plasma into a vacuum. The effect of electron cooling was found to be proportional to (m i =m e ) 1=2 [21, 23]. In particular, Denavit’s [21] kinetic simulation showed that the cooling factor = 1, where denotes the coefficient in a polytropic relation between the temperature and density (or pressure), was as small as 0:011 when m i =m e = 1600. These conclusions justified the use of isothermal electron model in the semi-infinite plasma expansion problems as m i =m e is at least 1836:15 so electron cooling is negligible. For a finite plasma expansion, assuming the electrons are isothermal 5.3 Results and Discussions 87 0 200 400 600 800 1000 1200 1400 1600 ˜ x 0.0 0.5 1.0 1.5 2.0 ˜ T e ˜ T e,x ˜ T e,y R beam =20λ D0 R beam =40λ D0 Figure 5.10: One-dimentional electron temperature profiles of H + beams for beam sizes R 0 = 20 D0 and 40 D0 at t! pi = 50. The profiles are extracted along the beam direction at the position 20 D0 away from the beam edge, namely along ~ y = 0 (beam center) for R 0 = 20 D0 and ~ y = 20 for R 0 = 40 D0 , respectively. is inappropriate because only a finite amount of electron energy is available for driving the expansion. It is pointed out that the electron temperature can decrease as fast as proportional to t 2 in a finite plasma if no external energy is supplied [26, 28, 30]. Consequently, electron cooling plays an important role in the expansion of a finite plasma. The full PIC results shown in Fig. 5.5, which shows a near constant ~ T e;x and a decrease in ~ T e;y , supports the above argument. The distinct ion acceleration processes along the plume and transverse directions, e.g. smaller ~ v max i;y than ~ v max i;x , seen in the full PIC simulations are due to such anisotropic characteristics of electron temperatures. According to our hypothesis, it is expected that the electron cooling along the transverse direction should be less effective when the beam size is increased. To further demonstrate this, we have carried out another full PIC simulation for the H + plume with the beam size doubled (R 0 = 40 D0 ). All other simulation setups are the same as those in Sec. 5.2. Fig. 5.10 shows the comparison of the 1D electron temperature profiles along x-axis between the two different beam sizes att! pi = 50. The comparison has clearly shown that the two ~ T e;x profiles are almost identical if we ignore the statistical noise, but ~ T e;y for R beam = 40 D0 is higher and decreases slower than that for R beam = 20 D0 . 5.4 Summary 88 These observations suggest that, in order for an electron fluid model to be valid in modeling the plasma plume expansion, one needs to at least take the strong anisotropy of electron thermodynamics into account, namely the electron temperature in the plume direction and the transverse direction needs to be modeled separately. This is needed not only for the simplest Boltzmann electron model, but also for the more sophisticated approach of solving the electron energy equation. 5.4 Summary The expansion of a collisionless, hypersonic plasma plume into a vacuum is simulated using both the hybrid and fully kinetic PIC models. Both the hybrid and full PIC simulations reveal that the plume structure exhibits four distinct expansion regions. These regions are identified as, from the plasma emission plane to the downtream along the plume direction, the unperturbed, quasi-steady expansion, self-similar expansion and electron front regions, respectively. The formation of the unperturbed region, bounded by the Mach cone, is similar to the process of a supersonic neutral gas expansion. The establishment of the quasi-steady propagation region is due to the combined effect of the forward plume propagation and the backward rarefaction wave propagation. The self-similar expansion mechanism in the hypersonic plasma plume is recognized the same as that in the 1D semi-infinite plasma expansion. The electron temperatures are shown to be highly anisotropic in the plume. T e;x is almost unchanged from the source until the self-similar region, while T e;y drops drastically immediately downstream of the injection plane. A consequence of the anisotropic electron temperatures is the different accelerations that the fastest ions can achieve. A higher electron temperature generally results in a larger ion acceleration. Besides, a smaller ion mass leads to a more rapid electron temperature decrease and thus a lower electron temperature in the transverse direction. The saturated accelerations of the fastest ions are obtained by 5.4 Summary 89 extrapolating the present data to the time at infinity. The extrapolated saturated ion accelerations scale with p m i =m e . The anisotropic electron thermodynamics is attributed to the different expansion features along the plume direction and the transverse direction. On the one hand, the plasma expansion in the plume direction is similar to the 1D semi-infinite plasma model, so the electrons are able to maintain isothermal in that direction. On the other, the expansion process transverse to the plume can be analogous to the 1D finite plamsa situation, where the thermal energy is limited and thereby the electron cooling effect on the expansion is important. A direct comparison between the hybrid and full PIC simulations demonstrates that the commonly used assumption that the electrons may be considered as an isothermal fluid is over-simplified and causes not only quantitative but also qualitative differences in modeling the collisionless, hypersonic plasma expansion. For an electron fluid model to be applicable, the model must at least have the capability of dealing with the strong anisotropic properties of thermodynamics. 90 Chapter 6: Expansion of a Collisionless Hypersonic Plume into Vacuum: 3D Effects In the last chapter, a two-dimensional (2D) collisionless hypersonic plasma plume is study using both both the hybrid and fully kinetic particle-in-cell (PIC) simulations. It is found that the 2D expansion of a collisionless plasma plume is qualitatively different from those well stduied 1D expansion models. This chapter presents the simulations of an axi-symmetric plume and investigates the 3D effects on the plasma expansion dynamics. 6.1 Introduction The expansion of a collisionless plasma into vacuum is a very generic problem that arises in many applications relevant to plasma physics. This problem has been extensively studied with various 1D models theoretically in the literature. For real applications, which primarily happen in 2D or 3D, numerical simulations are the main tools. The last chapter presents a comprehensive study of a 2D collisionless hypersonic plasma plume expanding into a vacuum. The simulations are set to mimic the real plasma jets from an electric propulsion (EP) thruster. It is generally considered that the differences between 1D and 2D flows are qualitative. Hence, the results obtained from the theories derived based on the 1D models cannot be directly exntended to the situations in higher dimensions without modifications. The simulations presented in the last chapter shows that the electron kinetic effects are significant and exhibits strong ansotropic characteristics in the 2D collisionless plasma expansion, qualitatively differing from those well established theory based on 1D models. Moreover, our investigation reveals that the expansions of a 2D collisionless hypersonic plasma plume into a vacuum in 6.2 Simulation Model 91 Figure 6.1: Schematics of the simulation setup for a collisionless hypersonic plume expansion into a vacuum. its plume direction and transverse direction are closely related to two well studied 1D models, namely, the expansions of a semi-infinite plasma and a finite plasma, respectively. The expansion of an EP thruster plume is indeed a 3D or axi-symmetric process, though it is generally considered that the results from 3D flow models are only quantitatively different from those 2D results. This chapter presents a comparative study of the collisionless hypersonic plasma plume expansion between the 2D and 3D (axi-symmetric) models. The 3D effects on the electron kinetics as well as the plasma plume expansion dynamics are the primary focus of this study. 6.2 Simulation Model Both the hybrid PIC and fully kinetic PIC simulations are carried out. Same as the last chapter, the hybrid PIC modeling employs the Boltzmann electron fluid model Eq. (2.20). For the full particle PIC simulations, only Xe + is considered because it is the the most commonly used propellant of an ion thruster and the hybrid PIC results agree better with the results based on Xe + compared with other lighter ions. The normalization for the hybrid and full PIC modeling follow the same strategies presented in the last chapter. The plume conditions is the same as those in the last chapter, listed in Table 6.1. The geometry of the simulation domain and other simulation setups, shown in Fig. 6.1, are the same as those in the last chapter. Note that ~ r = 0 is the axis of symmetry for 6.2 Simulation Model 92 Table 6.1: Ion-to-electron mass ratio, plasma plume speed, and electron-to-ion plasma frequency ratio at emission surface Model m i =m e v 0 =C s0 v 0 =v te0 ! pe0 =! pi0 Hybrid PIC 1 20 - - Full PIC (Xe + ) 241073 20 0.040734 491 the axi-symmetric simulations. For an easy reading, the detailed simulation parameters are summarized again here. The simulation domain is a rectangular box with L x L r = 2000 D0 1000 D0 . The emission source is positioned in the region of (~ x; ~ r) = (0 6; 0 20). Particles are injected into the simulation domain through the injection plane located at (~ x; ~ r) = (6; 0 20). The electric potential of the source is ~ s = 0, and a zero-electric field condition is used at the simulation domain boundary. Hence, ~ s is floating with respect to the ambient. Any particles (mainly the electrons in the fully kinetic PIC simulations) attracted back to the emission source are removed from the simulation domain. At the boundaries of ~ x = 0 and 2000, and ~ r = 1000, an absorbing boundary condition for particles is also applied. The simulation is considered to be symmetric with respect to ~ r = 0, so the specular reflection condition for the particles and the symmetric condition for the Poisson’s solver are imposed at ~ r = 0. Because of the symmetric boundary with respect to ~ r = 0, the emission source is considered to have a radius of R 0 = 20 D0 . The simulation domain is discretized with a uniform mesh. The cell size is ~ x = ~ r = 1. The time step adopted to push the particles is ~ t = t! pe0 = 0:05 in all the full PIC simulations, and ~ t = t! pi0 = 0:01 in the hybrid PIC simulation, respectively. The simulations are run sufficiently long (t! pi0 50) such that the expansion of plasma plume is well developed, but are terminated when the plasma front is still far away from the outflow boundary to minimize the boundary effects on the expansion. One notable difference between the 2D and axi-symmetric simulations is the distribution of particles. For 2D simulations, the number of particles is distributed uniformly along the transverse or radial (y-) direction if the plasma density keeps the same in that direction. A 6.3 Results and Discussions 93 uniform particle weight is used so the number of macro-particles in a cell scales with the local plasma number density (~ n). For axi-symmetric simulations, the number of particles, instead, increases quadratically along the radial direction as moving away from the center because the volume is proportional to r 2 . A consequence of the quadratic increase of particle number in the radial direction is that a relatively larger number of macro-particles has to be used. The reason is if we want to keep the same number of particle per cell in the cells near the axis of symmetric to maintain a low statistical noise, the total number of particles the entire simulation end up with will be large. A remedy for this is to use a non-uniform particle weight so that the number of macro-particles per cell can be maintained almost the same across different r values. A drawback of the non-uniform particle weight scheme is that the split or merge of macro-particles when they cross cells will introduce additional errors, even with the state-of-art strategies [144, 145]. Typically, over 160 and 360 million macro-particles are tracked at the end of the 2D and axi-symmetric plume simulations, respectively. 6.3 Results and Discussions 6.3.1 Evolution of plasma plume The simulation results up to t! pi = 50, the same as that in the last chapter, are presented. Longer simulations (up to t! pi0 = 60) were run, but not shown, to ensure the conclusions are still valid at later time. Fig. 5.2 shows the time evolution of the 1D electric potential profiles along the x-axis (plume direction) inside the core plume region at t! pi = 30, 40 and 50. The comparisons between the 2D and axi-symmetric plumes from both the hybrid and full PIC simulations are presented. It is interesting to see that the same evolution trend as well as the plume structures, namely, the four distinct plasma regions discussed in detail for the 2D plume in the last chapter exist for the axi-symmetric plume. 6.3 Results and Discussions 94 0 200 400 600 800 1000 1200 1400 1600 ˜ x −14 −12 −10 −8 −6 −4 −2 0 2 4 ˜ Φ ˜ k= −1 30 tω pi =30 ˜ x=30(v0/Cs0−1)+ ˜ xs ˜ k= −1 40 tω pi =40 ˜ x=40(v0/Cs0−1)+ ˜ xs ˜ k= −1 50 tω pi =50 ˜ x=50(v0/Cs0−1)+ ˜ xs (a) ˜ y=0 two-dim axi-sym 0 200 400 600 800 1000 1200 1400 1600 ˜ x −14 −12 −10 −8 −6 −4 −2 0 2 4 ˜ Φ ˜ k= −1 30 tω pi =30 ˜ x=30(v0/Cs0−1)+ ˜ xs ˜ k= −1 40 tω pi =40 ˜ x=40(v0/Cs0−1)+ ˜ xs ˜ k= −1 50 tω pi =50 ˜ x=50(v0/Cs0−1)+ ˜ xs (b) ˜ y=10 two-dim axi-sym 0 200 400 600 800 1000 1200 1400 1600 ˜ x −14 −12 −10 −8 −6 −4 −2 0 2 4 ˜ Φ ˜ k= −1 30 tω pi =30 ˜ x=30(v0/Cs0−1)+ ˜ xs ˜ k= −1 40 tω pi =40 ˜ x=40(v0/Cs0−1)+ ˜ xs ˜ k= −1 50 tω pi =50 ˜ x=50(v0/Cs0−1)+ ˜ xs (c) ˜ y=18 two-dim axi-sym 0 200 400 600 800 1000 1200 1400 1600 ˜ x −14 −12 −10 −8 −6 −4 −2 0 2 4 ˜ Φ ˜ k= −1 30 tω pi =30 ˜ x=30(v0/Cs0−1)+ ˜ xs ˜ k= −1 40 tω pi =40 ˜ x=40(v0/Cs0−1)+ ˜ xs ˜ k= −1 50 tω pi =50 ˜ x=50(v0/Cs0−1)+ ˜ xs (d) ˜ y=0 two-dim axi-sym 0 200 400 600 800 1000 1200 1400 1600 ˜ x −14 −12 −10 −8 −6 −4 −2 0 2 4 ˜ Φ ˜ k= −1 30 tω pi =30 ˜ x=30(v0/Cs0−1)+ ˜ xs ˜ k= −1 40 tω pi =40 ˜ x=40(v0/Cs0−1)+ ˜ xs ˜ k= −1 50 tω pi =50 ˜ x=50(v0/Cs0−1)+ ˜ xs (e) ˜ y=10 two-dim axi-sym 0 200 400 600 800 1000 1200 1400 1600 ˜ x −14 −12 −10 −8 −6 −4 −2 0 2 4 ˜ Φ ˜ k= −1 30 tω pi =30 ˜ x=30(v0/Cs0−1)+ ˜ xs ˜ k= −1 40 tω pi =40 ˜ x=40(v0/Cs0−1)+ ˜ xs ˜ k= −1 50 tω pi =50 ˜ x=50(v0/Cs0−1)+ ˜ xs (f) ˜ y=18 two-dim axi-sym Figure 6.2: One-dimensional electric potential profiles along the x-axis (plume direction) inside the core region for t! pi = 30 (blue), 40 (red) and 50 (black). The upper (a, b, c) and lower (d, e, f) panles show the results from the hybrid PIC and full PIC simulations, respectively. The profiles are extracted along ~ r = (a, d) 0, (b, e) 10, and (c, f) 18. In the unperturbed region (1), the electric potential is almost constant for both the 2D and axi-symmetric plumes. Followed by the unperturbed region is the quasi-steady expansion region (2), where the potential decreases along the plume direction with a moderate rate. The boundary between the quasi-steady region (2) and the self-similar expansion region (3) propagates downstream with a constant speed. It is surprising to see that propagation speed of the boundary between region (2) and (3) for the axi-symmetric plume is also v 0 C s0 , the same as that for the 2D plume. The comparison between the 2D and axi-symmetric plumes in the quasi-steady expansion region shows that the potential drops with a faster rate in the axi-symmetric plume. The faster potential decrease in the axi-symmetric plume is signified by the hybrid PIC modeling. The 2D and axi-symmetric plume have almost the identical pattern of the self-similar expansion region. The slope of the potential drop in the this region is well described by the theory based on the 1D semi-infinite plasma model. Though the self-similar expansion characteristics is the same for both the 2D and axi-symmetric plumes, the self-similar expansion region ends or the pure electron region (4) appears much earlier 6.3 Results and Discussions 95 Figure 6.3: (Color) Iso-contour lines of the ion (“ ”) and electron (“ ”) number density at the steady state: (a) two-dimensional hybrid PIC, (b) axi-symmetric hybrid PIC, (c) two-dimensional full PIC and (d) axi-symmetric full PIC. in the latter one. A consequence of this is that the axi-symmetric case has a smaller total potential drop across the entire domain compared with the 2D plume. 6.3.2 Steady state plasma plume In real EP applications, people may be more interested in the steady state operation not far away from the thruster. The plume evolution shown in Fig. 6.2 has provided us with solid evidence about when and where the plume reaches a steady state. For instance, the region with ~ x 800 is guaranteed to reach a steady state after t! pi 50. Fig. 6.3 and Fig. 6.4 show the steady state (t! pi = 50) iso-density lines and 1D profiles along plume direction, respectively. The 2D and axi-symmetric plume structures from both the hybrid and fully kinetic PIC simulations are compared. It is clear that, in the steady state region, the plume is quasi-neutral, namely ~ n i ' ~ n e , for the models. The pattern of the 2D plasma plume agrees qualitatively with the axi-symmetric plasma plume, despite that the density in the 2D plasma decreases with a slower rate outside the unperturbed region. Such 6.3 Results and Discussions 96 0 100 200 300 400 500 600 700 800 10 -2 10 -1 10 0 hybrid PIC (two-dim) hybrid PIC (axi-sym) full PIC (two-dim) full PIC (axi-sym) 0 100 200 300 400 500 600 700 800 10 -2 10 -1 10 0 hybrid PIC (two-dim) hybrid PIC (axi-sym) full PIC (two-dim) full PIC (axi-sym) 0 100 200 300 400 500 600 700 800 10 -2 10 -1 10 0 hybrid PIC (two-dim) hybrid PIC (axi-sym) full PIC (two-dim) full PIC (axi-sym) 0 100 200 300 400 500 600 700 800 10 -2 10 -1 10 0 hybrid PIC (two-dim) hybrid PIC (axi-sym) full PIC (two-dim) full PIC (axi-sym) 0 100 200 300 400 500 600 700 800 10 -2 10 -1 10 0 hybrid PIC (two-dim) hybrid PIC (axi-sym) full PIC (two-dim) full PIC (axi-sym) 0 100 200 300 400 500 600 700 800 10 -2 10 -1 10 0 hybrid PIC (two-dim) hybrid PIC (axi-sym) full PIC (two-dim) full PIC (axi-sym) Figure 6.4: One-dimensional plasma density profiles along the x-axis (plume direction) inside the core region for the steady state results, a comparison between the two-dimensional and axi-symmetric results from both the hybrid and full PIC simulations. The upper (a, b, c) and lower (d, e, f) panles show the ion and electron densities, respectively. The profiles are extracted along ~ r = (a, d) 0, (b, e) 10, and (c, f) 18. 0 100 200 300 400 500 600 700 800 0 20 40 60 80 100 hybrid PIC (two-dim) hybrid PIC (axi-sym) full PIC (two-dim) full PIC (axi-sym) Figure 6.5: Ion streamlines originated from the injection plane ~ x = 6 for ~ r = 2, 10, and 18, a comparison between the two-dimensional and axi-symmetric results from both the hybrid and full PIC simulations. quantitative difference between the 2D and axi-symmetric plumes is more significant in the hybrid PIC modeling than in the fully kinetic PIC modeling. The difference in the density profiles is also consistent with the observation on the potential profiles shown in Fig. 6.2. Fig. 6.5 compares the ion streamlines for different models. The ion streamlines are originated from the plasma injection plane. Three different positions, ~ r = 2, 10, and 18, along the injection plane are shown in Fig. 6.5. In general, the ion trajectories are deflected more when the origin is further away from the center of the plume. The hybrid PIC simulations predict much larger divergence angle for the ion streamlines than the full PIC simulations. In terms of the 3D effects, the ions are deflected more in the axi-symmetric plumes. A larger 6.3 Results and Discussions 97 Figure 6.6: (Color) Iso-contour lines of the ion (“ ”) and electron (“ ”) number density at the steady state: (a) two-dimensional hybrid PIC, (b) axi-symmetric hybrid PIC, (c) two-dimensional full PIC and (d) axi-symmetric full PIC. 0 200 400 600 800 0 0.5 1 1.5 2 two-dim axi-sym 0 200 400 600 800 0 0.5 1 1.5 2 two-dim axi-sym 0 200 400 600 800 0 0.5 1 1.5 2 two-dim axi-sym Figure 6.7: One-dimensional electron temperature profiles along the x-axis (plume direction) inside the core region for the steady state full PIC results, a comparison between the two- dimensional and axi-symmetric results. The profiles are extracted along ~ r = (a) 0, (b) 10, and (c) 18. divergence is expected in the axi-symmetric plume because the density drop is more rapid and the mass should be conserved. The steady state electron temperature contours are shown in Fig. 6.6. The ion streamlines originated from several representative locations are also overlaid. ~ T e;x is observed to maintain the source temperature inside the core plume region (~ r . ~ R 0 ) as the plume propagates downstream, for not only the 2D but also the axi-symmetric case. This also explains why the boundary between region (2) and (3) propagates with a constant speed v 0 C s0 , which 6.3 Results and Discussions 98 requires an isothermal condition to be valid, seen in Fig. 6.2. ~ T e;y inside the plume core region is observed to decrease quickly as the plasma propagates downstream. Fig. 6.7 shows a comparison of the 1D electron temperature profiles extracted along the plume direction inside the core region between the 2D and axi-symmetric cases. It is surprising to see a negligible difference in the electron temperature profiles inside the plume core region between the 2D and axi-symmetric cases. The electron temperatures outside the plume core region exhibits some abnormal char- acteristics. Instead of cooling, the temperature of electron may undergo an increase in the expansion region as its density decreases. In particular, both ~ T e;x and ~ T e;y are considerably enhanced in the region outside the outermost ion streamline. The abnormal electron temper- ature enhancement is mainly due to the kinetic behaviors of the electrons in the low density region, similar to what has been discussed in the plasma wake of plate discussed in Chapter 3. 6.3.3 Electron thermodynamics The comparison between the hybrid and full PIC results have provided sound evidence that the electron fluid appproximation is inadequate for studying the expansion of a collisionless plasma plume into a vacuum, especially in the region outside the plume core region where the electron kinetic effects are important and abnormal electron thermodynamic properties are seen. However, the fully kinetic approaches, though very accurate, are much more computationlly expensive than the fluid based methods. As a result, the fluid based modelings of the plasma plume expansion, including full fluid and hybrid fluid-kinetic methods, are still heavily used. A popular method is to employ the polytropic law, Eq. 1.14, to describe the electron thermodynamics mainly due to its simplicity but better generalizability than the isothermal assumption. The proper choice of the polytropic coefficient is the key for accurate modelings. However, lacking reliable high fildelity data, the selection of is arbitary in practice [75, 76, 114]. Fig. 6.8 shows the log 10 ~ T e log 10 ~ n e relation along several representaive ion streamlines for both the 2D and axi-symmetric plumes in the steady state region. Note that the polytropic 6.3 Results and Discussions 99 -1.5 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 two-dim axi-sym -1.5 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 two-dim axi-sym -1.5 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 two-dim axi-sym -1.5 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 two-dim axi-sym -1.5 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 two-dim axi-sym -1.5 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 two-dim axi-sym -1.5 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 two-dim axi-sym -1.5 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 two-dim axi-sym -1.5 -1.25 -1 -0.75 -0.5 -0.25 0 0.25 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 two-dim axi-sym Figure 6.8: log 10 ( ~ T e ) log 10 (~ n e ) along the ion streamlines originated from the injection plane ~ x = 6 at various ~ r positions: (a) 2, (b) 4, (c) 6, (d) 8, (e) 10, (f) 12, (g) 14, (h) 16, and (i) 18. coefficient can be obtained from the slopes (k) of the lines on the log 10 ~ T e log 10 ~ n e plot, that is = + 1 = k + 1. The polytropic relations is plotted along the ion streamlines because we consider the ion dynamics is the only concern in a fluid bases modeling. First of all, the strong anisotropic themodynamic properties of electrons are clearly seen for both the 2D and axi-symmetric plumes. Along the ion streamline whose origin ~ r is within 12, ( ~ T e;x ~ n e ) is approximately 0, indicating isothermal expansions. The log 10 ~ T e;y log 10 ~ n e plots show almost linearly relations along the same set of ion streamlines for both the 2D and axi-symmetric cases. The data is fitted to a linear relation to get the slope (k) of each set of log 10 ~ T e;y log 10 ~ n e , and thus the polytropic coefficient ( ~ T e;y ~ n e ). It is found that ( ~ T e;y ~ n e ) slightly increases as the origin of the ion streamline moves away from the center. Overall, ( ~ T e;y ~ n e )' 2.1 and 1.6 for the 2D and axi-symmetric plumes, respectively. It is noted that, for a monotonic gas like electron, = 2 and 5/3 for the 2D and 3D adiabatic 6.4 Summary 100 processes, respectively. ( ~ T e;y ~ n e ) obtained in the present full PIC simulations is close to the values for adiabatic processes when the ion streamlines are close to the plume center. However, it should be emphasized that this does not mean the expansion process is an adiabatic, or even a hydrodynamic process due to the strong anisotropic electron kinetics. In addition, for the ion streamlines that are closer to the edge of the plasma source (~ r 14), the simple polytropic relation completely fails to describe the electron thermodynamics, as shown in Figs. 6.8 (g)(i). 6.4 Summary This chapter presents a comparative study of the expansion of a collisionless hypersonic plasma plume into a vacuum in 2D and 3D. Both the hybrid and fully kinetic PIC simulations are carried out. The comparison shows that the 2D plume is qualitatively similar to the 3D plume, but quantitative differences exist. The plasma density drops more rapidly in the 3D plume expansion than in the 2D plume. The more rapid decrease in the plasma density profile in the 3D plume is signified by the hybrid PIC modeling. The electrons exhibit strong anisotropic behaviors in both the 2D and 3D plume expansions. ~ T e;x is almost constant and equal to the source temperature inside the plume core region, but shows abnormal enhancement in the outer expansion region. ~ T e;y decreases as the plasma expands downstream inside the plume core region, but is elevated unusually in the outside region as well. However, there are small diferences in the spatial variations of both ~ T e;x and ~ T e;y inside the plume core region between the 2D and 3D plumes. As a result, the smaller ( ~ T e;y ~ n e ) for the 3D plume is due to the more rapid decrease of the plasma density. The fully kinetic PIC results show that the fluid models are generally inadequate for the collisionless plasma plume expansion because the electron kinetic behaviors are mainly non-equilibrium, especially in the outer expansion region. 101 Chapter 7: Fully Kinetic Simulation of Collisionless, Mesothermal Plasma Emission: Macroscopic Plume Structure and Microscopic Electron Characteristics This chapter presents a fully kinetic particle-in-cell (PIC) simulation study on the emission of a two-dimensional collisionless plasma plume consisting of cold beam ions and thermal electrons. All the important processes, including the ion-electron coupling, beam neutralization, and plume expansion, relevant to the beam plasma emission from an ion thruster are taken into account in the present model. The main focus is to understand the connections between the macroscopic plume structure and the microscopic electron kinetics in a two-dimensional analysis. It is found that the macroscopic plume structure exhibits several distinctive regions, including an undisturbed core region, an electron cooling expansion region, and an electron isothermal expansion region. The properties of each region are determined by microscopic electron kinetic characteristics. The division between the undisturbed region and the cooling expansion region approximately matches the Mach line generated at the edge of the emission surface, and that between the cooling expansion region and the isothermal expansion region approximately matches the potential well established in the beam. The interactions between electrons and the potential well lead to a new near-equilibrium state different from the initial distribution for the electrons in the isothermal expansion region. The electron kinetic characteristics in the plume are also very anisotropic. 7.1 Introduction 102 7.1 Introduction In the last two chapters, we considered the expansion of a pre-neutralized plasma plume into a vacuum. However, the emitted plasma is indeed non-neutral near the source due to the huge difference in the electron and ion emitting velocities. Though the ion-electron coupling and beam neutralization processes take place very quickly, their effects on the subsequent plasma plume expansion, especially the electron thermodynamics, is unknown. Previous studies of collisionless plasma plume have mostly focused on the expansion dynamics occuring at the ion time scale. It is thought that the primary role of the electrons is to simply provide a charge neutralization mechanism. Hence, a common approach is to assume that the electrons may be modeled as an equilibrium, massless fluid so to simplify the electron dynamics. Although the electron fluid assumption and thus the hybrid PIC models have been widely used to study the expansion of a collisionless mesothermal plasma plume, the ion-electron coupling and the neutralization processes can only be captured by the full particle PIC model. Recently, several 2D and 3D PIC simulations were carried out to investigate the emission of a finite size mesothermal plasma beam from an ion thruster type source with the beam neutralization process considered [79–84, 146]. In particular, the PIC simulation results from Wang et al. [83] reveal that, during the emission of a plasma beam, a potential well is established between the emission source and the propagating beam front along the beam direction which traps thermal electrons. This chapter extends the work of Wang et al. [83] and presents a more complete full particle PIC simulation study of the emission of a collisionless, mesothermal plasma plume from an ion thruster type source. The focus is to understand the underlying connection between the macroscopic plume structure and microscopic electron kinetic characteristics. 7.2 Simulation Model 103 7.2 Simulation Model The full particle PIC methods described in Sec. 2.2 is used to conduct the simulation study in this chapter. We consider the emission of a collisionless, mesothermal plasma beam into vacuum. At the emission surface, the plasma beam is current-free, J i0 +J e0 = 0, but has a different initial number density, n i0 6=n e0 . J i0 , J e0 , n i0 and n e0 denote the initial ion and electron current density, and ion and electron number density, respectively. Both the electrons and ions are modeled as macro-particles. Fig. 7.1 shows the simulation setup. Note that in this study, the horizontal axis is marked as the z-axis, while the vertical axis is named the x-aixs. In the setup, the simulation domain includes a fraction of the plasma source body. Macro-particles representing cold beam ions and thermal electrons are emitted along the z direction into the simulation domain, which is initially a vacuum. At the emission surface, located at a couple of cells downstream of the source body, the equal electron and ion current density condition for the emitted plasma, J i0 =J e0 , is imposed. The emitted ions are sampled from a drifting Maxwellian distribution with a drifting velocity v 0 and a temperature T i0 for the source. The electrons are thermal and stationary, and follow a stationary Maxwellian distribution. The thermal electrons are injected into the simulation domain through a process similar to the free molecular effusion. Hence, the initial velocity distribution function of the electrons at the emission surface is given by f(v x ) = r m e 2k b T e0 exp m e 2k b T e0 v 2 x ; v x 2 (1; +1) f(v z ) = m e k b T e0 v z exp m e 2k b T e0 v 2 z ; v z 0 (7.1) Note that f(v z ) in Eq. (7.1) is the flux distribution function. Since the electrons and ions have different initial drifting velocities, the condition of J i0 =J e0 leads to an initially non-neutral beam at the emission surface. The full particle ES-PIC simulation described in Sec. 2.2 is carried out. The simulation domain is symmetric with respect to the x = 0 boundary to save computational time. The symmetric boundary condition is applied at the x = 0 surface for the field as well as 7.2 Simulation Model 104 Figure 7.1: Simulation setup for the fully kinetic PIC modeling of collisionless, mesothermal plasma emission the particles. The Neumann boundary condition for the electric field is applied at other boundaries of the simulation domain. Argon and xenon are the most commonly used propellent for EP thrusters. Though we are trying to model the plasma emission from an ion thruster, the main concern in this chapter is indeed the neutralization and the subsequent expansion processes of the plume, phenomena occurring at the ion time scale. We model the proton as the ion species in order to reduce the simulation time needed for ion time scale phenomena but not to deviate too far away from the real applications in engineering. The real proton-to-electron mass ratio m i =m e = 1836:15 is used. It is critical to maintain the v ti0 v 0 v te0 velocity order so the marco-particles injected into the simulation domain represent a mesothermal plasma beam. The temperature ratio of ion to electron at the source exit is taken to be T i0 =T e0 = 0:01. This corresponds to a normalized ion thermal velocity of v ti0 =v te0 = 0:0023. The macroscopic beam velocity for the ions is taken to be v 0 = 0:1v te0 . The Mach number of the beam at the emission exit is M 0 =v 0 =C s0 4:29, where C s0 = p T e0 =m i is the ion acoustic velocity. We note that the values of v 0 =v te0 and v ti0 =v te0 are similar to that of a plasma beam emitted from a typical ion thruster. The center of the plasma emission surface is located at (x;z) = (0; 6). For a typical plasma thruster, the radius of the plasma beam emitted, R b , is much larger than the Debye 7.3 Results and Discussion 105 length at the emission source D0 . In a simulation, the steady state is controlled by the time scale of ion residence time in the simulation domain. In order to obtain a steady state plume expansion, the simulation needs to run for a time duration much longer than the ion plasma period. The beam transient length, L b =v 0 t, also needs to be much larger than the beam radius, L b R b so the results are meaningful for the beam emission problem. These requirements determine the simulation domain size and simulation time. In the simulation, the size of the initial beam size in x is taken to be R b = 20 D0 . The simulation domain used is 600 D0 1000 D0 , with a mesh resolution d = D0 . Note that the domain size in xz corresponds to 30R b 50R b . The source body is represented by a box at (0 ~ x 30; 0 ~ z 4) with a fixed potential ~ body (set to be zero in the simulation). Through the use of the zero Neumann boundary condition for the electric field, the potential at the outer boundary of the simulation domain is floating with respect to ~ body . This simulates plasma emission in space. At steady state, the potential at the outer domain boundary represents the ambient potential relative to the thruster. In the rest of the paper, we use the potential at the ~ z max boundary, ~ (0; 1000), as the reference potential. The simulation was run using a time step resolution of t! pe0 = 0:05 (t! pi0 ' 1:16 10 3 ) for a duration of t! pi0 > 40. At the end of the simulation, the beam transient length is L b 10R b . We note that the plasma plume expansion is still well within the simulation domain at the end of the simulation. Hence, no effects due to momentum and energy transfer at the simulation domain boundary on the simulation result were observed. At each step, 2000 simulation particles (1000 for each species) are injected. The total number of macro-particles is about 80 million at the end of simulation. 7.3 Results and Discussion The simulation results are presented and discussed in this section. The macroscopic flow field properties, such as the electric potential, charge density, electron temperatures, etc., are first 7.3 Results and Discussion 106 ˜ z ˜ x ˜ Φ (a) 0 100 200 300 400 500 600 700 800 0 100 200 300 400 500 600 0 1 2 3 4 5 6 7 8 9 10 ˜ z ˜ x ˜ ρ tot (b) 0 100 200 300 400 500 600 700 800 0 100 200 300 400 500 600 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 Figure 7.2: Electric potential contours (a) and total charge density contours (b) at t! pi0 = 40 (t! pe0 = 1713:6). The values of ~ plotted are with respect to the ambient reference potential. ˜ z ˜ x 0 50 100 150 200 250 300 350 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 ˜ n i,e /˜ n i0 M 0 =4.29 tω pi =40 Figure 7.3: Number density contours (normalized by ~ n i0 ) for ions (colored “ ”) and electrons (colored “ ”) at t! pi0 = 40 (t! pe0 = 1713:6) for m i =m e = 1836 and M 0 = 4:29. The black “ ” represents the Mach line based on M 0 originating from the edge of the emission source, (~ x 0 ; ~ z 0 ) = (20; 6). examined. Then, the microscopic information of electrons, the velocity distribution, in the plume is discussed, and how it relates to the macroscopic plume structure is analyzed. 7.3.1 Macroscopic Plume Structure We first examine the macroscopic plume properties. Fig. 7.2 shows the electric potential contours and total charge density contours at t! pi0 = 40 (t! pe0 = 1713:6). At this time step, the ion beam transient length is L b ' 8:6R b (at z' 177). Fig. 7.3 shows a zoomed-in plot of the ion number density and electron number density contours. The results show that, although the beam is initially non-neutral at the emission surface (z = 6), a well-neutralized beam is formed starting at about 3 D0 downstream. The 7.3 Results and Discussion 107 neutralization process for this emission setup was discussed in detail in Wang et al. [83]. A potential well is established between the beam exit and the beam front in the beam core region which traps the majority of the emitted electrons, and beam neutralization is achieved via the interactions of the trapped electrons and the potential well. Figs. 7.2 and 7.3 also show that the plume exhibits an expansion fan structure. During the mesothermal plasma emission, the expansion outside the core region does not influence the neutralization process inside the plume but determines the potential difference between the beam and the ambient. In a quasi-neutral meosthermal plasma flow expansion, the boundary that separates the undisturbed plasma flow from expansion region would be the first characteristic line of the expansion fan, which is also the Mach line computed based on the Mach number inside the undisturbed region. The relation between the slope of the Mach line and the the Mach number is k = tan(); where' sin 1 1 M 0 (7.2) In Fig. 7.3, the black dashed line is the Mach line intersecting the edge of the source exit (~ x 0 ; ~ z 0 ) = (20; 6) computed using the Mach number at the emission exit, M 0 . We note that, while the plume emission exhibits a similar expansion fan structure in the radial direction, there is a small discrepancy between this computed Mach line and the expansion fan boundary in the simulation result. This small discrepancy is because the plasma plume is non-neutral at the source exit in our simulation setup. It also suggests that the electrons in the plume may not be represented by an isothermal fluid. As the drifting velocity v d v ti , the ion trajectories follow clearly identified streamlines. Fig. 7.4 shows several selected ion streamlines. As v te v d , the electrons follow a mostly thermal motion. Fig. 7.4 also shows the electron temperature contour. The electron temperature T e is computed in the same manner as that presented in Sec. 3.4.2. Note the notation in this chapter is slightly differnt, that is, x denotes the vertical direction and z represents the horizental direction. 7.3 Results and Discussion 108 ˜ z ˜ x (2, 6) (6, 6) (12, 6) (16, 6) ˜ T e,x (a) 0 50 100 150 200 250 300 350 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 B C D ˜ Φ= ˜ Φ inj ˜ z ˜ x (2, 6) (6, 6) (12, 6) (16, 6) ˜ T e,z (b) 0 50 100 150 200 250 300 350 0 20 40 60 80 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 B C D ˜ Φ= ˜ Φ inj Figure 7.4: Contours of electron temperature ~ T e;x (a) and ~ T e;z (b) and four selected ion streamlines at t! pi0 = 40 (t! pe0 = 1713:6). The selected ion streamlines originated at (~ x 0 ; ~ z 0 ) = (2; 6); (6; 6); (12; 6); (16; 6). The boundaries for different plume regions are also shown. The iso-potential line with a value equal to the average at the plume emission plane ( 9:125) is also shown. 0 50 100 150 200 250 300 350 400 0 2 4 6 8 10 12 B C D (a) 0 50 100 150 200 250 300 350 400 -3 -2 -1 0 1 B C D (b) 0 50 100 150 200 250 300 350 400 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 B C D (c) Figure 7.5: 1-D profiles of potential (a), electron number density (b), and electron temperature (c) along the z direction for ~ x = 0, 1 2 R b , and R b at t! pi0 = 40 (t! pe0 = 1713:6). The division between the different regions is based on the ~ x = 0 profile. Fig. 7.4 shows that, starting from the emission source, the ~ T e;j contours exhibit first a very slow cooling region, then a rapid cooling region, and finally a near-isothermal region. The statistical errors in these regions are typically with a few percentage. Further away from the source, some random pattens start to show up in the ~ T e;j contours. Due to the small N p at the far away locations, the statistical error in Fig. 7.4 at these locations can be as high as about 10% to 20%. It is also interesting to note that T e;x and T e;z are different. Hence, the electrons are not only non-uniform but also anisotropic. Fig. 7.5 shows several 1-D profiles (at ~ x = 0, 1 2 R b , and R b ) for the potential, electron density, and electron temperature along the emission direction. These 1-D profiles and the 2-D contour plots shown in Fig. 7.4 suggest that the plume structure consists of several distinct regions (see Fig. 7.4). Immediate adjacent to the beam emission surface is the sheath region. For the simulation setup considered here, the sheath thickness is about 3 4 D0 . Inside the 7.3 Results and Discussion 109 sheath, there is a rapid rising in both the number density and temperature of the electrons. We note that the sheath region is determined by the specific design of plasma emission source, and thus is not the focus of this paper. The focus here is the general macroscopic structure of a collisionless, mesothermal plasma plume downstream of the sheath. We find that the plume structure downstream of the sheath may be divided into three regions: B, C and D. The sheath region is followed by Region B, where the potential, electron density, and electron temperature are nearly constant. Region B is approximately bounded by the expansion boundary shown in Fig. 7.3. The expansion boundary intersects the ~ x = 0 axis at ~ z 50). Outside the Region B is a plasma expansion region. We find that the expansion region may be further divided into a cooling expansion region (Region C) and an isothermal expansion region (Region D) based on the electron temperature characteristics. On the ~ x = 0 axis, Region C is at approximately 50< ~ z < 150 where both the number density and temperature drop quickly. In contrast, in Region D (150< ~ z < 300 on the ~ x = 0 axis),the electron temperature is almost constant while the number density keeps decreasing. Note that Region D is the transition region from the core plume region consisting of both ions and electrons to the outer expansion region consisting of only electrons. The boundaries between Region C and D in Fig. 7.4 can be identified from the electron temperature variation along the ion streamlines. Fig. 7.6 shows the T e -n e relation for the electrons plotted along ion streamlines. (The T e vs. n e points in Region B are clustered together at the end of the curve.) The two distinct expansion regions, the cooling and isothermal expansions, can be clearly observed. In Fig. 7.6, we set log 10 (~ n e ) = 0:25 as the boundary between region C and D, and log 10 (~ n e ) =2 as the boundary between region D and the nearly vacuum region. Cooling is a signature feature in classical gas expansion where the gas is in equilibrium due to inter-molecular collisions. Here, the existence of an electron cooling region shows that the macroscopic properties of the electrons in Region C are somewhat similar to that of an equilibrium fluid even though the plasma is collisionless. This is because a potential well is established inside the beam core region and the electrons are thermalized by the interactions 7.3 Results and Discussion 110 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −1.6 −1.2 −0.8 −0.4 0 0.4 log 10 (˜ n e ) log 10 ( ˜ T e ) (2, 6) (6, 6) (12, 6) (16, 6) C D ˜ T e,x ˜ T e,z Figure 7.6: log 10 (T e )log 10 (n e ) relation along selected ion streamlines at t! pi0 = 40 (t! pe0 = 1713:6). The T e vs. n e points in Region B are clustered together at the end of the curve. with the potential well. The electrons eventually reach a new, near equilibrium state as they exit Region C. This suggests that a fluid model assumption of electrons might be appropriate in this region. Region C is followed by an isothermal region further downstream, Region D. The electrons in Region D are dominated by those electrons at the tail of distribution who can escape the trap created by the potential well in the beam core region. These electrons do not undergo particle-particle collisions because the flow is collisionless. There are also no significant interactions between the electrons and the electric field because the kinetic energy of these electrons is larger than the electric potential field. As a result, these electrons experience little energy dissipation. Hence, they behave almost like the neutral molecules in a free molecular flow and are isothermal. In Fig. 7.4, we also overlay the iso-potential line with a value equal to average beam potential outside the emission source with respect to the ambient potential, = inj . We find that this iso-potential line almost matches the boundary that separates the region C from the region D. This further suggests that the electrons in Region D are those that have overcome the potential well in the beam core region. It is interesting to note that an isothermal expansion controlled by the energetic particles far downstream of the was also seen in the quasi-1D collisionless plasma expansion [78]. 7.3 Results and Discussion 111 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −1.6 −1.2 −0.8 −0.4 0 0.4 log 10 (˜ n e ) log 10 ( ˜ T e ) k =−0.0082 k = 0.5831 k = 0.0157 k = 1.3245 C D (a) (˜ x 0 ,˜ z 0 ) = (2,6) ˜ T e,z ˜ T e,x −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −1.6 −1.2 −0.8 −0.4 0 0.4 log 10 (˜ n e ) log 10 ( ˜ T e ) k =−0.0011 k = 0.6112 k = 0.0294 k = 1.2867 C D (b) (˜ x 0 ,˜ z 0 ) = (4,6) ˜ T e,z ˜ T e,x −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −1.6 −1.2 −0.8 −0.4 0 0.4 log 10 (˜ n e ) log 10 ( ˜ T e ) k =−0.0029 k = 0.6241 k = 0.0715 k = 1.2904 C D (c) (˜ x 0 ,˜ z 0 ) = (6,6) ˜ T e,z ˜ T e,x −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −1.6 −1.2 −0.8 −0.4 0 0.4 log 10 (˜ n e ) log 10 ( ˜ T e ) k =−0.0233 k = 0.6201 k = 0.0891 k = 1.2856 C D (d) (˜ x 0 ,˜ z 0 ) = (8,6) ˜ T e,z ˜ T e,x −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −1.6 −1.2 −0.8 −0.4 0 0.4 log 10 (˜ n e ) log 10 ( ˜ T e ) k =−0.0262 k = 0.5917 k = 0.1059 k = 1.2449 C D (e) (˜ x 0 ,˜ z 0 ) = (10,6) ˜ T e,z ˜ T e,x −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −1.6 −1.2 −0.8 −0.4 0 0.4 log 10 (˜ n e ) log 10 ( ˜ T e ) k =−0.0141 k = 0.5503 k = 0.1033 k = 1.1609 C D (f) (˜ x 0 ,˜ z 0 ) = (12,6) ˜ T e,z ˜ T e,x −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −1.6 −1.2 −0.8 −0.4 0 0.4 log 10 (˜ n e ) log 10 ( ˜ T e ) k = 0.0126 k = 0.5287 k = 0.1620 k = 1.0020 C D (g) (˜ x 0 ,˜ z 0 ) = (14,6) ˜ T e,z ˜ T e,x −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −1.6 −1.2 −0.8 −0.4 0 0.4 log 10 (˜ n e ) log 10 ( ˜ T e ) k = 0.0326 k = 0.5021 k = 0.1700 k = 0.7369 C D (h) (˜ x 0 ,˜ z 0 ) = (16,6) ˜ T e,z ˜ T e,x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 δγ( ˜ Te,x,C)=−4.19( ˜ x0 ˜ Rb ) 3 +3.71( ˜ x0 ˜ Rb ) 2 −1.11( ˜ x0 ˜ Rb )+1.40 δγ( ˜ Te,z,C)= 0.56±0.06 δγ( ˜ Te,x,D)= 0.22( ˜ x0 ˜ Rb )−0.01 δγ( ˜ Te,z,D)= 0.00±0.03 (i) ˜ x 0 / ˜ R b δγ Figure 7.7: log 10 (T e )log 10 (n e ) relation (a) - (h) and variations of polytropic coefficients (i) with different ion streamlines at t! pi0 = 40 (t! pe0 = 1713:6). The T e vs. n e points in Region B are clustered together at the end of the curve. The correlation between the electron temperatures and the electron density along the selected ion streamlines are further investigated quantitatively. Fig. 7.7 plots log 10 ( ~ T e ) log 10 (n e ) relation plus the best fitted the polytropic relations in Region C and D. As T e and n e are plotted in the log scale, the slope (k) of each correlation can be used directly to construct the polytropic coefficient in Eq. (1.14) = 1. The fitted polytropic coefficients in region C and D are shown in Fig. 7.7(i) as a function of the origins of ion streamlines followed. ~ x 0 denotes the distance between the ion streamlines’ origins to the beam center. In both Region C and D, for ~ T e;x is larger than that for ~ T e;z , suggesting that ~ T e;x drops more rapidly as the plasma expands. This is understandable because the expansion outside the beam core is mainly in the x (radial) direction. The polytropic coefficients for ~ T e;z do not 7.3 Results and Discussion 112 significantly change as the streamline origin varies while ( ~ T e;x ) shows strong dependence on the starting point of the streamline at the emission plane, espeically in Region C. ( ~ T e;z ) is given as ( ~ T e;z ) = 8 > > < > > : 0:56 0:06; in Region C 0:00 0:03; in Region D (7.3) The relation between ( ~ T e;x ) and ~ x 0 = ~ R b can be fitted as ( ~ T e;x ) = 8 > > < > > : 4:19 ~ x 0 ~ R b 3 + 3:71 ~ x 0 ~ R b 2 1:11 ~ x 0 ~ R b + 1:40; in Region C 0:22 ~ x 0 ~ R b 0:01; in Region D (7.4) 7.3.2 Microscopic Electron Characteristics The full PIC results are further analyzed to investigate microscopic electron kinetic character- istics. Fig.7.8 shows the instant, local electron velocity distribution functions (EVDFs) along the central axis of the emission source. The EVDFs of both ~ v x and ~ v z are considered. ^ f(~ v x ) and ^ f(~ v z ) are the velocity distributions calculated by direct sampling of the macro-particles contained in one cell at the specified location. Because the far down stream region only contain a small number of macro-particles per cell, we apply the kernel density estimation method with the Gaussian kernel [147] to evaluate the EVDFs. For comparison, we also plot the Maxwellian velocity distribution fit, ^ f M (~ v x ) and ^ f M (~ v z ), using the velocity moments of the electrons at each location in Fig. 7.8. The Maxwellian distribution fits are computed as ^ f M (~ v j ) = 1 q 2 ~ T e;j exp (~ v j h~ v j i) 2 2 ~ T e;j ! (7.5) where the subscript j again denotes x and z, and ~ T e;j andh~ v j i are calculated from Eqs.(3.6). In Fig. 7.8, the local electron kinetic characteristics is shown by the EVDFs near the emission surface (~ z = 6), in the middle of the B region (~ z = 25), near the Region B/C boundary (~ z = 50), in the middle of the C region (~ z = 105), near the Region C/D boundary 7.3 Results and Discussion 113 −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (0∼ 1,6∼ 7) −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (0∼ 1,25∼ 26) −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (0∼ 1,50∼ 51) −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (0∼ 1,105∼ 106) −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (0∼ 1,140∼ 141) −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (0∼ 1,220∼ 221) Figure 7.8: EVDFs at selected locations along the central axis of the emitting source at t! pi = 40. ^ f(~ v x ) (“ ”) and ^ f(~ v z ) (“ ”) are the instant electron velocity distributions sampled from the macro-particles at that time step. The EVDFs are calculated from the macro-particles contained in one cell. The Maxwellian velocity distribution fits using the temperatures and velocities from the macro-particles at each location, ^ f M (~ v x ) (“4”) and ^ f M (~ v z ) (“”), are also plotted for comparison. (~ z = 140), and in the middle of the D region (~ z = 225), respectively. Recall that the electrons are emitted as a Maxwellian distribution in the transverse direction x and as an effusion from a stationary Maxwellian along the beam direction z. We find that, starting from the sheath, both ^ f(~ v z ) and ^ f(~ v x ) first transition from their initial distribution to non-Maxwellian distributions in Region B. They then transition to new, near-Maxwellian distributions at the Region C/Region D boundary, and maintain their respective new, near-equilibrium distribution functions in Region D. These changes in the local EVDFs are consistent with the macroscopic plume properties discussed in the last section. For a collisionless plasma, the kinetic energy associated with the bulk velocity can be dissipated through the collective interactions. Some previous studies of a collisionless plasma beam suggested that a non-equilibrium electron distribution in the beam can be restored to the near-equilibrium state very rapidly through collective interactions such as instabilities[148]. However, for the collisionless plasma plume considered here, no 7.3 Results and Discussion 114 −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (20∼ 21,25∼ 26) −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (20∼ 21,105∼ 106) −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (20∼ 21,220∼ 221) −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (35∼ 36,25∼ 26) −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (35∼ 36,105∼ 106) −4 −3 −2 −1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 ˜ v x,z ˆ f (˜ x,˜ z) = (35∼ 36,220∼ 221) Figure 7.9: EVDFs at selected locations along two off-center axis, ~ x = 20 and ~ x = 35, at t! pi = 40. ^ f(~ v x ) (“ ”) and ^ f(~ v z ) (“ ”) are the instant electron velocity distributions sampled from the macro-particles at that time step. The EVDFs are calculated from the macro-particles contained in one cell. The Maxwellian velocity distribution fits using the temperatures and velocities from the macro-particles at each location, ^ f M (~ v x ) (“4”) and ^ f M (~ v z ) (“”), are also plotted for comparison. instabilities were observed in the simulation. Rather, the primary collective interaction process to the new equilibrium is through the interactions between the electrons trapped inside beam core region and the potential well[83]. Such interactions dissipate the electron kinetic energy associated with the bulk velocity along thez direction and the electron thermal energy along the x direction. The boundary of the potential well is outlined in Fig. 7.4, which approximately matches the boundary between Region C and Region D. In Region D, there is relatively little energy dissipation on electrons because the electric field is weak compared to the electron energy, and the electrons are collisionless. Hence, the electron expansion in Region D is approximately isothermal. The electrons there also maintain approximately the near-equilibrium state established at the Region C/Region D boundary. In the region beyond Region D (~ z > 350), the electron density becomes so low that the region is almost vacuum. Hence, the kinetic characteristics of electrons in this region is 7.4 Summary 115 similar to that of a collisionless neutral gas absorbed by a sink, and the distribution deviates from the near-equilibrium state in Region D. Fig. 7.9 shows the instant, local EVDFs at selected locations along two off-center axes, ~ x = 20 = R b and ~ x = 35. The EVDFs plotted at ~ x = 20 mostly show that they are close to a Maxwellian distribution. From Fig. 7.4, we find that the locations of (~ x; ~ z) = (20; 25) and (~ x; ~ z) = (20; 105) are approximately at the Region C/D boundary and the location of (~ x; ~ z) = (20; 220) is inside Region D. Hence, the EVFDs at these locations are consistent with the discussions presented above. The locations along ~ x = 35 are further outside where the plume density is very low and hence the EVDFs start to deviate from the near-equilibrium. In particular, we note that the EVDF at (~ x; ~ z) = (30; 25) is highly non-Maxwellian. The patten shown in the T e;j contours at that location in Fig. 7.4 could be just a reflection of the second velocity moment of a non-Maxwellian distribution. 7.4 Summary A full particle PIC simulation is carried out to investigate the emission and subsequent expansion of a collisionless, mesothermal plasma plume in this chapter. The plume is initially non-neutral at the emission source. The focus is to understand the connections between the two-dimensional macroscopic plume structure and the microscopic electron kinetics and the effects of beam neutralization process on the expansion. We find that the microscopic electron kinetic characteristics plays an important role in determining the macroscopic plume structure. The plume exhibits several distinctive regions including an undisturbed region immediately downstream of the sheath at the emission source and a plume expansion region. The plume expansion region further includes an inner cooling expansion region and an outer isothermal expansion region where the plume transitions from a quasi-neutral plasma to a rarefied electron-only gas. The inner electron cooling expansion region approximately matches the potential well established by the beam. The interactions between the electrons and the potential well dissipate the electron kinetic energy associated with the bulk velocity 7.4 Summary 116 along the z direction and the electron thermal energy along the x direction and establish a new, near-equilibrium state different from that at the emission source for the electrons. Subsequently, those electrons escaping the potential well in the beam undergo an almost isothermal expansion outside the cooling expansion region until they escape the plume. We find that the electron expansion in the plume is anisotropic and non-equilibrium in most of the plume regions, and the local electron velocity distribution function varies at different plume locations. Hence, the commonly used approach of utilizing a single equilibrium fluid model for electrons in collisionless mesothermal plasma plume expansion in general is not valid. 117 Chapter 8: Summary and Conclusions The collisionless, mesothermal plasma flow dynamics is of fundamental importance for plasma physics and many applications. In the field of space engineering, the plasma-spacecraft interactions and the expansion of a plasma plume from an EP thruster are two very impor- tant subjects closely related to the collisionless, mesothermal plasma flow dynamics. This dissertation presents a comprehensive investigation on these two subjects by using the fully kinetic PIC modeling as well as the hybrid fluid-particle PIC methods. 8.1 Conclusions 8.1.1 Collisionless, Mesothermal Plasma Flow over a Large Con- ducting Plate The interactions between the collisionless, mesothermal plasma and spacecraft are studied by using both hybrid fluid-electron/particle-ion PIC simulations and full particle PIC simulations. The spacecraft is modeled as a large conducting plate, where the term “large” indicates the plate length L p is much larger than the ambient plasma Debye length D0 . The collisionless, mesothermal plasma flow over a conducting plate with a surface potential the same as the ambient plasma’s is first studied. In this study, the most widely used isothermal electron fluid model is first incorporated into the hybrid PIC simulation. The validity of this model is investigated by directly comparing the results between the hybrid and full PIC simulation. We find that the plasma wake may be characterized into two regions based on electron behavior, a fluid expansion region between the initial Mach line and the expansion fan angle at ' 2 0 ' 2=M 0 , and a kinetic electron expansion region beyond ' 2 0 ' 2=M 0 . In the fluid electron expansion region, the approaches based on 8.1 Conclusions 118 the electron fluid approximation and the fully kinetic model lead to a similar result, and thus a fluid approximation for electrons may be applied. The perturbation generated by the plate in this region is analogous to that of the Prandtl-Meyer expansion of a supersonic gas flow. In the kinetic electron expansion region, the fluid approximation for electrons breaks down and a fully kinetic approach must be adopted. The full particle PIC simulations show that the local electron VDFs in this region are strongly non-Maxwellian and the electron “temperature” is strongly anisotropic and exhibits a complex pattern of “cooling” and “heating”. A potential well structure is seen in the wake of the plate. This potential well structure is found to be responsible for the complex kinetic “cooling” and “heating” phenomena associated with the electrons. In addition to the most widely used isothermal electron approximation, the more general polytropic electron fluid model is also incorporated into the hybrid PIC simulations. A set of simulations with different values is carried out and compared to the full particle PIC simulation as well as the isothermal electron fluid based hybrid PIC simulation. The comparison shows that very little improvement is achieved by taking advantage of the polytropic relation for the electron. The failure of the hybrid PIC modelings indicates the fluid approximations for the electron break down, especially in the kinetic expansion region. To quantify the electron kinetic effects on macroscopic plasma flow, a non-equilibrium parameter is proposed. This non-equilibrium parameter takes into account both the electron temperature anisotropy and the deviation of the electron VDF from the Maxwellian VDF. It is shown that the non-equilibrium parameter with a properly chosen threshold accurately predicts where the hybrid PIC simulation fails. Furthermore, the hybrid PIC prediction errors against the full PIC results are found to scale with . The influence of the plate surface potential on the plasma-plate interactions is then studied. The plate surface potential varies in a very wide range such that the plasma-plate interactions cover the quasi-neutral, “thin-sheath”, and “thick-sheath” regimes. Both the fully kinetic PIC and the isothermal electron based hybrid PIC simulations are conducted. It is found that the hybrid PIC modeling is able to describe the ion dynamics in the vicinity of 8.1 Conclusions 119 the plate very well. The difference bewteen the hybrid and full PIC predictions on the ion current collection is very small for all the surface potential values considered in this study. The electron kinetic effects on macroscopic plasma flow in the expansion region are quantified using both the hybrid PIC error and the non-equilibrium parameter . We find that the area where hybrid PIC errors are significant (> 10%) first decreases and then increases as the surface potential ~ p changes from2 to2000. Together with this trend, the enhancement of ~ T e;x or the dropping of ~ T e;y , compared with ambient ~ T e , in the deep expansion region is observed to first decrease and then increase. The failure of the hybrid PIC modeling is primarily due to the breakdown of electron fluid approximation. This is confirmed by the local non-equilibrium degree of electron flow measured by the non-equilibrium parameter . The results also show that the parameter to predict the non-equilibrium degree of a flow has a very good generalizability. 8.1.2 Expansion of a collisionless plasma plume into a vacuum The expansion of a plasma plume from an EP thruster is an important problem involving many fundamental plasam physics issues. It also plays a critical role in many engineering problems, e.g. the contamination or damage to spacecraft facilities caused by plume-spacecraft interactions. Most of the previous modeling studies of EP thruster plume have focused on the physics driven by the ion dynamics. The hybrid PIC modeling dominates this field. However, there is few discussion on the applicability of the electron fluid approximation for modeling the plasma plume expansion, especially validated by the high fidelity fully kinetic approaches. First, the expansion of a 2D collisionless, hypersonic plasma plume into a vacuum is simulated with both the hybrid and fully kinetic PIC models. The plasma is injected from a pre-neutralized source with the paramters similar to an EP plume. The hybrid PIC model considers the electrons as a massless isothermal fluid. The full PIC simulations are carried out using the real ion-to-electron mass ratios for H + , Ar + and Xe + . A direct comparative study of the 2D plasma plume using the two models is carried out. Both the hybrid and full PIC simulations reveal that the plume structure exhibits four distinct expansion regions. 8.1 Conclusions 120 These regions are identified as, from the plasma emission plane to the downtream along the plume direction, the unperturbed, quasi-steady expansion, self-similar expansion and electron front regions, respectively. The formation of the unperturbed region, bounded by the Mach cone, is similar to the process of a supersonic neutral gas expansion. The establishment of the quasi-steady propagation region is due to the combined effect of the forward plume propagation and the backward rarefaction wave propagation. The self-similar expansion mechanism in the hypersonic plasma plume is recognized the same as that in the 1D semi-infinite plasma expansion. The electron temperatures are shown to be highly anisotropic in the plume. T e;x is almost unchanged from the source until the self-similar region, while T e;y drops drastically immediately downstream of the injection plane. A consequence of the anisotropic electron temperatures is the different accelerations that the fastest ions can achieve. A higher electron temperature generally results in a larger ion acceleration. Besides, a smaller ion mass leads to a more rapid electron temperature decrease and thus a lower electron temperature in the transverse direction. The saturated accelerations of the fastest ions are obtained by extrapolating the present data to the time at infinity. The extrapolated saturated ion accelerations scale with p m i =m e . The anisotropic electron thermodynamics leads to very different expansion features along the plume emission direction and the transverse direction. The results show that the expansion along the plume direction is similar to that described by the 1D semi-infinite plasma model and that the electrons are mostly isothermal in that direction. However, the expansion transverse to the plume is analogous to a 1D finite plasma expansion situation, where the thermal energy is limited and thereby the electron cooling effect on the expansion is important. The fully kinetic PIC is used to benchmark the commonly used hybrid PIC simulation of plume simulation. The comparison between the hybrid and full PIC results shows that modeling electrons either as an isothermal fluid or as a polytropic fluid with an isotropic temperatureisover-simplifiedandleadstonotonlyquantitativebutalsoqualitativedifferences. 8.1 Conclusions 121 In order for an electron fluid approximation to be applicable, the model must at least have the capability of resolving the strong anisotropic thermodynamics for electrons. Next, the 3D effects on the collisionless hypersonic plasma plume expansion is examined. The study shows that the 2D plume is qualitatively similar to the 3D plume, but quantitative differences exist. The plasma density drops more rapidly in the 3D plume expansion than in the 2D plume. The more rapid decrease in the plasma density profile in the 3D plume is signified by the hybrid PIC modeling. The electrons exhibit strong anisotropic behaviors in both the 2D and 3D plume expansions. ~ T e;x is almost constant and equal to the source temperature inside the plume core region, but shows abnormal enhancement in the outer expansion region. ~ T e;y decreases as the plasma expands downstream inside the plume core region, but is elevated unusually in the outside region as well. However, there are negligible differences in the spatial profiles of both ~ T e;x and ~ T e;y inside the plume core region between the 2D and 3D cases. A smaller ( ~ T e;y ~ n e ) is observed for the 3D plume due to its faster decrease of the plasma density but similar decrease rate of ~ T e;y . The fully kinetic PIC results show that the fluid models are generally inadequate for the collisionless plasma plume expansion because the electron kinetic behaviors are mainly non-equilibrium, especially in the outer expansion region. Finally, a fully kinetic PIC simulation is carried out to study the emission of a 2D collisionless plasma plume. All the important processes, including the ion-electron coupling, beam neutralization, and plume expansion, relevant to the beam plasma emission from an ion thruster are taken into account in the present model. This study aims at understanding the connections between the macroscopic plume structure and the microscopic electron kinetics in a 2D analysis when the beam neutralization is taken into account. It is found that the macroscopic plume structure exhibits several distinctive regions, including an undisturbed core region, an electron cooling expansion region, and an electron isothermal expansion region. The properties of each region are determined by microscopic electron kinetic characteristics. The division between the undisturbed region and the cooling expansion region approximately matches the Mach line generated at the edge of the emission surface, and that between the 8.2 Recommendations for Further Study 122 cooling expansion region and the isothermal expansion region approximately matches the potential well established in the beam. The interactions between electrons and the potential well lead to a new near-equilibrium state different from the initial distribution for the electrons in the isothermal expansion region. The electron kinetic characteristics in the plume are also very anisotropic. The polytropic coefficient is not a constant but depend on the specific positions as well as the direction of the “equivalent” electron temperatures. 8.2 Recommendations for Further Study This dissertation carries out one of the first direct comparative studies for collisionless plasma flow problems using the hybrid particle-ion/fluid-electron PIC and fully kinetic PIC simulations. One of the main purposes is to better understand the fundamental physical process of the collisionless, mesothermal plasma flows, which is the most important plasma physics problem in space engineering field, by using the high fidelity, fully kinetic PIC method. Moreover, the limitations of commoly used fluid approximations for electrons is also assessed, and the directions of possible future development of more accurate electron fluid models are elucidated by the present high fidelity, fully kinetic results. According to the scope of this dissertation and what has been done, there are a few important works recommended to follow. Effects of Magnetic Field. The current studies ignore the presence of magnetic fields. The magetic field plays an important role in some applications in the field of space engineering. One example is the interactions between the spacecraft and the ionospheric plasma, where the Earth of magnetic field might have an important influence on the electron dynamics, but not the ion dynamics. Another example where the magnetic field have potential effects is the plasma plume from a Hall thruster, the most widely used EP device at present. To include an external magnetic field in the code used for this dissertation has no difficulty and the function is already available as the PIC method is still the electrostatic PIC. 8.2 Recommendations for Further Study 123 Studies of the influence on the electron kinetics, especially the validity of applying the commonly used electron fluid approximations, will be very meaningful for the fundamental plasma physics as well as the engineering applications. Ion temperature. For all the simulations in this dissertation, the ion to electron temperature ratio T i =T e is chosen to be 0.01. Such a treatment is because we hoped to exclude the influence of ion kinetics and to focus only on the electron kinetics. However, in many real applications, such as the LEO and solar wind plasma environment as well as the EP plume plasma, T i =T e varies between0.1 to1, seeing Tables 1.1 and 1.2. Though the dominant parameter in the collisionless, mesothermal plasma flows of concern here is the Mach number M =v d =C s =v d = p k b T e =m i , it is recommended to choose a more realistic ion temperature to better model the space plasma environment or the EP device and to investigate the influence of ion temperatures as well as ion kinetics in the future. Multiple Ion Species in Wake Study. The space environmental plasmas seen by a spacecraft may contain multiple ion species. For example, as shown in Table 1.1, the plasma at the altitude of 100 400 km consists of a majorO + species and a minonH + . Hence, it will be important to consider the space plasma environment in a more realistic way. More importantly, some previous studies based on simple 1D models have shown that there present interesting new phenomena in the plasmas composed of multiple ion species with different masses and thermal energies, depending on their relative concentrations [149–152]. It is necessary to investigate how such new phenomena appear in the more realistic 2D or 3D configurations. Upstream condition in Plume Study. In this dissertation, the upstream condition of the plume is considered to be uniform with zero angle of divergence. This is a highly ideal situation. For plumes, the plasma at the exit surface of the thruster always has a non-uniform density distribution as well as an angle of divergence due to the design of the EP device. The studies of the EP plume 8.2 Recommendations for Further Study 124 that rigorously take the design relevant issues into consideration are recommended for the future. Electron Mixing in Plume Study. One very important problem associated with the plasma plume from an EP thruster is the mixing of electrons with exhaust ion beam. This is the pre-process of the ion- electron coupling and the subsequent propagation of a neutralized plasma beam. Such a problem is not only interesting to the fundamental plasma physics but also useful to real applications. For example, the initial temperature of the electrons emitted from a cathode neutralizer is at most a few tenths eV, but this temperature can reach as high as several eVs. There are a lot of debates on what is the heating mechanism for the electrons. Hence, this is a crucial problem deserving great efforts. However, a proper fully kinetic PIC simulation with an aim of resolving all the physics is extremely computationally demanding. The emission surface at the cathode neutralizer is small, and the electron density there is very high. As a result, one needs to employ a very fine mesh near the cathode in order to properly resolve the electron dynamics. The fine mesh not only increase the total number of cells also limits the time step to be used in a PIC simulation. Accordingly, better strategies of handling the mesh, time step and parallelization of the code are required in order to carry out such a study accurately and efficiently. Further evaulations on current practice in modeling spacecraft charging and EP plume. The ultimate goal of the present study is to provide valuable guidance on the engineering design with the help of high fidelity, fully kinetic view. As a first step, the importance of the electron kinetics in the spacecraft charging and EP plume are studied in great detail by using the fully kinetic method in this dissertation. In particular, the errors induced by the Boltzmann and polytropic relation based electron fluid models, which are widely adopted in engineering oriented researches, are quantitatively evaluated with 8.2 Recommendations for Further Study 125 the simple geometric models of spacecraft and EP thruster. To better help engineers and their desgins, fully kinetic insights into the influence of electron kinetics on spacecraft charging and EP plume problems with more complex but realistic geometries deserve futher studies. BIBLIOGRAPHY 126 Bibliography [1] J. Wang and D. Hastings, “Ionospheric plasma flow over large high-voltage space platforms. II: The formation and structure of plasma wake,” Physics of Fluids B– Plasma Physics, vol. 4, pp. 1615–1629, Jun 1992. [2] P. Mora, “Plasma expansion into a vacuum,” Physical Review Letters, vol. 90, no. 18, p. 185002, 2003. [3] Y. L. Al’pert, A. V. Gurevich, and L. P. Pitaevskii, Space physics with artificial satellites. Consultants Bureau, 1965. [4] Y. L. Al’pert, The Near-Earth and Interplanetary Plasma. Cambridge University Press, 1983. [5] D. Han, Particle-in-Cell Simulations of Plasma Interactions with Asteroidal and Lunar Surfaces. PhD thesis, University of Southern California, 2015. [6] W. A. Hoskins, R. J. Cassady, O. Morgan, R. M. Myers, F. Wilson, D. Q. King, and K. DeGrys, “30 years of electric propulsion flight experience at aerojet rocketdyne,” in 33rd International Electric Propulsion Conference, (Washington, D.C., USA), IEPC- 2013-439, October 2013. [7] I. D. Boyd, “Review of Hall thruster plume modeling,” Journal of Spacecraft and Rockets, vol. 38, no. 3, pp. 381–387, 2001. [8] J. L. Polansky, Laboratory Investigations of the Near Surface Plasma Field and Charging at the Lunar Terminator. PhD thesis, University of Southern California, 2013. [9] J. E. Polk, R. Kakuda, J. Anderson, J. R. Brophy, V. Rawlin, M. Patterson, J. Sovey, and J. Hamley, “Validation of the NSTAR ion propulsion system on the Deep Space One mission: Overview and initial results,” in 35th AIAA/ASME/SAE/ASEE Joint Conference & Exhibit, (Los Angeles, CA, USA), AIAA 1999-2274, June 1999. [10] J. Wang, D. E. Brinza, D. T. Young, J. E. Nordholt, J. E. Polk, M. D. Henry, R. Goldstein, J. J. Hanley, D. J. Lawrence, and M. Shappirio, “Deep Space One investigations of ion propulsion plasma environment,” Journal of Spacecraft and Rockets, vol. 37, no. 5, pp. 545–555, 2000. [11] J. S. Sovey, V. K. Rawlin, and M. J. Patterson, “Ion propulsion development projects in us: Space electric rocket test i to Deep Space 1,” Journal of Propulsion and Power, vol. 17, no. 3, pp. 517–526, 2001. BIBLIOGRAPHY 127 [12] F. F. Chen, Introduction to Plasma Physics and Controlled Fusion, vol. 1. Springer, 2006. [13] S. Chapman and T. G. Cowling, The Mathematical Theory of Non-uniform Gases. 3 ed., 1991. [14] I. Langmuir, “Scattering of electrons in ionized gases,” Physical Review, vol. 26, no. 5, pp. 585–613, 1925. [15] D. Gabor, E. A. Ash, and D. Dracott, “Langmuir’s paradox,” Nature, vol. 176, no. 4489, pp. 916–919, 1955. [16] A. V. Gurevich, L. V. Pariiskaya, and L. P. Pitaevskii, “Self-similar motion of rarefied plasma,” Soviet Physics JETP, vol. 22, no. 2, pp. 449–454, 1966. [17] J. E. Allen and J. G. Andrews, “A note on ion rarefaction waves,” Journal of Plasma Physics, vol. 4, no. 1, pp. 187–194, 1970. [18] Y. Huang, Y. Bi, X. Duan, X. Lan, N. Wang, X. Tang, and Y. He, “Self-similar neutral- plasma isothermal expansion into a vacuum,” Applied Physics Letters, vol. 92, no. 3, p. 031501, 2008. [19] J. S. Pearlman and R. L. Morse, “Maximum expansion velocities of laser-produced plasmas,” Physical Review Letters, vol. 40, no. 25, pp. 1652–1655, 1978. [20] J. E. Crow, P. L. Auer, and J. E. Allen, “The expansion of a plasma into a vacuum,” Journal of Plasma Physics, vol. 14, no. 1, pp. 65–76, 1975. [21] J. Denavit, “Collisionless plasma expansion into a vacuum,” Physics of Fluids, vol. 22, no. 7, pp. 1384–1392, 1979. [22] J. Denavit, “Numerical simulation of plasmas with periodic smoothing in phase space,” Journal of Computational Physics, vol. 9, no. 1, pp. 75–98, 1972. [23] P. Mora and R. Pellat, “Self-similar expansion of a plasma into a vacuum,” Physics of Fluids, vol. 22, no. 12, pp. 2300–2304, 1979. [24] C. Sack and H. Schamel, “Evolution of a plasma expanding into vacuum,” Plasma Physics and Controlled Fusion, vol. 27, no. 7, p. 717, 1985. [25] C. Sack and H. Schamel, “Plasma expansion into vacuum – a hydrodynamic approach,” Physics Reports, vol. 156, no. 6, pp. 311–395, 1987. [26] G. Manfredi, S. Mola, and M. Feix, “Rescaling methods and plasma expansions into vacuum,” Physics of Fluids B: Plasma Physics, vol. 5, no. 2, pp. 388–401, 1993. [27] A. V. Baitin and K. M. Kuzanyan, “A self-similar solution for expansion into a vacuum of a collisionless plasma bunch,” Journal of Plasma Physics, vol. 59, no. 1, pp. 83–90, 1998. [28] D.S.DorozhkinaandV.E.Semenov, “Exactsolutionofvlasovequationsforquasineutral expansion of plasma bunch into vacuum,” Physical Review Letters, vol. 81, no. 13, p. 2691, 1998. BIBLIOGRAPHY 128 [29] V. F. Kovalev and V. Y. Bychenkov, “Analytic solutions to the vlasov equations for expanding plasmas,” Physical Review Letters, vol. 90, no. 18, p. 185004, 2003. [30] P. Mora, “Thin-foil expansion into a vacuum,” Physical Review E, vol. 72, no. 5, p. 056401, 2005. [31] P. Mora, “Collisionless expansion of a Gaussian plasma into a vacuum,” Physics of Plasmas, vol. 12, no. 11, p. 112102, 2005. [32] T. Grismayer and P. Mora, “Influence of a finite initial ion density gradient on plasma expansion into a vacuum,” Physics of Plasmas, vol. 13, no. 3, p. 032103, 2006. [33] P. Mora and T. Grismayer, “Rarefaction acceleration and kinetic effects in thin-foil expansion into a vacuum,” Physical Review Letters, vol. 102, no. 14, p. 145001, 2009. [34] T. Grismayer, P. Mora, J. Adam, and A. Héron, “Electron kinetic effects in plasma expansion and ion acceleration,” Physical Review E, vol. 77, no. 6, p. 066407, 2008. [35] M. Murakami and M. M. Basko, “Self-similar expansion of finite-size non-quasi-neutral plasmas into vacuum: Relation to the problem of ion acceleration,” Physics of Plasmas, vol. 13, no. 1, p. 012105, 2006. [36] M. A. Morgan, C. Chan, D. L. Cooke, and M. F. Tautz, “The dynamics of charged particles in the near wake of a very negatively charged body–Laboratory experiment and numerical simulation,” IEEE Transactions on Plasma Science, vol. 17, no. 2, pp. 220–227, 1989. [37] É. Coggiola and A. Soubeyran, “Mesothermal plasma flow around a negatively wake side biased cylinder,” Journal of Geophysical Research: Space Physics, vol. 96, no. A5, pp. 7613–7621, 1991. [38] J. Wang and D. Hastings, “Ionospheric plasma flow over large high-voltage space platforms. I: Ion-plasma-time scale interactions of a plate at a zero-angle of attack,” Physics of Fluids B–Plasma Physics, vol. 4, pp. 1597–1614, Jun 1992. [39] P. Guio and H. Pécseli, “Phase space structures generated by an absorbing obstacle in a streaming plasma,” Geophysical Research Letters, vol. 31, no. 3, 2004. [40] G. Murphy, J. Pickett, N. D’angelo, and W. S. Kurth, “Measurements of plasma parameters in the vicinity of the space shuttle,” Planetary and Space Science, vol. 34, no. 10, pp. 993–1004, 1986. [41] G. B. Murphy, D. L. Reasoner, A. Tribble, N. D’Angelo, J. S. Pickett, and W. S. Kurth, “The plasma wake of the shuttle orbiter,” Journal of Geophysical Research: Space Physics, vol. 94, no. A6, pp. 6866–6872, 1989. [42] W. A. Oran, U. Samir, N. H. Stone, and E. G. Fontheeim, “Laboratory observations of electron temperature in the wake of a sphere in a streaming plasma,” Planetary and Space Science, vol. 23, no. 7, pp. 1081–1083, 1975. [43] V. A. Shuvalov, “Structure of the near wake behind a sphere in the flow of a nonequi- librium rarefied plasma,” Geomagnetism and Aeronomy, vol. 19, pp. 440–443, 1979. BIBLIOGRAPHY 129 [44] V. A. Shuvalov, “Structure of the near wake behind a cylinder in a nonequilibrium rarefied plasma flow,” Geomagnetism and Aeronomy, vol. 20, pp. 425–429, 1980. [45] M. Morgan, C. Chan, and R. C. Allen, “A laboratory study of the electron temperature in the near wake of a conducting body,” Geophysical Research Letters, vol. 14, no. 11, pp. 1170–1173, 1987. [46] N. Singh, U. Samir, K. H. Wright, and N. H. Stone, “A possible explanation of the electron temperature enhancement in the wake of a satellite,” Journal of Geophysical Research: Space Physics, vol. 92, no. A6, pp. 6100–6106, 1987. [47] N. Singh, K. H. Wright, N. H. Stone, U. Samir, and K. S. Hwang, “On the interpretation of measured ion streams in the wake of the shuttle orbiter in terms of plasma expansion processes,” Journal of Geophysical Research: Space Physics, vol. 94, no. A9, pp. 12075– 12080, 1989. [48] T. Umeda, T. Kimura, K. Togano, K. Fukazawa, Y. Matsumoto, T. Miyoshi, N. Terada, T. K. Nakamura, and T. Ogino, “Vlasov simulation of the interaction between the solar wind and a dielectric body,” Physics of Plasmas, vol. 18, no. 1, p. 012908, 2011. [49] T. Nakagawa, “Ion entry into the wake behind a nonmagnetized obstacle in the solar wind: Two-dimensional particle-in-cell simulations,” Journal of Geophysical Research: Space Physics, vol. 118, no. 5, pp. 1849–1860, 2013. [50] M. I. Zimmerman, W. M. Farrell, and A. R. Poppe, “Grid-free 2D plasma simulations of the complex interaction between the solar wind and small, near-Earth asteroids,” Icarus, vol. 238, pp. 77–85, 2014. [51] M. Martinez-Sanchez and J. E. Pollard, “Spacecraft electric propulsion - an overview,” Journal of Propulsion and Power, vol. 14, no. 5, pp. 688–699, 1998. [52] E. Y. Choueiri, “A critical history of electric propulsion: The first 50 years (1906-1956),” Journal of Propulsion and Power, vol. 20, no. 2, pp. 193–203, 2004. [53] R. G. Jahn, Physics of Electric Propulsion. Dover Publications, 2006. [54] D. M. Goebel and I. Katz, Fundamentals of Electric Propulsion: Ion and Hall Thrusters. John Wiley & Sons, 2008. [55] E. Y. Choueiri, “New dawn for electric rockets,” Scientific American, vol. 300, no. 2, pp. 58–65, 2009. [56] D. Winske, L. Yin, N. Omidi, H. Karimabadi, and K. Quest, “Hybrid simulation codes: Past, present and future – a tutorial,” in Space Plasma Simulation, pp. 136–165, Springer, 2003. [57] J. Wang and J. Brophy, “3-D Monte-Carlo Particle-in-Cell simulations of ion thruster plasma interactions,” in 32nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, (San Deigo, CA, USA), AIAA 1995-2826, July 1996. BIBLIOGRAPHY 130 [58] J. Wang, J. Brophy, and D. Brinza, “3-D simulations of NSTAR ion thruster plasma envi- ronment,” in 32nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, (Lake Buena Vista, FL, USA), AIAA 1996-3202, July 1996. [59] J. Wang, D. Brinza, and M. Young, “Three-dimensional particle simulations of ion propulsion plasma environment for Deep Space 1,” Journal of Spacecraft and Rockets, vol. 38, no. 3, pp. 433–440, 2001. [60] R.KafafyandJ.Wang, “AhybridgridimmersedfiniteelementParticle-in-Cellalgorithm for modeling spacecraft–plasma interactions,” IEEE Transactions on Plasma Science, vol. 34, no. 5, pp. 2114–2124, 2006. [61] J. Wang, Y. Cao, R. Kafafy, J. Pierru, and V. K. Decyk, “Simulations of ion thruster plume–spacecraft interactions on parallel supercomputer,” IEEE Transactions on Plasma Science, vol. 34, no. 5, pp. 2148–2158, 2006. [62] R. I. S. Roy and D. E. Hastings, “Three-dimensional modeling of dual ion-thruster plumes for spacecraft contamination,” Journal of Spacecraft and Rockets, vol. 33, no. 4, pp. 519–524, 1996. [63] R. I. S. Roy, D. E. Hastings, and S. Taylor, “Three-dimensional plasma Particle-in-Cell calculations of ion thruster backflow contamination,” Journal of Computational Physics, vol. 128, no. 1, pp. 6–18, 1996. [64] M. J. Mandell, V. A. Davis, D. L. Cooke, A. T. Wheelock, and C. J. Roth, “Nascap-2k spacecraft charging code overview,” IEEE Transactions on Plasma Science, vol. 34, no. 5, pp. 2084–2093, 2006. [65] T. Muranaka, S. Hosoda, J.-H. Kim, S. Hatta, K. Ikeda, T. Hamanaga, M. Cho, H. Usui, H. O. Ueda, K. Koga, et al., “Development of multi-utility spacecraft charging analysis tool (muscat),” IEEE Transactions on Plasma Science, vol. 36, no. 5, pp. 2336–2349, 2008. [66] J.-F. Roussel, F. Rogier, G. Dufour, J.-C. Mateo-Velez, J. Forest, A. Hilgers, D. Rodgers, L. Girard, and D. Payan, “Spis open-source code: Methods, capabilities, achievements, and prospects,” IEEE Transactions on Plasma Science, vol. 36, no. 5, pp. 2360–2368, 2008. [67] B. Korkut and D. A. Levin, “Three-dimensional simulations of backflows from ion thruster plumes using unstructured grid refinement,” Journal of Propulsion and Power, vol. 33, no. 1, pp. 264–275, 2017. [68] B. Korkut and T. O. Levin, Deborah A, “Simulations of ion thruster plumes in ground facilities using adaptive mesh refinement,” Journal of Propulsion and Power, vol. 33, no. 3, pp. 681–696, 2017. [69] R. I. S. Roy, D. E. Hastings, and N. A. Gatsonis, “Ion-thruster plume modeling for backflow contamination,” Journal of Spacecraft and Rockets, vol. 33, no. 4, pp. 525–534, 1996. BIBLIOGRAPHY 131 [70] R. I. S. Roy, D. E. Hastings, and N. A. Gatsonis, “Numerical study of spacecraft contamination and interactions by ion-thruster effluents,” Journal of Spacecraft and Rockets, vol. 33, no. 4, pp. 535–542, 1996. [71] D. B. VanGilder, I. D. Boyd, and M. Keidar, “Particle simulations of a Hall thruster plume,” Journal of Spacecraft and Rockets, vol. 37, no. 1, pp. 129–136, 2000. [72] S.-W. Kim, J. E. Foster, and A. D. Gallimore, “Very-near-field plume study of a 1.35 kw spt-100,” in 32nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, (Lake Buena Vista, FL, USA), AIAA 1996-2972, July 1996. [73] I. D. Boyd and J. T. Yim, “Modeling of the near field plume of a Hall thruster,” Journal of applied physics, vol. 95, no. 9, pp. 4575–4584, 2004. [74] M. Merino and E. Ahedo, “Influence of electron and ion thermodynamics on the mag- netic nozzle plasma expansion,” in 33rd International Electric Propulsion Conference, (Washington, D.C., USA), IEPC-2013-247, Oct 2013. [75] M. Merino and E. Ahedo, “Influence of electron and ion thermodynamics on the magnetic nozzle plasma expansion,” IEEE Transactions on Plasma Science, vol. 43, no. 1, pp. 244–251, 2015. [76] M. Merino, F. Cichocki, and E. Ahedo, “A collisionless plasma thruster plume expansion model,” Plasma Sources Science and Technology, vol. 24, no. 3, p. 035006, 2015. [77] F. Cichocki, M. Merino, E. Ahedo, Y. Hu, and J. Wang, “Fluid vs PIC modeling of a plasma plume expansion,” in 34th International Electric Propulsion Conference, IEPC-2015-IEPC, (Kobe, Hyogo, Japan), IEPC-2015-420, Oct 2015. [78] M. Martinez-Sanchez, J. Navarro-Cavallé, and E. Ahedo, “Electron cooling and finite potential drop in a magnetized plasma expansion,” Physics of Plasmas, vol. 22, no. 5, p. 053501, 2015. [79] C. Othmer, K. Glassmeier, U. Motschmann, J. Schüle, and C. Frick, “Three-dimensional simulations of ion thruster beam neutralization,” Physics of Plasmas, vol. 7, no. 12, pp. 5242–5251, 2000. [80] C. Othmer, K. Glassmeier, U. Motschmann, and I. Richter, “Numerical parameter studies of ion-thruster-beam neutralization,” Journal of propulsion and power, vol. 19, no. 5, pp. 953–963, 2003. [81] A. Wheelock, D. Cooke, and N. A. Gatsonis, “Investigation of ion beam neutralization processes with 2D and 3D PIC simulations,” Computer Physics Communications, vol. 164, pp. 336–343, 2004. [82] L. Brieda and J. Wang, “Modeling of ion thruster beam neutralization using a fully kinetic ES-PIC code,” in 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, (Tucson, AZ, USA), AIAA 2005-4045, July 2005. [83] J. Wang, O. Chang, and Y. Cao, “Electron-ion coupling in mesothermal plasma beam emission: Full particle PIC simulations,” IEEE Transactions on Plasma Science, vol. 40, no. 2, pp. 230–236, 2012. BIBLIOGRAPHY 132 [84] H. Usui, A. Hashimoto, and Y. Miyake, “Electron behavior in ion beam neutralization in electric propulsion: Full Particle-in-Cell simulation,” in Journal of Physics: Conference Series, vol. 454, p. 012017, 2013. [85] C. Z. Cheng and G. Knorr, “The integration of the Vlasov equation in configuration space,” Journal of Computational Physics, vol. 22, no. 3, pp. 330–351, 1976. [86] E. Sonnendrücker, J. Roche, P. Bertrand, and A. Ghizzo, “The semi-Lagrangian method for the numerical resolution of the Vlasov equation,” Journal of computational physics, vol. 149, no. 2, pp. 201–220, 1999. [87] T. Nakamura and T. Yabe, “Cubic interpolated propagation scheme for solving the hyper-dimensional VlasovâĂŤPoisson equation in phase space,” Computer Physics Communications, vol. 120, no. 2-3, pp. 122–154, 1999. [88] J. A. Carrillo and F. Vecil, “Nonoscillatory interpolation methods applied to Vlasov- based models,” SIAM Journal on Scientific Computing, vol. 29, no. 3, pp. 1179–1206, 2007. [89] J. Qiu and A. Christlieb, “A conservative high order semi-Lagrangian WENO method for the Vlasov equation,” Journal of Computational Physics, vol. 229, no. 4, pp. 1130–1149, 2010. [90] J. Qiu and C. Shu, “Conservative semi-Lagrangian finite difference WENO formulations with applications to the Vlasov equation,” Communications in Computational Physics, vol. 10, no. 4, pp. 979–1000, 2011. [91] X. Cai, J. Qiu, and J. Qiu, “A conservative semi-Lagrangian HWENO method for the Vlasov equation,” Journal of Computational Physics, vol. 323, pp. 95–114, 2016. [92] F. Filbet, E. Sonnendrücker, and P. Bertrand, “Conservative numerical schemes for the Vlasov equation,” Journal of Computational Physics, vol. 172, no. 1, pp. 166–187, 2001. [93] F. Filbet and E. Sonnendrücker, “Comparison of Eulerian Vlasov solvers,” Computer Physics Communications, vol. 150, no. 3, pp. 247–266, 2003. [94] N. Crouseilles, M. Mehrenberger, and E. Sonnendrücker, “Conservative semi-Lagrangian schemes for Vlasov equations,” Journal of Computational Physics, vol. 229, no. 6, pp. 1927–1953, 2010. [95] T. Umeda, Y. Nariyuki, and D. Kariya, “A non-oscillatory and conservative semi- Lagrangian scheme with fourth-degree polynomial interpolation for solving the Vlasov equation,” Computer Physics Communications, vol. 183, no. 5, pp. 1094–1100, 2012. [96] J.-M.Qiu andC.-W.Shu, “Positivitypreserving semi-LagrangiandiscontinuousGalerkin formulation: theoreticalanalysisandapplicationtotheVlasov–Poissonsystem,” Journal of Computational Physics, vol. 230, no. 23, pp. 8386–8409, 2011. [97] W. Guo and J. Qiu, “Hybrid semi-Lagrangian finite element-finite difference methods for the Vlasov equation,” Journal of Computational Physics, vol. 234, pp. 108–132, 2013. BIBLIOGRAPHY 133 [98] S. Vadlamani, S. E. Parker, Y. Chen, and C. Kim, “The particle-continuum method: an algorithmic unification of particle-in-cell and continuum methods,” Computer Physics Communications, vol. 164, no. 1-3, pp. 209–213, 2004. [99] N. Crouseilles, T. Respaud, and E. Sonnendrücker, “A forward semi-Lagrangian method for the numerical solution of the Vlasov equation,” Computer Physics Communications, vol. 180, no. 10, pp. 1730–1745, 2009. [100] J. A. Rossmanith and D. C. Seal, “A positivity-preserving high-order semi-Llagrangian discontinuous Galerkin scheme for the Vlasov–Poisson equations,” Journal of Computa- tional Physics, vol. 230, no. 16, pp. 6203–6232, 2011. [101] A. Christlieb, W. Guo, M. Morton, and J. Qiu, “A high order time splitting method based on integral deferred correction for semi-Lagrangian Vlasov simulations,” Journal of Computational Physics, vol. 267, pp. 7–27, 2014. [102] N. V. Elkina and J. Büchner, “A new conservative unsplit method for the solution of the Vlasov equation,” Journal of Computational Physics, vol. 213, no. 2, pp. 862–875, 2006. [103] J. W. Banks and J. A. F. Hittinger, “A new class of nonlinear finite-volume methods for Vlasov simulation,” IEEE Transactions on Plasma Science, vol. 38, no. 9, pp. 2198–2207, 2010. [104] J. A. Hittinger and J. W. Banks, “Block-structured adaptive mesh refinement algorithms for Vlasov simulation,” Journal of Computational Physics, vol. 241, pp. 118–140, 2013. [105] R. E. Heath, I. M. Gamba, P. J. Morrison, and C. Michler, “A discontinuous Galerkin method for the Vlasov–Poisson system,” Journal of Computational Physics, vol. 231, no. 4, pp. 1140–1174, 2012. [106] X. Cai, W. Guo, and J. Qiu, “A high order semi-Lagrangian discontinuous Galerkin method for Vlasov–Poisson simulations without operator splitting,” Journal of Compu- tational Physics, vol. 354, pp. 529–551, 2018. [107] R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles. CRC Press, 1988. [108] C. K. Birdsall and A. B. Langdon, Plasma physics via computer simulation. CRC press, 2004. [109] J. M. Dawson, “Particle simulation of plasmas,” Reviews of modern physics, vol. 55, no. 2, p. 403, 1983. [110] C. K. Birdsall, “Particle-in-cell charged-particle simulations, plus Monte Carlo collisions with neutral atoms, PIC-MCC,” IEEE Transactions on Plasma Science, vol. 19, no. 2, pp. 65–85, 1991. [111] J. P. Verboncoeur, “Particle simulation of plasmas: review and advances,” Plasma Physics and Controlled Fusion, vol. 47, no. 5A, p. A231, 2005. BIBLIOGRAPHY 134 [112] D. Y. Oh, D. E. Hastings, C. M. Marrese, J. M. Haas, and A. D. Gallimore, “Modeling of stationary plasma thruster-100 thruster plumes and implications for satellite design,” J. Propul. Power, vol. 15, no. 2, pp. 345–357, 1999. [113] F. Taccogna, D. Pagano, F. Scortecci, and A. Garulli, “Three-dimensional plume simulation of multi-channel thruster configuration,” Plasma Sources Sci. Technol., vol. 23, no. 6, p. 065034, 2014. [114] F. Cichocki, A. Domínguez-Vázquez, M. Merino, and E. Ahedo, “Hybrid 3d model for the interaction of plasma thruster plumes with nearby objects,” Plasma Sources Sci. Technol., vol. 26, no. 12, p. 125008, 2017. [115] C. Cai, “Numerical studies on plasma plume flows from a cluster of electric propulsion devices,” Aerosp. Sci. Technol., vol. 41, pp. 134–143, 2015. [116] A. L. Ortega, I. Katz, I. G. Mikellides, and D. M. Goebel, “Self-consistent model of a high-power hall thruster plume,” IEEE Trans. Plasma Sci., vol. 43, no. 9, pp. 2875–2886, 2015. [117] M. Ikegawa and J. Kobayashi, “Development of a rarefield gas flow simulator using the direct-simulation Monte Carlo method,” JSME International Journal Ser. 2, vol. 33, no. 3, pp. 463–467, 1990. [118] R. P. Nance, D. B. Hash, and H. Hassan, “Role of boundary conditions in monte carlo simulation of microelectromechanical systems,” Journal of Thermophysics and Heat Transfer, vol. 12, no. 3, pp. 447–449, 1998. [119] C. Cai, I. D. Boyd, J. Fan, and G. V. Candler, “Direct simulation methods for low- speed microchannel flows,” Journal of Thermophysics and Heat Transfer, vol. 14, no. 3, pp. 368–378, 2000. [120] J.-S. Wu and K.-C. Tseng, “Analysis of micro-scale gas flows with pressure boundaries using direct simulation Monte Carlo method,” Computers & Fluids, vol. 30, no. 6, pp. 711–735, 2001. [121] Y. Fang and W. W. Liou, “Computations of the flow and heat transfer in microdevices using DSMC with implicit boundary conditions,” Journal of Heat Transfer, vol. 124, no. 2, pp. 338–345, 2002. [122] M. Wang and Z. Li, “Simulations for gas flows in microgeometries using the direct simulation Monte Carlo method,” International Journal of Heat and Fluid Flow, vol. 25, no. 6, pp. 975–985, 2004. [123] E. Farbar and I. D. Boyd, “Subsonic flow boundary conditions for the direct simulation Monte Carlo method,” Computers & Fluids, vol. 102, pp. 99–110, 2014. [124] T. A. Davis, Direct Methods for Sparse Linear Systems. SIAM, 2006. [125] T. Tanaka, “Configurations of the solar wind flow and magnetic field around the planets with no magnetic field: Calculation by a new MHD simulation scheme,” J. Geophys. Res.: Space Phys., vol. 98, no. A10, pp. 17251–17262, 1993. BIBLIOGRAPHY 135 [126] S. Ledvina, Y.-J. Ma, and E. Kallio, “Modeling and simulating flowing plasmas and related phenomena,” in Comparative Aeronomy, pp. 143–189, 2008. [127] H. Shimazu, “Three-dimensional hybrid simulation of solar wind interaction with unmagnetized planets,” J. Geophys. Res.: Space Phys., vol. 106, no. A5, pp. 8333–8342, 2001. [128] S. Simon, T. Bagdonat, U. Motschmann, and K.-H. Glassmeier, “Plasma environment of magnetized asteroids: a 3-D hybrid simulation study,” in Annales Geophysicae, vol. 24, pp. 407–414, 2006. [129] J.Wang, P.Leung, H.Garrett, andG.Murphy, “Multibody-plasmainteractions-charging in the wake,” J. Spacecr. Rockets, vol. 31, no. 5, pp. 889–894, 1994. [130] R. E. Ergun, D. M. Malaspina, S. D. Bale, J. P. McFadden, D. E. Larson, F. S. Mozer, N. Meyer-Vernet, M. Maksimovic, P. J. Kellogg, and J. R. Wygant, “Spacecraft charging and ion wake formation in the near-Sun environment,” Phys. Plasmas, vol. 17, no. 7, p. 072903, 2010. [131] U. Samir, K. H. Wright, and N. H. Stone, “The expansion of a plasma into a vac- uum: Basic phenomena and processes and applications to space plasma physics,” Rev. Geophys., vol. 21, no. 7, pp. 1631–1646, 1983. [132] E. M. Lifshitz and L. P. Pitaevskii, Course of Theoretical Physics: Physical Kinetics, vol. 10. 1981. [133] U. Samir and G. L. Wrenn, “Experimental evidence of an electron temperature enhance- ment in the wake of an ionospheric satellite,” Planet. Space Sci., vol. 20, no. 6, pp. 899– 904, 1972. [134] J. M. Illiano and L. R. O. Storey, “Apparent enhancement of electron temperature in the wake of spherical probe in a flowing plasma,” Planet. Space Sci., vol. 22, no. 5, pp. 873–878, 1974. [135] B. E. Troy, E. J. Maier, and U. Samir, “Electron temperatures in the wake of an ionospheric satellite,” J. Geophys. Res., vol. 80, no. 7, pp. 993–997, 1975. [136] U. Samir, N. H. Stone, and K. H. Wright, “On plasma disturbances caused by the motion of the space shuttle and small satellites: A comparison of in situ observations,” J. Geophys. Res.: Space Phys., vol. 91, no. A1, pp. 277–285, 1986. [137] G. A. Bird, “Breakdown of translational and rotational equilibrium in gaseous expan- sions,” AIAA J., vol. 8, no. 11, pp. 1998–2003, 1970. [138] W.-L. Wang and I. D. Boyd, “Predicting continuum breakdown in hypersonic viscous flows,” Phys. Fluids, vol. 15, no. 1, pp. 91–100, 2003. [139] P.-H. Chen, I. Boyd, and J. Camberos, “Assessment of entropy generation rate as a predictor of continuum breakdown,” in 36th AIAA Thermophysics Conference, p. 3783, 2003. BIBLIOGRAPHY 136 [140] Y. Hu, S. Chen, and Q. Sun, “Hypersonic aerodynamics of a flat plate: Bridging formula and wall temperature effects,” in AIP Conference Proceedings, vol. 1501, pp. 1493–1499, 2012. [141] J. Wang, D. Han, and Y. Hu, “Kinetic simulations of plasma plume potential in a vacuum chamber,” IEEE Transactions on Plasma Science, vol. 43, no. 9, pp. 3047–3053, 2015. [142] J. E. Allen and M. Perego, “On the ion front of a plasma expanding into a vacuum,” Phys. Plasmas, vol. 21, no. 3, p. 034504, 2014. [143] M. Pfeiffer, S. Copplestone, T. Binder, S. Fasoulas, and C.-D. Munz, “Comparison of plasma plume expansion simulations using fully kinetic electron treatment and electron fluid models,” in AIP Conference Proceedings, vol. 1786, p. 130005, 2016. [144] M. Pfeiffer, A. Mirza, C.-D. Munz, and S. Fasoulas, “Two statistical particle split and merge methods for particle-in-cell codes,” Computer Physics Communications, vol. 191, pp. 9–24, 2015. [145] M. Vranic, T. Grismayer, J. L. Martins, R. A. Fonseca, and L. O. Silva, “Particle merging algorithm for pic codes,” Computer Physics Communications, vol. 191, pp. 65–73, 2015. [146] L. Brieda, D. VanGilder, and J. Wang, “Modeling ion beam neutralization and near- thruster plume interactions,” in 29th International Electric Propulsion Conference, (Princeton, NJ, USA), IEPC-2005-270, October 2005. [147] T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer Series in Statistics, Springer, 2 ed., 2009. [148] O. Buneman, “Maintenance of equilibrium by instabilities,” Journal of Nuclear Energy. Part C, Plasma Physics, Accelerators, Thermonuclear Research, vol. 2, no. 1, p. 119, 1961. [149] A. V. Gurevich, L. V. Pariiskaya, and L. P. Pitaevskii, “Ion acceleration upon expansion of a rarefield plasma,” Soviet Physics JETP, vol. 36, p. 274, 1973. [150] A. V. Gurevich and A. P. Meshcherkin, “Ion acceleration in an expanding plasma,” Soviet Physics JETP, vol. 53, pp. 937–945, 1981. [151] N. Singh and R. W. Schunk, “Numerical calculations relevant to the initial expansion of the polar wind,” J. Geophys. Res., vol. 87, no. A11, pp. 9154–9170, 1982. [152] V. Y. Bychenkov, V. N. Novikov, D. Batani, V. T. Tikhonchuk, and S. G. Bochkarev, “Ion acceleration in expanding multispecies plasmas,” Phys. Plasmas, vol. 11, no. 6, pp. 3242–3250, 2004.

Abstract (if available)

Abstract The collisionless, mesothermal plasma flow dynamics is a generic plasma physics problem having many potential applications. In space engineering, the plasma-spacecraft (S/C) interactions and the expansion of a plume from an electric propulsion (EP) thruster are two subjects in which the collisionless, mesothermal plasma flow dynamics plays an very important role. To understand the electron kinetics and validate the commonly used electron fluid approximations in such problems, this dissertation carries out comprehensive fully kinetic particle-in-cell (PIC) simulations as well as hybrid fluid-electron/particle-ion PIC simulations. ❧ The plasma-S/C interactions are first studied by considering a collisionless, mesothermal plasma flow over a large unbiased conducting plate. The hybrid PIC modeling first employs the most widely used Boltzmann relation, in which the electron is approximated by a massless, isothermal fluid, to describe the electron dynamics. The validity of this model is investigated by directly comparing the results between the hybrid and full PIC simulations. The comparison shows that the plasma wake may be characterized into two regions based on electron behavior, a fluid expansion region and a kinetic electron expansion region. In the fluid electron expansion region, the approaches based on the electron fluid approximation and the fully kinetic model lead to a similar result. In the kinetic electron expansion region, the fluid approximation for electrons breaks down and the fully kinetic approach must be adopted. The results also show that the error of the hybrid PIC simulation strongly correlates with a new non-equilibrium parameter (Ψ) based on the weighted deviation of the actual electron velocity distribution from the one at the equilibrium. The capability of the hybrid PIC modeling that incorporates the more general polytropic electron fluid approximation is then examined. The results predicted by this model exhibit very little improvement over the Boltzmann electron based hybrid PIC method. ❧ The influence of surface charging on the plasma-plate interactions is investigated next. In particular, the plasma wake structures, ion impact and current collection on the plate are studied with the Boltzmann electron based hybrid PIC and the full PIC simulations. The hybrid PIC modeling is shown to accurately describe the ion dynamics in the vicinity of the plate, but to fail downstream. The hybrid PIC failure is primarily due to the breakdown of electron fluid approximation, confirmed by the local non-equilibrium degree of electron flow measured by the non-equilibrium parameter Ψ. ❧ The EP plume problem is studied first with a 2D plume injected from a pre-neutralized source by using both the hybrid and full PIC simulations. The full PIC simulations are carried out using the real ion-to-electron mass ratios of H⁺, Ar⁺ and Xe⁺, while the hybrid PIC model employ the Boltzmann relation for the electrons. The simulations reveal that the plume structure exhibits four distinct expansion regions. These regions are identified as, from the plasma emission plane to the downstream along the plume direction, the unperturbed, quasi-steady expansion, self-similar expansion and electron front regions, respectively. The behavior of electrons is strongly anisotropic, causing considerably different expansion characteristics between the plume direction and the transverse direction. The 3D effects on the pre-neutralized plume expansion is further examined. The results show that the 2D plume is qualitatively similar to the 3D plume, but quantitative differences exist. The comparison between the hybrid and full PIC results shows that modeling electrons either as an isothermal fluid or as a polytropic fluid with an isotropic temperature is over-simplified and leads to not only quantitative but also qualitative differences. In order for an electron fluid approximation to be applicable, the model must at least have the capability of resolving the strong anisotropic thermodynamics for electrons. ❧ Finally, a fully kinetic PIC simulation is carried out to study the emission of a 2D collisionless plasma plume. All the important processes, including the ion-electron coupling, beam neutralization, and plume expansion, relevant to the beam plasma emission from an ion thruster are taken into account in this model. The macroscopic plume structure is shown to have several distinctive regions, including an undisturbed core region, an electron cooling expansion region, and an electron isothermal expansion region. The properties of each region are determined by the microscopic electron kinetic characteristics.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Grid-based Vlasov method for kinetic plasma simulations

PDF

Numerical and experimental investigations of ionic electrospray thruster plume

PDF

Particle-in-cell simulations of kinetic-scale turbulence in the solar wind

PDF

Three-dimensional kinetic simulations of whistler turbulence in solar wind on parallel supercomputers

PDF

Laboratory investigations of the near surface plasma field and charging at the lunar terminator

PDF

Experimental investigation on dusty surface charging in plasma

PDF

Numerical and experimental investigations of dust-plasma-asteroid interactions

PDF

The cathode plasma simulation

PDF

The space environment near the Sun and its effects on Parker Solar Probe spacecraft and FIELDS instrumentation

PDF

Experimental and numerical investigations of charging interactions of a dusty surface in space plasma

PDF

Particle simulations of colloid thruster beam

PDF

An experimental investigation of cathode erosion in high current magnetoplasmadynamic arc discharges

PDF

Dynamic modeling and simulation of flapping-wing micro air vehicles

PDF

An investigation into magnetically induced extractor-less electrospray propulsion devices

PDF

Advanced simulation models and concepts for positron and electron acceleration in plasma wakefield accelerators

PDF

A Boltzmann model for tracking aerosols in turbulent flows

PDF

Development and applications of a body-force propulsor model for high-fidelity CFD

PDF

Interactions of planetary surfaces with space environments and their effects on volatile formation and transport: atomic scale simulations

PDF

Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

Asset Metadata

Creator Hu, Yuan (author)

Core Title Kinetic studies of collisionless mesothermal plasma flow dynamics

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Astronautical Engineering

Publication Date 07/27/2018

Defense Date 06/19/2018

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag electric propulsion,kinetic theory,non-equilibrium plasma physics,OAI-PMH Harvest,particle-in-cell (PIC),plasma-spacecraft interaction

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Wang, Joseph J. (committee chair), Erwin, Daniel A. (committee member), Nakano, Aiichiro (committee member)

Creator Email yuanhu@usc.edu,yuanhu09@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-37382

Unique identifier UC11671923

Identifier etd-HuYuan-6543.pdf (filename),usctheses-c89-37382 (legacy record id)

Legacy Identifier etd-HuYuan-6543.pdf

Dmrecord 37382

Document Type Dissertation

Format application/pdf (imt)

Rights Hu, Yuan

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA