University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content THREE-DIMENSIONAL KINETIC SIMULATIONS OF WHISTLER TURBULENCE IN SOLAR WIND ON PARALLEL SUPERCOMPUTERS by Ouliang Chang A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ASTRONAUTICAL ENGINEERING) December 2013 Copyright 2013 Ouliang Chang Abstract The objective of this dissertation is to study the physics of whistler turbulence evolution and its role in energy transport and dissipation in the solar wind plasmas through computational and theoretical investigations. This dissertation presents the first fully three-dimensional (3D) particle-in-cell (PIC) simulations of whistler turbulence forward cascade in a homo- geneous, collisionless plasma with a uniform background magnetic fieldB o , and the first 3D PIC simulation of whistler turbulence with both forward and inverse cascades. Such computationally demanding research is made possible through the use of massively paral- lel, high performance electromagnetic PIC simulations on state-of-the-art supercomputers. Simulations are carried out to study characteristic properties of whistler turbulence under variable solar wind fluctuation amplitude ( e ) and electron beta ( e ), relative contri- butions to energy dissipation and electron heating in whistler turbulence from the quasilin- ear scenario and the intermittency scenario, and whistler turbulence preferential cascading direction and wavevector anisotropy. The 3D simulations of whistler turbulence exhibit a forward cascade of fluctuations into broadband, anisotropic, turbulent spectrum at shorter wavelengths with wavevectors pref- erentially quasi-perpendicular toB o . The overall electron heating yieldsT k >T ? for all e and e values, indicating the primary linear wave-particle interaction is Landau damping. But linear wave-particle interactions play a minor role in shaping the wavevector spectrum, whereas nonlinear wave-wave interactions are overall stronger and faster processes, and ultimately determine the wavevector anisotropy. ii Simulated magnetic energy spectra as function of wavenumber show a spectral break to steeper slopes, which scales ask ? e ' 1 independent of e values, where e is electron inertial length, qualitatively similar to solar wind observations. Specific spectral indices from simulated wavevector energy spectra do not match the frequency spectral indices from observations due to the inapplicability of Taylor’s hypothesis. In contrast, the direct comparison of simulated frequency energy spectra and solar wind observations shows a closer similarity. Electron density fluctuations power spectra also exhibit a close similarity to solar wind observations and MHD predications, both qualitatively and quantitatively. Linear damping represents an intermediate fraction of the total dissipation of whistler turbulence over a wide range of e and e . The relative importance of linear damping by comparison to nonlinear dissipation increases with increasing e but decreases with increasing e . Correlation coefficient calculations imply that the nonlinear dissipation pro- cesses in our simulation are primarily associated with dissipation in regions of intermittent current sheet structures. The simulation results suggest that whistler fluctuations could be the substantial constituent of solar wind turbulence at higher frequencies and short wave- lengths, and support the magnetosonic-whistler interpretation of the quasilinear scenario. An even larger scale 3D whistler turbulence simulation exhibits both a forward cascade to shorter wavelengths with wavevectors preferentiallyk ? >k k , and an inverse cascade to longer wavelengths with wavevectorsk k & k ? . The inverse cascade process is primarily driven by the nonlinear wave-wave interaction. It is shown that the energy inverse cas- cade rate is similar to the energy forward cascade rate at early times although the overall energy in the two cascades is very different. The presence of inverse cascade process does not affect qualitative conclusions established from the whistler turbulence forward cascade simulations. iii Dedication ... to my family and most dear iv Acknowledgments First, I would like to express my gratitude to my Ph.D. advisor Dr. Joseph Wang, who picked me from hundreds of applicants six year ago and give me a chance to fulfill my dream. I believe this dissertation is perfect example for his insight in choosing talents. I learned a lot from his persistence in research, his broad knowledge on space science and astronautical engineering, and most importantly, the computational method of particle-in- cell technique. Without his direction and support, I would not be able to complete this work. We have a very good work relationship, and I enjoy working with him. I hope we could continue to corporate after my graduation. I would like to express my special gratitude to Dr. S. Peter Gary, who serves as my mentor and collaborator for the past three years on this work. His valuable advice, guid- ance and encouragement help me a lot in finishing this work. I am really impressed by his enthusiasm and passion in solving plasma physics problem. His research attitude of focusing on details and fundamental physics help us overcome many problems. The way that he like writing papers and documenting every new discoveries inspire me a lot. I really appreciate him for encouraging me and eventually forcing me to write. We have argued on interpretations of many new results, but I really enjoying the process in such a collective way to do plasma physics research. I also hope we could continue to corporate in future. I would like to thank Dr. Aiichiro Nakano for leading me to the world of high perfor- mance computing and giving me general computational advices. I thank him for providing me his space filling curve routines. I have learned a lot from his courses and am able to v apply the knowledge into my actual research work. I also like to thank him for creating the dual-degree program in high performance computing and simulations, which not only helps me during my graduate study, but also would bring me a lot of opportunities in my future professional career. I would like to thank Dr. Mike Gruntman for the continued inspiration and advice, and his passionate lectures. This department and my role within it could not have existed without him. I would also like to thank Dr. Kaijun Liu and Dr. Shinji Saito for their suggestions on the turbulence spectrum initialization in the simulations. I additionally thank members of Dr. Wang’s research group, Dr. Ning Ding, Dr. John Polansky, Daoru Han, Will Yu, Kevin Chou and R. Scott Hughes for their support in my graduate study in USC. I also like to thank Thata Suksila, Dr. Doug Codron, Vaibhav Gupta, Sam Barbour and all other current Astronautical Engineering graduate students to make this department better and my life more joyfully. I additionally thank Dell Cuason, Marrietta Penoliar, Ana Olivares for making this department a more desirable place to stay and help my work and life on a daily basis. I am sincerely thankfully and indebted to those closest to me. I thank my parents and grandparents for their sacrifice and love, and their spiritually comfort for the last four years. I thank my sibling, especially Mike, for having great time together both in the United States and China. Finally and most importantly, I am deeply indebted to my girl friend, Viviana Huan Wei, for her love in the past 10 years since we met in college. I could not have done this without you, and I wish to make greater things with you in the next chapter of our life. This work is supported by National Science Foundation Grant AGS-1202603 and fund- ing from USC Viterbi School of Engineering. Additionally, this work benefited from my research performed during the 2011 Los Alamos Space Weather Summer School, sup- ported by the Institute of Geophysics and Planetary Physics, the Science, Technology and Engineering Directorate and the Global Security Directorate at Los Alamos National Lab- oratory. Computational resources supporting this work were provided by the National vi Aeronautics and Space Administration (NASA) High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center, by the Yellowstone supercomputer provided by National Center for Atmospheric Research (NCAR)’s Computational and Information Systems Laboratory (CISL) sponsored by the National Science Foundation, as well as by the USC High-Performance Computing and Communications (HPCC). Ouliang Chang August 1, 2013 vii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Background overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 Whistler turbulence and its role in short-wavelength solar wind tur- bulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Energy dissipation at electron scales . . . . . . . . . . . . . . . . . 8 1.2 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.1 Recent simulations of short-wavelength solar wind turbulence . . . 13 1.2.2 Recent observations of short-wavelength solar wind turbulence . . 18 1.3 Motivation and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.4 Dissertation outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Chapter 2: Computational Approach . . . . . . . . . . . . . . . . . . . . . . . . 25 2.1 Linear Vlasov dispersion solver . . . . . . . . . . . . . . . . . . . . . . . 25 2.2 Electromagnetic particle-in-cell modeling . . . . . . . . . . . . . . . . . . 27 2.2.1 Electromagnetic finite-difference time-domain method . . . . . . . 27 2.2.2 Parallelism on distributed-memory architecture . . . . . . . . . . . 29 2.2.3 Parallelism on hybrid shared- and distributed-memory architecture . 31 2.2.4 Adaptive particle sorting . . . . . . . . . . . . . . . . . . . . . . . 33 2.2.5 Parallel I/O and fast Fourier transform . . . . . . . . . . . . . . . . 37 2.3 Numerical techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.3.1 Turbulent spectrum initialization . . . . . . . . . . . . . . . . . . . 39 2.3.2 Spectral filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 viii 2.3.3 Wavelet transform . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.4 Computational requirements and facilities . . . . . . . . . . . . . . . . . . 43 2.4.1 Computational challenge . . . . . . . . . . . . . . . . . . . . . . . 43 2.4.2 State-of-the-art supercomputing facilities . . . . . . . . . . . . . . 44 Chapter 3: Three-Dimensional Particle-in-Cell Simulations of Whistler Turbulence Forward Cascade: Comparison with Two-Dimensional Simulations . . 45 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2 Initial conditions and simulation setup . . . . . . . . . . . . . . . . . . . . 46 3.3 Particle-in-cell simulations and discussions . . . . . . . . . . . . . . . . . 49 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Chapter 4: Effects of Initial Fluctuation Amplitude on Whistler Turbulence For- ward Cascade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.3 Particle-in-cell simulations and discussions . . . . . . . . . . . . . . . . . 63 4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Chapter 5: Effects of Electron Beta on Whistler Turbulence Forward Cascade . . 74 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.2 Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.3 Linear analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.4 Particle-in-cell simulations and discussions . . . . . . . . . . . . . . . . . 82 5.5 Electron-scale spectral break scaling . . . . . . . . . . . . . . . . . . . . . 106 5.6 Wave-wave coupling and the forward cascade . . . . . . . . . . . . . . . . 109 5.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Chapter 6: Energy Dissipation by Whistler Turbulence Forward Cascade . . . . . 118 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.2 Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.3 Particle-in-cell simulations and discussion . . . . . . . . . . . . . . . . . . 121 6.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 Chapter 7: Three-Dimensional Particle-in-Cell Simulations of Whistler Turbulence: Forward Cascade vs. Inverse Cascade . . . . . . . . . . . . . . . . . . 137 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 7.2 Initial conditions and simulation setup . . . . . . . . . . . . . . . . . . . . 139 7.3 Particle-in-cell simulations and discussions . . . . . . . . . . . . . . . . . 140 7.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 Chapter 8: Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . 153 8.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 8.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 ix Chapter A: Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 A.1 Background in plasma waves and instabilities . . . . . . . . . . . . . . . . 161 A.2 3D-EMPIC code sequential algorithm and normalization . . . . . . . . . . 163 A.3 Some definitions used in turbulence diagnostics . . . . . . . . . . . . . . . 164 Reference List . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 x List of Symbols ! angular wave frequency . . . . . . . . . . . . . . . . . . . . . . . . . 2 e electron cyclotron frequency . . . . . . . . . . . . . . . . . . . . . . 2 lh lower hybrid frequency . . . . . . . . . . . . . . . . . . . . . . . . . 2 p proton cyclotron frequency . . . . . . . . . . . . . . . . . . . . . . . 2 k wavenumber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 c speed of light . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 ! r the real part of wave frequency . . . . . . . . . . . . . . . . . . . . . 2 ! e electron plasma frequency . . . . . . . . . . . . . . . . . . . . . . . 2 e electron beta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 B o background magnetic field . . . . . . . . . . . . . . . . . . . . . . . 2 angle of mode propagation . . . . . . . . . . . . . . . . . . . . . . . 2 p proton inertial length . . . . . . . . . . . . . . . . . . . . . . . . . . 3 ! p proton plasma frequency . . . . . . . . . . . . . . . . . . . . . . . . 3 f wave frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 spectral index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 MHD magnetohydrodynamic . . . . . . . . . . . . . . . . . . . . . . . . . 4 AU astronomical unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 jBj 2 magnetic fluctuation energy . . . . . . . . . . . . . . . . . . . . . . 5 KAW kinetic Alfv´ en wave . . . . . . . . . . . . . . . . . . . . . . . . . . 5 EMHD electron magnetohydrodynamic . . . . . . . . . . . . . . . . . . . . 7 xi 2D two-dimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 PIC particle-in-cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 e electron inertial length . . . . . . . . . . . . . . . . . . . . . . . . . 8 k wavevector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 PDF probability density function . . . . . . . . . . . . . . . . . . . . . . 12 B i magnetic fluctuation field ini direction . . . . . . . . . . . . . . . . 12 P power per unit volume . . . . . . . . . . . . . . . . . . . . . . . . . 12 J current density vector . . . . . . . . . . . . . . . . . . . . . . . . . . 12 ? perpendicular to the background magnetic fieldB o . . . . . . . . . . 13 k parallel to the background magnetic fieldB o . . . . . . . . . . . . . 13 3D three-dimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 e electron gyroradius . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 v te electron thermal speed . . . . . . . . . . . . . . . . . . . . . . . . . 18 V sw solar wind velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 v phase velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 De Debye length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 v teo initial electron thermal speed . . . . . . . . . . . . . . . . . . . . . . 46 o dimensionless fluctuation energy by normalizing toB 2 o =8 . . . . . . 46 E(k ? ) k ? magnetic fluctuation energy spectrum . . . . . . . . . . . . . . . 47 B average wavevector anisotropy angle for forward cascade . . . . . . . 47 t Simulation time step . . . . . . . . . . . . . . . . . . . . . . . . . . 47 growth or damping rate . . . . . . . . . . . . . . . . . . . . . . . . . 63 e dimensionless fluctuation energy by normalizing ton e k B T e;t=0 . . . . 76 K e electron kinetic energy density . . . . . . . . . . . . . . . . . . . . . 77 PSD(k ? ) k ? power spectrum of electron density fluctuations . . . . . . . . . . 92 ^ e unit vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 r j B i fluctuating magnetic field increment . . . . . . . . . . . . . . . . . . 121 xii <= 4 ij > kurtosis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 inv average wavevector anisotropy angle for inverse cascade . . . . . . . 144 E(jkj) jkj magnetic fluctuation energy spectrum . . . . . . . . . . . . . . . 144 xiii List of Tables 2.1 3D-EMPIC parallel efficiency for some benchmark cases . . . . . . . . . . 31 2.2 Parallel I/O efficiency on NCAR’s GPFS file system . . . . . . . . . . . . . 37 3.1 Plasma parameters in the simulations . . . . . . . . . . . . . . . . . . . . . 48 3.2 Computational parameters in the simulations . . . . . . . . . . . . . . . . . 49 5.1 Different plasma parameters in three variable e simulations . . . . . . . . 76 5.2 A detailed listing of wave mode data from the e = 0:1 and e = 2:0 simulation that satisfies the wave-wave coupling Eqs. 1.7 and 1.8. All data are picked from the simulation as illustrated in the middle panel of Fig. 5.15. 104 5.3 Computational parameter difference in the new run simulation . . . . . . . 106 6.1 Correlation coefficients between n e T e and J 2 . All grid points in the 3D simulation space are taken into account. The largest coefficients are listed here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 7.1 Computational parameter difference in Run 1 and Run 2 simulations . . . . 139 xiv List of Figures 1.1 Solar wind observations of short-wavelength magnetic fluctuation spectra from Cluster spacecrafts’ magnetometers are fit with both two power laws frequency dependance in the left panel [76] and an exponential wavenum- ber dependence in the right panel [3]. . . . . . . . . . . . . . . . . . . . . . 4 1.2 (a) Linear theory results of kinetic Alfv´ en wave and whistler wave from [38]. (b) Left panel: Solar wind observations from four Cluster spacecrafts show fluctuations (green diamond symbols) with ! r > p at k p > 1 consistent with whistler mode dispersions [62]. Right panel: Solar wind observations from four Cluster spacecrafts show all fluctuations with! r p consistent with kinetic Alfv´ en wave dispersions [76] . . . . . . . . . . . 6 1.3 2D PIC simulation results of whistler turbulence from [80]. (a) The mag- netic energy spectrum is anisotropic, withk ? k k . (b)jB(k)j 2 steepen with 3<< 4 atjk ? e j 1 . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.4 (a) 2D configuration used in [79, 80] and pseudo-3D configuration used in [32]. (b) Fully 3D simulation configuration used in this dissertation research 14 1.5 2D PIC simulation results with a pseudo-3D configuration from [32]. (a) Frequency spectra of the wave magnetic field forB o at = 60 . The initial whistler decays into a daughter whistler by radiating a lower frequency magnetosonic / lower hybrid wave. (b) The magnetic wave vector spectrum for forB o at = 78 atj! e jt = 1000 andj! e jt = 2000 indicates that the mother wave experienced a large angle scattering. . . . . . . . . . . . . . . 15 1.6 2D PIC simulation results of sheared flow driven turbulence from [46]. (a) Current sheet formation: a contour plot ofjJj. (b) Frequency vsk y spec- trum of magnetic fluctuations computed at the edge of the simulation away from the vortex. Superimposed on the spectrum are lines corresponding to dispersion of compressionalk(U 0 V A ) and shear Alfv´ en modeskU 0 k k V A . 16 xv 1.7 3D gyrokinetic simulation results of antenna driven kinetic Alfv´ en turbu- lence from [93]. Collisionless damping via the Landau resonance with the electrons is sufficient to account for the measured heating as a function of scale in the simulation, without the need for significant Ohmic dissipa- tion. Measured heating of the electrons from the simulation by perpendic- ular wavenumber, Q e (k ? ) (black solid), an estimate of the electron heat- ing based on linear wave-particle dampingQ wp (k ? ) (blue dotted), and the Ohmic heating rate Q (k ? ) (red dashed). Landau damping sufficient to account for electron heating in simulation. . . . . . . . . . . . . . . . . . . 17 1.8 (a) Power spectra of electron density fluctuations (PSD) from calibrated data with artificial spikes removed and smoothed [23]. (b) Normalized transverse and parallel probability density functions of fluctuations in the inertial and dissipation ranges with Poisonnian error bars attached [47]. . . 20 2.1 The Yee lattice shows the location of the field components on the staggered, finite-difference mesh, where grid points are at the corners of the cell [107]. 28 2.2 Iteration flow chart for the parallel 3D-EMPIC. Routines in the rounded boxes require no inter-processor communication, whereas routines in the rectangular boxes involve communication [100]. . . . . . . . . . . . . . . . 30 2.3 (a) The physical volume is spatially decomposed into processes,P n . Each process consists of a number of blocks of computational cells in a pseudo- linked-list order (see section 2.2.4 for definition). Cells inside each block are traversed concurrently by threads (denoted by red dots with arrows) to compute the particle motion and solve the electromagnetic field. (b) Comparison of scalability of the hybrid parallel implementation and the purely distributed parallel implementation. All cases are performed on dual quad-core AMD Opteron 2:3 GHz processor nodes. For hybrid cases, 8 threads are spawned mapping 8 physical cores of a cluster node. . . . . . . 32 2.4 Comparison of simulation time per step for cases with and without the sort- ing mechanism, and with different space-filling curve sorting orders. . . . . 36 2.5 Cold plasma dispersion curves,!(k)=! p = cos 1=2 kx, showing compen- sation and smoothing fora 1 = 0:5 and variousa 2 values. The dashed line is the cold plasma dispersion for no compensation and no attenuation [7]. . 41 3.1 Upper panels: Reduced magnetic fluctuation energy wavevector spectra at t = 0: (left panel) Run 2 in the plane perpendicular toB o , (center panel) Run 3 in the plane perpendicular toB o , and (right panel) Run 3 in thek y -k z plane containingB o . Lower panel: 3D magnetic fluctuation energy spectra of Run 3 att = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.2 Perpendicular magneticB ? and electricE ? field polarization as a function of time from all three runs. . . . . . . . . . . . . . . . . . . . . . . . . . . 52 xvi 3.3 Reduced magnetic fluctuation energy wavevector spectra atj e jt = 447 from Run 1 (upper left panel), Run 2 (center panels) and Run 3 (right pan- els). The upper panels represents spectra in the k y -k z plane, whereas the lower panels represents spectra in the plane perpendicular to B o . Lower left panel is the corresponding 3D magnetic fluctuation energy spectra of Run 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.4 k ? reduced magnetic fluctuation energy spectra atj e jt = 447 from the 2D Run 1 (red line), and the 3D Run 2 (green line) and Run 3 (blue line) simulations. The blue dashed lines represent power law functions as labeled for comparison against the simulation results. . . . . . . . . . . . . . . . . 54 3.5 Time histories of the spectral wavevector anisotropy factor tan 2 B from the 2D Run 1 (red line) and the 3D Run 2 (green line) and Run 3 (blue line) simulations. Evaluations here are taken on spectra over 0:65jk e j 3:0. 56 3.6 Time histories of the (a) total fluctuating magnetic energy density and (b) the parallel (solid lines) and perpendicular (dashed lines) components of the total electron kinetic energy a function of time, for the 2D Run 1 (red line) and the 3D Run 2 (green line) and Run 3 (blue line), respectively. . . . 57 3.7 Reduced electron velocity distributions as functions ofv k (top panels) and v ? (bottom panels) for Run 1 (left panels) and Run 3 (right panels). The grey lines represent initial velocity distributions at t = 0, values, the red and blue solid lines correspond to Run 1 and Run 3 atj e jt = 447, and the dashed lines represent the difference of the two distributions normalized by the initial value of that distribution. . . . . . . . . . . . . . . . . . . . . . . 58 4.1 Simulation time histories of (a) the total fluctuating magnetic energy den- sityjB(t)j 2 =jBj 2 t=0 , (b) the parallel (solid lines) and perpendicular (dashed lines) components of electron temperature T ke =T keo and T ?e =T ?eo where T eo T e;t=0 , (c) the normalized magnetic energy damping ratedjBj 2 =dt=jBj 2 and (d) the normalized electron parallel temperature heating ratedT ke =dt=jBj 2 . Here the light blue lines represent the o = 0:02 run, red lines repre- sent the o = 0:05 run, dark blue lines represent the o = 0:1 run, green lines represent the o = 0:2 run, black lines represent the o = 0:5 run, and magenta lines represent the o = 1:0 run. The dashed line in (c) is djBj 2 =dt=jBj 2 =0:004. . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.2 Simulation time histories of (a) the average wavevector anisotropy factor tan 2 B , and (b) the average electron temperature anisotropyT ke =T ?e from the various o runs. The colors correspond to values of o are as labeled. . . 65 4.3 k z k y cuts from 3D wavevector iso-surfaces of magnetic fluctuation energy density at four times as labeled during the simulation of o = 0.2. . . . . . . 67 4.4 Time series of 3D wavevector iso-surfaces of magnetic fluctuation energy density as labeled during the simulations of o = 0:05; 0:1 and 0.2. . . . . . 68 xvii 4.5 k ? reduced magnetic fluctuation energy spectra atj e jt = 134 from the various o runs. The dashed lines represent power law functions as labeled for comparison against the simulation results. The colors correspond to values of o are as labeled. . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.6 Reduced electron velocity distributions as functions ofv k (top panel) and v ? (bottom panel) atj e jt = 447 from the various o runs. The colors of the solid lines correspond to values of o as stated in the caption to Fig. 4.1, with the grey line representing the initial distributions att = 0. The dashed lines represent the difference between the late-time and the initial distribution for the case o = 0:1 normalized by the initial distribution. . . . 71 5.1 Linear theory results for whistler dispersion as a function of wavenumber jk e j for three propagation angles at three e as labeled. Green, blue, red lines denote = 0 o ; 55 o , and 80 o or 70 o for e = 1:0 cases. Left column shows dispersion solutions for real frequency! and damping rate . Solid, dotted, dash dotted, dashed lines represent solutions for! 2 e = 2 e = 2500, 50, 5 and 0.5. Right column shows dispersion solutions for vs.!. The values ! 2 e = 2 e correspond to those listed in Table 5.1. . . . . . . . . . . . . . . . . 78 5.2 2D contours for whistler dispersion as functions of wavenumberjk k j and jk ? j at e = 0.01, 0.1 and 1.0 as labeled. The left column’s color represents the real frequency!, and the right column’s color represents the damping rate absolute valuej j. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 5.3 Simulation results as functions of time for three runs with initial values of e as labeled. (a) Total energy density (dashed lines) and total magnetic fluctuation energy density (solid lines). (b) Cascaded magnetic fluctuation energy density (summed over 0:65jk e j 3:0) normalized to initial electron kinetic energy density. (c) An enlarged (b) in linear scale of both axes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.4 Time histories of the normalized magnetic energy damping ratedjBj 2 =dt=jBj 2 for the three simulation sets with initial values of e and e as labeled. The dashed lines in upper, middle and lower panels are djBj 2 =dt=jBj 2 = - 0.00002, -0.004, -0.02, respectively. . . . . . . . . . . . . . . . . . . . . . 85 5.5 Wavevector anisotropy factor of the cascaded magnetic fluctuation energy as a function of time during the three simulations with initial e = 0.01, 0.10, and 1.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.6 (a) Normalized electron kinetic energy density and (b) normalized electron kinetic energy anisotropy as functions of time during the three simulations with initial e = 0.01, 0.10, and 1.0. . . . . . . . . . . . . . . . . . . . . . 86 5.7 Reduced electron velocity distributions as functions ofv k (top panels) and ofv ? (bottom panels) at late times during the three simulations with initial e = 0.01, 0.10, and 1.0. The distributions are shown in linear scales in the left panels, and logarithmic scales in the right panels. . . . . . . . . . . . . 88 xviii 5.8 k k k ? reduced magnetic fluctuation energy spectra, log[jB(k)j 2 =8n e k B T e;t=0 ], at four times for the three simulations with initial values of e as labeled. . . 90 5.9 k ? reduced magnetic fluctuation energy spectra at corresponding early time from the three simulation sets with initial values of e and e as labeled. The blue and red dashed lines represent, respectively, the break model’s power law functions and the exponential model’s exponential rollover functions, as labeled for comparison against the simulation results. . . . . . . . . . . . 91 5.10 k ? reduced power spectrum of electron density fluctuations PSD(k ? ) at corresponding early time from the three simulation sets with initial values of e and e as labeled. The blue and red dashed lines represent, respec- tively, the break model’s power law functions and the exponential model’s exponential rollover functions, as labeled for comparison against the simu- lation results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.11 Magnetic energy frequency spectra with initial values of e = 2:0 and e as labeled. The solid blue, green and red lines represent the spectra at e = 0.01, 0.1 and 1.0., the dashed blue lines represent the power law functions as labeled for comparison against the simulation results, and the two vertical black dashed lines in each panel bound the frequency range of the initially imposed whistler fluctuations. The spectrum is averaged over 4096 series of data points computed over the full temporal evolution of each simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.12 Continuous wavelet transform of magnetic fluctuation energy as functions of both timet and real frequency! from a single data series in the simula- tion at three different values of e as labeled. The 2D color contours repre- sent the value of the coefficient log[CWT(B x ) 2 +CWT(B y ) 2 +CWT(B z ) 2 ], where CWT denotes continuous wavelet transform function. Here the mother wavelet chosen is Morlet wavelet for all three cases. The two parallel white dashed lines in each panel indicate the frequency range of the ini- tially imposed whistler fluctuations. . . . . . . . . . . . . . . . . . . . . . 98 5.13 Fluctuation dispersion (real frequency! as a function of parallel wavenum- berjk k j) computed over the full temporal evolution of each simulation at three different values of e as labeled. The color contours represent simula- tion results of magnetic fluctuation energy, log[jB(k k ;k y )j 2 =8n e k B T e;t=0 ], and the black dashed lines represent numerical solutions of the full kinetic dispersion equation for whistler fluctuations. Herejk ? e j = 0.61 for all three values of e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.14 3D iso-surfaces of magnetic fluctuation dispersion (real frequency ! as functions of parallel wavenumber k k and perpendicular wavenumber k ? ) computed over the full temporal evolution of the e = 0:1 simulation. The color contours represent simulation results of magnetic fluctuation energy, log[jB(k k ;k y )j 2 =8n e k B T e;t=0 ]. Each subfigure represents a different view perspective as labeled. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 xix 5.15 k k k ? reduced magnetic fluctuation energy spectra, log[jB(k k ;k y )j 2 =8n e k B T e;t=0 ], at four frequency bands as labeled for the three simulations with initial val- ues of e as labeled. These spectra, which are as functions of frequencies, are computed over the same time intervals as those used in Fig. 5.3(a) . . . 103 5.16 Wavevector anisotropy factor of the cascaded magnetic fluctuation energy as a function of frequency during the three simulations with initial e = 0.01, 0.10, and 1.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.17 k ? reduced magnetic fluctuation energy spectra at corresponding early time from the various simulations at e = 0:1 and e = 2:0. The blue dashed lines represent power law functions as labeled for comparison against the simulation results, and the vertical black dashed lines indicates the scaling as labeled. The left panel, normalized tojk ? e j, shows the spectra of the new run as defined before. The right panel, normalized to k ? e , shows the spectrum of the old run as defined before in the green sold line, the spectral filtered spectrum of the old run in the black solid line with the filter parameters as labeled, and the spectrum of the modified old run in the red sold line with only largest initial fluctuation modesk k=? e =0:1227 retained. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.18 Linear theory results for whistler dispersion properties as functions of wavenum- ber. (a) = 0 o . (b) = 75 o . Solid lines represent the real frequency ! and dotted lines show the damping rate . The blue, green and red lines correspond the three simulation sets in section 5.4 at e = 0.01, 0.1, and 1.0, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.19 Linear theory results for the frequency mismatch ! for whistler fluctu- ations over 0.20 kc=! e 1.0 as a function of propagation angle. The blue, green and red lines correspond the three simulation sets in section 5.4 at e = 0.01, 0.1, and 1.0, respectively. . . . . . . . . . . . . . . . . . . . . 113 6.1 Comparison of total damping rate t (solid lines) and linear damping rate l (dashed lines). (a) as functions of e andj e jt at e = 2:0; (b) as functions of e andj e jt at e = 0:1; (c) as functions of e and e at their respective early times, that is for e = 0:01 cases,j e jt = 707:1, for e = 0:1 cases, j e jt = 111:8, and for e = 1:0 cases,j e jt = 35:4. . . . . . . . . . . . . . 122 6.2 e = 0:01 cases: Parallel current density, J k , for a perpendicular plane (xy) at timej e jt = 707:1, with initial value of e and different band-pass filters applied as labeled. The color contours represent simulation values of J k = p K e;t=0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.3 e = 0:1 cases: Parallel current density,J k , for a perpendicular plane (x y) at timej e jt = 134:2, with initial value of e and different band-pass filters applied as labeled. The color contours represent simulation values of J k = p K e;t=0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 xx 6.4 e = 1:0 cases: Parallel current density,J k , for a perpendicular plane (x y) at timej e jt = 35:4, with initial value of e and different band-pass filters applied as labeled. The color contours represent simulation values of J k = p K e;t=0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 6.5 Current density square J 2 , electron temperature n e T e , and electric field energy E 2 from a perpendicular plane (xy). The color contours rep- resent normalized simulation values of corresponding quantities. The left column isJ 2 , the middle column isn e T e , and the left column isE 2 . The upper panels are from e = 0:1 and e = 2:0 case at timej e jt = 134:2, and the correlation coefficient ratio (neTe;J 2 ) (neTe;E 2 ) = 4:83. The lower panels are from e = 1:0 and e = 5:0 case at timej e jt = 14:1, and the correlation coefficient ratio (neTe;J 2 ) (neTe;E 2 ) = 1:92. Note that the correlation coefficients are computed from all grid points in the 3D space, not just the 2D planes shown here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 6.6 Normalized PDF of fluctuating magnetic field increments from e = 0:1 and e = 2:0 case at timej e jt = 134:2 for different spatial separation lengthr as labeled. (a) the perpendicular component ry B x ; (b) the parallel component ry B z . Here the PDFs are normalized to the standard deviation of their respective magnetic field increments. The black circled lines indicate the Gaussian distribution. . . . . . . . . . . . . . . . . . . . . . . 131 6.7 Kurtosis of the magnetic field increments r j B i as a function of spatial separation length r for different e and e as labeled. The left panel is from e = 0:01 cases at timej e jt = 707:1; The middle panel is from e = 0:1 cases at timej e jt = 134:12; The right panel is from e = 1:0 cases at timej e jt = 35:4. Data are low pass filtered atk x=y=z e = 5 to remove the statistical noise due to finite superparticle number of the PIC simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 7.1 Reduced magnetic fluctuation energy wavevector spectra, log(jBj 2 =8n e k B T e;t=0 ), from Run 2. The left column isk x k y spectra, the right column isk k k y spectra. Upper panel: att = 0; Lower panel: at atj e jt = 78:3. . . . . . . . 141 7.2 (a-b) Ak k k y cut atk x = 0 from 3D magnetic energy spectrum, log[jB(k k ;k y ;k x = 0)j 2 =8n e k B T e;t=0 ], at t = 0 andj e jt = 78:3 from Run 2. (c) The anisotropy factor tan 2 B and tan 2 inv of the cascaded energy as a func- tion of time. The green line is tan 2 B from Run 1. The black and blue lines are tan 2 B from the forward and inverse cascade of Run 2. The red line is tan 2 inv from Run 2 inverse cascade. Evaluations are taken on spectra over 0:65jk e j 3:0 for Run 1, 0:52jk e j 3:0 for forward cascade of Run 2, andjk e j 0:22 for inverse cascade of Run 2. . . . . . . . . . . . . 143 xxi 7.3 jkj reduced magnetic fluctuation energy spectraE(jkj) atj e jt = 67:1. The red and black lines are from Run 1 and Run 2, respectively. The green line is k ? reduced spectrumE(k ? ) from Run 1. The grey dashed line represents the initial spectrum of Run 2 at t = 0. The blue dashed lines represent power law functions as labeled for comparison against the simulation results.145 7.4 Magnetic fluctuation energy as functions of time from both Run 1 and Run 2 as labeled. (a) Total fluctuation energy density. (b) Cascaded fluctuation energy density normalized to initial electron kinetic energy density. (c) The rate of cascaded fluctuation energy density. The cascaded range is described in the caption of Fig. 7.2(c). . . . . . . . . . . . . . . . . . . . . 147 7.5 (a) Normalized electron kinetic energy density and (b) normalized electron kinetic energy anisotropy as functions of time from both Run 1 and Run 2 as labeled. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 7.6 Reduced electron velocity distributions as functions ofv k (top panels) and ofv ? (bottom panels) atj e jt = 67:1 from both Run 1 and Run 2 as labeled. The distributions are shown in linear scales in the left panels, and logarith- mic scales in the right panels. . . . . . . . . . . . . . . . . . . . . . . . . . 150 8.1 A schematic picture of a magnetic energy spectrum of short wavelength solar wind turbulence based on the magnetosonic-whistler interpretation under the framework of quasilinear scenario. . . . . . . . . . . . . . . . . . 159 xxii CHAPTER 1: INTRODUCTION Understanding the Sun and its various interactions in the solar system is a major goal in many space science research programs. For instance, a top science objective of NASA’s heliophysics program is to understand the fundamental physical processes of the space environment from the Sun to Earth, to other planets, and beyond to the interstellar medium [1]. This dissertation describes a research in the solar wind, to investigate the physics of one of fundamental plasma turbulence processes, whistler turbulence. Its overall objective is to understand the role of whistler turbulence in energy transport and dissipation in the solar wind. This research falls into the broad research area concerning the flow of energy and particles from the Sun and its influences throughout the interconnected Sun-heliosphere- Earth system. In this chapter, an overall background review is presented first. The concepts of space plasma turbulence and microinstability are introduced, followed by discussions of whistler turbulence’s role in short wavelength energy cascade and dissipation and the sig- nificance of whistler turbulence research. Energy dissipation mechanisms in electron-scale short-wavelength turbulence and current scenarios in describing plasma turbulence is then described. Next, we review recent research results of relevant theoretical, computational and observational studies of whistler turbulence evolution and solar wind energy dissipa- tion at electron scale wavelengths. Based on the previous works and present understanding of whistler turbulence and its relationship to energy flow in the solar wind, the motivation and the specific objectives of this dissertation research are presented. Finally, an outline of the work described in this dissertation is presented. 1 1.1 Background overview 1.1.1 Whistler turbulence and its role in short-wavelength solar wind turbulence Whistler waves are low frequency, right-hand circular / elliptical polarized electromagnetic waves, which not only exist in the Earth’s magnetospheric and ionospheric plasmas [28], but also are often observed in the solar wind [6, 51]. The whistler wave frequency resides between lh < ! < e , where! is the angular wave frequency, e is electron cyclotron frequency, and lh is the lower hybrid frequency defined as p e p , where p is the pro- ton cyclotron frequency. Whistlers’ propagation speed increase with the frequency, and since their frequency range could be easily translated into the audio band, the character- istic descending tone with time sounds like someone whistling and therefore gives rise to the name. The whistler wave dispersion relation is very complex and requires a numerical solution to the full kinetic linear dispersion equation. An approximate analytic expression using least-squares fits of an assumed form to the numerical results in the warm plasmas condition is given in Eq. 1.1 [80], ! r j e j = kk k c 2 =! 2 e 1 + (1 + p e )(k 2 c 2 =! 2 e ) (1.1) wherek is the scalar wavenumber,c is the speed of light,! r is real part of the frequency, ! e is the electron plasma frequency, e is electron beta, and the subscript k denotes the direction parallel to the background magnetic fieldB o . In the cold plasma approximation, the “classical quasi-parallel propagating whistler” dispersion equation is [33] ! r ' kc p ! p (1 + k 2 c 2 cos 2 ! 2 p ) 1=2 (1.2) where is the angle of mode propagation. 2 We define “plasma turbulence” as a broadband ensemble of incoherent field fluctuations in both frequency and wavenumber regimes that a continuous spectrum is present in an ion- ized medium. “whistler turbulence” is then a broadband ensemble of incoherent field fluc- tuations in a magnetized plasma at frequencies between the lower hybrid lh and electron cyclotron e frequencies and at wavelengths much shorter than the proton inertial length p . Note that the whistler wave linear dispersion of Eq. 1.1 are not necessarily complete to describe the whistler turbulence, and it is more appropriate to use fully nonlinear compu- tational approach in conjunction with linear dispersion theory results to study the whistler turbulence as we did in this dissertation. We also define “short-wavelength” as wavelengths at and shorter than the proton inertial length, p = c=! p , and “long-wavelength” as wave- lengths larger than p , where! p is the proton plasma frequency. In plasma turbulence, both the kinetic energy and the electromagnetic field energy are distributed over a broadband range of wavelengths, and the term “energy cascade” is used to describe energy transport between different scales of fluctuations. We define “forward cascade” to mean a trans- fer of fluctuating magnetic field energy from longer to shorter wavelengths, that is, from smallerk values to largerk values. Similarly, “inverse cascade” means a transfer of field energy from largerk to smallerk. Evolution of electromagnetic turbulence is a fundamen- tal plasma physics problem with implications to a wide range of areas, and therefore, is of great interest in solar wind plasmas as well as in magnetospheric and ionospheric plasmas. The fundamental picture of long-wavelength plasma turbulence in the solar wind fol- lows the classic Kolmogorov scenario of neutral fluid turbulence: a mechanism, such as a velocity shear instability, drives enhanced fluctuations at relatively long wavelengths. Through nonlinear interactions among these fluctuations, a forward cascade is established which transfers the energy of the oscillations to successively shorter wavelengths. How- ever, the solar wind is not a neutral fluid, but an almost collisionless, magnetized plasma which bears magnetic and electric field fluctuation spectra that exhibit turbulent properties in the observed frequency f. Spacecraft observations in the near-Earth solar wind show 3 (a) Sahraoui et al. 2010 (b) Alexandrova et al. 2012 Figure 1.1: Solar wind observations of short-wavelength magnetic fluctuation spectra from Cluster spacecrafts’ magnetometers are fit with both two power laws frequency dependance in the left panel [76] and an exponential wavenumber dependence in the right panel [3]. magnetic fluctuation energy spectra scaling approximately as f 5=3 from about 10 4 Hz down to about 0:2 0:5 Hz [9] as illustrated in Fig. 1.1; this is called the inertial range of solar wind turbulence, which corresponds to magnetic spectra observed at wavelengths long compared to p and frequencies smaller than the proton cyclotron frequency p . Iner- tial range spectra satisfy power-law dependences of spectral index'5=3 in frequencies and wavenumbers over several orders of magnitude in those variables. The forward cascade of the solar wind inertial range may be described by magneto-hydrodynamic (MHD) theory which treats the plasma as a magnetized fluid. In contrast to Kolmogorov isotropic neutral fluid turbulence, both solar wind observations and numerical simulations show MHD turbu- lence to be strongly anisotropic, with larger amplitude of magnetic fluctuations propagating perpendicular than parallel to the background magnetic fieldB o [54, 50, 86, 55, 42, 56]. At a frequency near 0.2 Hz 0.5 Hz, measurements at 1 astronomical unit (AU) show a spectral break, a distinct change to spectra that are of steeper power laws than those of the inertial range [50, 87, 2]. Recent analyses of data from the Cluster mission spacecraft have shown that magnetic spectra in the range 0.5 Hz < f < 50 Hz scale with 2.6 < < 4 2.8 [77, 76, 48, 4], and suggest that there is a second break at electron-scale with still more steeply decreasing spectra at approximately 50 Hz < f. We denote the first break as “proton-scale spectral break”, and the second as “electron-scale spectral break”. Mag- netic fluctuation spectra at frequencies above and wavelengths shorter than the proton-scale spectral break has been fit with different steeper power laws on two or three successively higher frequency ranges [76, 69] or an exponential wavenumber dependence [3], as shown in Fig. 1.1. Some of these observations [22, 76] demonstrated that turbulence in this regime is also anisotropic that has more fluctuation energy at propagation perpendicular toB o . Both of these spectral breaks imply failure of the MHD fluid description of turbulence at this distinct short wavelength and high frequency regime, and kinetic (i.e., velocity-space) physics must be used to describe plasma turbulence. Fluctuations above the first spectral break are relatively weak, such thatjBj 2 B 2 o , wherejBj 2 is the magnetic fluctua- tion energy. Turbulence in this regime is thus often described in terms of an ensemble of weakly interacting normal modes of the plasma which are individually described by lin- ear dispersion theory. In this so-called quasilinear framework, it is generally believed that the predominant constituents of solar wind turbulence at wavelengths immediately shorter than those of the proton-scale spectral break are likely to be kinetic Alfv´ en waves (KAWs) [50, 5, 76, 41, 81, 94], which propagate in directions quasi-perpendicular to B o and at frequencies! r < p , as shown in the left panel of Fig. 1.2(a). However, there is currently substantial debate about the properties of the very short wavelength turbulence that is far above the first proton-scale spectral break and extends down to and through the second electron-scale spectral break [88]. Within the quasilinear framework, there are two competing hypotheses as to the character of these fluctuations. One hypothesis is that KAWs remain the primary constituent of this electron-scale turbu- lence [44, 45, 77]. Recent solar wind observations [5, 77, 76] have been interpreted as consisting of KAWs. A fluid model of kinetic Alfv´ en turbulence, which does not include kinetic plasma effects, yields magnetic fluctuation energy spectra which scale ask 7=3 ? [43]. 5 (a) (b) Figure 1.2: (a) Linear theory results of kinetic Alfv´ en wave and whistler wave from [38]. (b) Left panel: Solar wind observations from four Cluster spacecrafts show fluctuations (green diamond symbols) with! r > p atk p > 1 consistent with whistler mode disper- sions [62]. Right panel: Solar wind observations from four Cluster spacecrafts show all fluctuations with! r p consistent with kinetic Alfv´ en wave dispersions [76] . A recent 3D gyrokinetic simulation of kinetic Alfv´ en turbulence has exhibited ak 2:8 ? mag- netic fluctuation spectrum similar to those of solar wind observations [45]. But other linear theory calculations and observations [38, 71, 70, 22] show that the angular propagation range of lightly damped KAWs becomes extremely narrow ( 1 ) atk p > 1, indicating that it is unlikely that KAWs could persist down to electron scale lengths, and thus not necessarily provide a complete description of short wavelength turbulence. 6 An alternate hypothesis is that some part of short-wavelength solar wind turbulence consists of magnetosonic-whistler fluctuations which propagate between p and the lower hybrid frequency lh , and whistler fluctuations which propagate between lh and e [40, 90, 49, 29, 37, 2, 85, 91, 20]. Recent observation [62] indicates many fluctuations are consistent with whistlers fluctuations. Simulations of whistler turbulence using 3D electron magnetohydrodynamic (EMHD) fluid models [10, 25, 26, 83] typically show a forward cas- cade to steep magnetic spectra withk 7=3 , and anisotropies preferring quasi-perpendicular propagation similar to those of the inertial range. Several recent two-dimensional (2D) particle-in-cell (PIC) simulations of whistler turbulence in low plasmas [37, 79, 80] have been in favor of this interpretation and demonstrated forward cascades of whistler fluc- tuations to broadband, relatively short wavelength turbulence with a preference to quasi- perpendicular propagation relative to B o and a very steep wavenumber dependences of magnetic spectra (e.g.,k 4 ? ). Observations of a pronounced quasi-perpendicular wavevector anisotropy [22, 76] are consistent with both the predictions of simulations for both kinetic Alfv´ en and whistler tur- bulence. However, an important difference between these two modes, involves dispersion relations. Fig. 1.2(a) compared the linear theory results of kinetic Alfv´ en and whistler waves. In the left panel of Fig. 1.2(b), Saharoui et al. [76] show all fluctuations with ! r p consistent with KAW dispersions, whereas in the right panel, Narita et al. [62] indicate many observations with! r > p atk p > 1, indicative of magnetosonic-whistler fluctuations. Resolving the difference between these two types of observations is critical in understanding the relative importance and contributions of whistler waves and KAWs to short-wavelength plasma turbulence. Therefore, to contribute to the knowledge of the character of short-wavelength solar wind turbulence, to enhance our understanding of the solar wind energy dissipation, which is one of fundamental processes that control our space environment and influence our 7 Earth’s atmosphere, this dissertation investigates the characteristics of whistler turbulence evolution and their effects on energy distribution in the homogeneous solar wind plasmas. Despite intense researches have been done to investigate the role of whistler waves in fast magnetic reconnection, relatively few studies focus on the whistler turbulence evolu- tion and its energy cascade in the solar wind. Whistler turbulence has been studied using both fluid models such as EMHD [30, 83] and kinetic models such as 2D PIC technique [37, 79, 80, 78, 32]. However, none of these provide a complete description of whistler turbulence characteristics in the an actual gyrotropic solar wind. EMHD models do not represent the velocity-space, or kinetic, properties of collisionless plasmas, thus lacking wave-particle interactions; 2D PIC simulations are intrinsically not gyrotropic which may lead to different whistler cascade process from those in realistic 3D systems. As a con- sequence, EMHD models generally yield magnetic spectra scaled as k 7=3 at k e < 1, whereas 2D PIC simulations give ak 4:0 ? , where e =c=! e is the electron inertial length. In contrast, 3D PIC simulations not only preserve plasma kinetic properties, but also represent the actual gyrotropic turbulence system, and therefore is the most accurate computational technique to study the full physics of whistler turbulence, and naturally become the major research approach in this dissertation. 1.1.2 Energy dissipation at electron scales Current discussions of short-wavelength plasma turbulence are often carried out in one of two distinct scenarios. In one viewpoint as we term the quasilinear scenario, at sufficiently small amplitudes, plasma turbulence may be treated as an ensemble of weakly interacting normal modes in a relatively homogeneous medium in which magnetic and electric field dissipation is described by linear damping theory, and plasma heating by quasilinear theory. Under the quasilinear scenario, the distribution and conversion of fluctuating field energy into plasma kinetic energy may be described in terms of the damping of linear plasma 8 modes (i.e. via linear wave-particle Landau and/or cyclotron damping), and by various higher order processes with or without dissipation within the weak turbulence theory (e.g., quasilinear theory, nonlinear Landau damping, wave-wave interaction). Linear Landau damping in unmagnetized plasmas is characterized by the expression ! =kv (1.3) where k is the wavevector. A particle is Landau resonant with a wave when its speed is equal to the phase speed of the wave in the direction of the particle’s velocity. In a magnetized plasma, there are two important linear dissipation mechanisms for whistler fluctuations. One is electron Landau damping, which occurs when parallel electric fields of obliquely propagating whistlers interact with electrons satisfying the resonance condition [80], !k k v ke = 0 (1.4) wherek k andv ke are the parallel wavenumber and the electron parallel velocity. Particles satisfying this condition move such thatEv ke does not change sign so that there can be a strong exchange of energy between the wave and the parallel motion of those particles. This resonance becomes strong when there are an appreciable number of electrons with parallel speed close to the parallel phase speed !=k k of whistlers. The associated wave-particle interaction scatters the electrons parallel to the background magnetic field, increasing their parallel temperature. The other is electron cyclotron damping, which is the magnetized analogue of Lan- dau damping. Physically, cyclotron damping occurs when the particle sees a wave whose Doppler shifted frequency is the cyclotron frequency or some harmonic thereof, !k k v ke =l e ; l =1;2; (1.5) 9 Particles satisfying this condition move such that Ev ?e does not change sign so that there is a strong exchange of energy between the wave and the perpendicular motion of those particles. The electrons propagating in the opposite direction to the whistler fluctua- tions can resonate with them because by definition whistler wave!<j e j. This resonance scatters the electrons perpendicular to the magnetic field, increasing or decreasing their perpendicular temperature. At sufficiently weak fluctuation amplitude, nonlinear wave-particle processes can be described by the quasilinear theory, which still involves the linear wave-particle interaction characterized by Eq. 1.3. In this theory, the wave amplitudes are still considered to be small enough that the wave propagation can be treated by the linear theory. The nonlinear part concerns the long-term effect of many waves on the background particles’ distribution functionf(v) [63]. Quasilinear effect thus describes diffusion of the distribution function f(v) in velocity space. It assumes that there is a sufficient spread in bandwidth in frequency and wavenumber, and its interaction involves a broader wave spectrum and is more of a random interaction between waves and particles. Another type of nonlinear wave-particle interaction is called nonlinear Landau damp- ing or induced scattering, which describes the nonlinear coupling of one wave to another wave through the background particles. In contrast to linear electron Landau damping in magnetized plasmas, the resonance is between two waves and one electron characterized by the expression ! 1 ! 2 = (k k1 k k2 )v ke (1.6) that is, the electron is resonant with the beat between two waves; its parallel speed is equal to the parallel phase speed of the beat between the two waves. Electrons could gain energy through the resonance with the beat and damp the waves just like the linear Landau damping mechanism. However, in the case that the second wave is a tiny wave with electric fieldE 2 compared to the first wave with a relatively large amplitudeE 1 , the second wave 10 could use the electrons to drain energy from the larger fieldE 1 and exciteE 2 ; hence the term induced scattering [63]. Since part of energy is transferred from waves to plasmas, the scattered second wave frequency! 2 is always slightly lower than the first wave! 1 , and a large change in the angle of propagation between wavevectorsk 1 andk 2 is possible [32]. There is a third important aspect of weak turbulence theory, which does not involve any wave-particle interaction, but is necessary for understanding whistler turbulence cas- cadezs, that is nonlinear three-wave interaction (generally called “wave-wave coupling / interaction”) [63]. Wave energy can be transferred via a three-wave interaction, a nonlin- ear process which must satisfy, the necessary but not sufficient, frequency and wavevec- tor matching conditions, representing conservation of fluctuation energy and momentum among the three interacting modes [103] ! 1 (k 1 )! 2 (k 2 ) =! 3 (k 3 ) (1.7) and k 1 k 2 =k 3 (1.8) where we assume!(k) is a solution of the linear kinetic dispersion equation, and indices 1 and 2 correspond to the two coupling wave modes, and 3 corresponds to the product wave mode. Note that this is a multidimensional set of equations with no preferred direction whereby forward and backward propagating waves at different angles can interact and sat- isfy these matching equations to enhance possible modes of higher or lower frequencies, shorter or longer wavelength. For wave-wave interactions, we assume that resonant parti- cles interactions (both linear and nonlinear) are not important, so that the conditions 1.7 and 1.8 do not involve a particle velocity, thus no plasma damping nor wave energy dissipation involved. The physical interaction of wave-wave coupling is the same as the physical inter- action that yields the parametric decay instability where the “parameter” is the amplitude of the single wave [63]. However, in that case, a single finite-amplitude wave decays into two 11 other waves. By contrast, wave-wave coupling in weak turbulence theory considers a broad spectrum of waves, and statistically predicts the evolution of the waves in an ensemble of system rather than in a single system. In section 5.4 and 5.6, we use these wave-wave cou- pling conditions to examine wavevector anisotropies for the forward cascade of whistler turbulence. The other viewpoint, which we call the intermittency scenario, is that, for sufficiently large amplitudes, plasma turbulence consists of fundamentally nonlinear, coherent struc- tures involving current sheet formation, which are the primary sites of magnetic field dis- sipation and plasma heating. In this picture the plasma is intrinsically inhomogeneous and intermittent, and conventional linear and quasilinear theories are inadequate for describ- ing the plasma dynamics. Recent solar wind research has demonstrated that magnetic field discontinuities can be observed all the way down to electron scales [68], and PIC simulations support this [97, 46, 20]. Analyses of intermittency can be framed in terms of probability density functions (PDFs) of the fluctuating magnetic fieldjB i j as a func- tion of spatial separation length, whereB i denotes magnetic fluctuation field ini direc- tion. At large separation scales, PDFs of short-wavelength inhomogeneous plasmas typi- cally exhibit Gaussian-like distributions, whereas as the scales decrease the PDFs develop enhanced “tails” indicating an excess of large amplitude values ofjB i j [46]. Cluster obser- vations by [2] show that these non-Gaussian tails increase as the frequency of the turbulent fluctuations increases. Kiyani et al. [48, 47] also demonstrate strongly non-Gaussian PDFs in the short-wavelength regime as measured by Cluster. The detailed kinetic physical mechanism by which dissipation occurs in current sheets has not been elucidated [93]. It has been suggested that the acceleration of electrons is due to interaction with parallel electric fieldE k associated with reconnections at current sheets [46]. Note that Joule / Ohmic heating does not exist in a collisionless plasma, andP =JE is dissipation rate for any dissipation processes in collisionless plasmas, whereP is power per unit volume,J is current density. 12 To identify the roles of and the relationships between the quasilinear scenario and the intermittency scenario in describing the dissipation of short-wavelength plasma turbulence is of fundamental importance in understanding the character of short-wavelength turbu- lence in the solar wind. 1.2 Literature review Theoretical, computational and observational studies of whistler and other short- wavelength kinetic turbulence relevant to the subject of this dissertation are reviewed in this section. 1.2.1 Recent simulations of short-wavelength solar wind turbulence This dissertation research is motivated by two sets of recent studies on nonlinear evolution of whistler turbulence / instability in the solar wind which use different approaches to initialize the simulation and focus on different aspects of plasma properties. The first set is a series of recent papers by [37, 79, 80] on cascading whistler turbulence. In these studies, 2D electromagnetic particle-in-cell (EMPIC) simulations were performed to study the energy cascade of whistler turbulence in a homogeneous plasma at relatively small amplitudes atjk e j 0:1 and electron plasma beta = 0:1=0:01. In these computa- tions, whistler turbulence evolved from an initial narrow-band spectrum of relatively long wavelength whistler modes. The results showed that the turbulence evolution is dominated by a forward cascade from long to short wavelengths, and that the magnetic energy spec- trum is anisotropic in which the fluctuation energy cascades preferentially perpendicular (?) than parallel (k to the background magnetic fieldB o , as illustrated in the upper panel Fig. 1.3(a). EMPIC simulations also showed that damping of the turbulence leads to heat- ing of electrons primarily in the directions parallel and antiparallel toB o , indicating that electron Landau damping at oblique propagation is the dominant dissipation factor in the 13 (a) (b) Figure 1.3: 2D PIC simulation results of whistler turbulence from [80]. (a) The magnetic energy spectrum is anisotropic, withk ? k k . (b)jB(k)j 2 steepen with 3 < < 4 at jk ? e j 1 Figure 1.4: (a) 2D configuration used in [79, 80] and pseudo-3D configuration used in [32]. (b) Fully 3D simulation configuration used in this dissertation research cascading whistler turbulence. Magnetic fluctuation spectrajB(k)j 2 steepen with a power law relation with4 atjk ? e j 1 (Fig. 1.3(b)). The second set is a recent study by [32] on the nonlinear evolution of instability- driven whistler fluctuations. The nonlinear theory developed in these studies suggested that 14 (a) (b) Figure 1.5: 2D PIC simulation results with a pseudo-3D configuration from [32]. (a) Frequency spectra of the wave magnetic field for B o at = 60 . The initial whistler decays into a daughter whistler by radiating a lower frequency magnetosonic / lower hybrid wave. (b) The magnetic wave vector spectrum for forB o at = 78 atj! e jt = 1000 and j! e jt = 2000 indicates that the mother wave experienced a large angle scattering. the whistler turbulence is fundamentally a nonlinear three-dimensional (3D) phenomena. Strong induced nonlinear scattering due to density perturbationn resulting from pondero- motive force can convert quasi-electrostatic waves into electromagnetic waves and vice versa, and thereby determines the character of the whistler turbulence. 2D EMPIC simu- lations were carried out for a pseudo-3D configuration (Fig. 1.4(a)), whereB o is placed at an inclination angle out of the 2D simulation plane to include the krn effect. Simulations triggered by a perpendicular ring velocity distribution instability show that the nonlinear scattering of whistler modes leads to inverse cascade toward lower frequencies and a large change in wavevectors propagation angle in the evolution of low e plasmas, as illustrated in panels Fig. 1.5(a) and (b), respectively. 15 (a) (b) Figure 1.6: 2D PIC simulation results of sheared flow driven turbulence from [46]. (a) Current sheet formation: a contour plot ofjJj. (b) Frequency vsk y spectrum of magnetic fluctuations computed at the edge of the simulation away from the vortex. Superimposed on the spectrum are lines corresponding to dispersion of compressional k(U 0 V A ) and shear Alfv´ en modeskU 0 k k V A . We note that it is difficult to compare the results based on the two set of studies from [37, 79, 80] and [32], since they use very different simulation setup and initial conditions. In [32], fluctuations are generated by initially loading a perpendicular ring distribution of “heavy electrons” in the plasma, whereas in [37, 79, 80], the turbulence is evolved from an ensemble of preloaded whistler modes. The initial wave modes in [32] are very anisotropic while that in [37, 79, 80] contain an almost isotropic initial spectrum, there- fore more general. The use of the ring distribution as the initial condition in [32] also introduced additional complications as the initial fluctuation could contain other electro- static/electromagnetic waves modes in addition to the whistler modes. More importantly, both set of studies were based on 2D simulations and thus do not allow the fluctuation spectra to be gyrotropic as what would be expected in a homogeneous system. Hence, to be able to draw a definite conclusion on the preferential cascade direction and anisotropy, a fully 3D PIC simulation performed for identical setup and initial conditions is necessary. 16 Figure 1.7: 3D gyrokinetic simulation results of antenna driven kinetic Alfv´ en turbulence from [93]. Collisionless damping via the Landau resonance with the electrons is sufficient to account for the measured heating as a function of scale in the simulation, without the need for significant Ohmic dissipation. Measured heating of the electrons from the sim- ulation by perpendicular wavenumber, Q e (k ? ) (black solid), an estimate of the electron heating based on linear wave-particle dampingQ wp (k ? ) (blue dotted), and the Ohmic heat- ing rateQ (k ? ) (red dashed). Landau damping sufficient to account for electron heating in simulation. Recently, Wan et al. and Karimabadi et al. [97, 46] performed PIC simulations using a large-scale sheared flow configuration which is subject to the Kelvin-Helmholtz instabil- ity; the consequent long-wavelength turbulence undergoes a forward cascade developing small-scale current sheets (Fig. 1.6(a)) which efficiently dissipate the fluctuation energy via associated parallel electric fieldsE k . The results of [46] shows that at ion gyro scales, turbulence exhibits both Alfv´ en (A) modes and magnetosonic-whistler (M) waves (Fig. 1.6(b)), and the parallel electron heating is consistent with both Landau damping of waves andE k generated by reconnecting current sheets. Their analytic estimations also suggest that current sheet heating due to reconnection in the electron layers is at least 100 times larger than that due to KAW wave-particle heating. A very recent gyrokinetic simulation of kinetic Alfv´ en turbulence [93] driven by an oscillating antenna also shows current sheet formation at electron gyroradius e scales. 17 They concluded that collisionless damping of KAWs via the electron Landau resonance is sufficient to account for the measured electron heating in their simulations. In Figure 1.7, the black solid line is measured heating of the electrons from the simulation by per- pendicular wavenumber, the blue dotted line is an estimate of the electron heating based on linear wave-particle damping, and the red dashed line is the predicated electron heating by the collisional Ohmic resistivity. It is clear that Ohmic heating rate is several orders of magnitude less than wave-particle damping, and their argument for the negligible Ohmic heating is that the electron-ion collisions are insufficient to significantly heat the electrons in a weakly collisional plasmas. However, it is trivial that ohmic heating is weak in a colli- sionless plasma. The more important part of this figure is the agreement between the black and the blue lines. 1.2.2 Recent observations of short-wavelength solar wind turbulence Recently, there have been an outburst of observation studies of short-wavelength solar wind turbulence [4, 77, 76, 69, 62, 64, 66, 65, 3, 11, 12, 22, 23, 41, 68, 72, 81, 88, 47]. Here we only describe several important studies that are used for direct comparison with our simulation results throughout this dissertation. Sahraoui et al. [76] analyzed recent data from the Cluster mission spacecrafts’ magne- tometers, and showed that the magnetic spectra at wavelengths shorter than the proton-scale spectral break in the frequency range 0.5 Hz < f < 50 Hz scale withjB(f)j 2 / (f) 2:8 , and suggested that there is a second break near the inverse electron gyroradius e , where e = p 2v te =j e j andv te is the electron thermal speed, with still more steeply decreasing spectra scaled as/ (f) 3:54:0 at approximately 50 Hz<f, as illustrated in the left panel of Fig. 1.1. We denote this fitting model for two power law dependance around the elec- tron scale as the “two-power-law break model”. Their results were further interpreted as consisting of KAWs, since all the data with ! r p are consistent with kinetic Alfv´ en 18 wave dispersion relations. Narita et al. [62] from another study using four Cluster space- crafts indicate many observed magnetic fluctuations with! r > p atk p > 1, indicative of magnetosonic-whistler fluctuations. These differences are compared in Fig. 1.2(b). Alexandrova et al. [4] from another Cluster spacecraft data study also showed a power law dependancef 2:8 for magnetic fluctuations at frequency larger than p . However, their results are more in favor of a exponential rollover at the electron scale instead of a distinct second spectral break to further steeper power law near inverse e or e as suggested by [76]. In their very recent work [3], a statistical analysis of 100 dissipation range measure- ments were performed to determine the form of the magnetic energy spectrum in the solar wind at the high frequencies corresponding to the scale of the electron gyroradius. They applied the Taylor hypothesis to get the wavenumber from the frequency,k ? = 2f=V sw , whereV sw is the solar wind velocity. Whistler waves intervals were discarded since the Taylor hypothesis is not applicable for v > V sw , where v is the phase velocity. The analysis suggested an empirical universal scaling form of the dissipation range given by / k ? exp(k ? e ), where the spectral index =2:63 0:15 with a typical value taken to be'8=3, as illustrated in Fig. 1.1(b). We denote this fitting model for exponential rollover in the short-wavelength dissipation range as the “exponential model”. Chen et al. [23] reported the first measurements of electron density fluctuation power spectra in solar wind turbulence at spacecraft frequencies from 0.1 Hz to approximately 40 Hz were using Artemis spacecraft data. The spectra scale asf 2:7 in the short-wavelength regime, which is consistent with previous observations of magnetic field spectra, as shown in Fig. 1.8(a). The results of [38] suggested that, the differences between the electron density fluctuation properties of KAWs and whistlers, such as short-wavelength electron density structures, may provide insight to help resolve the controversy concerning electron- scale solar wind turbulence. Recently, Kiyani et al. [47] investigated the magnetic field measurements from the Cluster and ACE spacecrafts and suggested that the probability density functions (PDFs) 19 (a) Chen et al. 2012 (b) Kiyani et al. 2013 Figure 1.8: (a) Power spectra of electron density fluctuations (PSD) from calibrated data with artificial spikes removed and smoothed [23]. (b) Normalized transverse and paral- lel probability density functions of fluctuations in the inertial and dissipation ranges with Poisonnian error bars attached [47]. observed in the solar wind at 1 AU show not only enhanced tails at large amplitudes, but also more strongly peaked than Gaussians at small amplitudes. This indicates that the observations are summed over many different small-scale measurements in inhomogeneous plasmas, so that the PDFs are summed over several different local Gaussians which result in a more peaked distribution at small amplitudes. They also found that the PDFs ofB ? andB k are similar for KAWs, but no observations of this anisotropy on the whistlers have been made, as illustrated in Fig. 1.8(b). 1.3 Motivation and objectives In summary, as described in the previous sections, there are three questions of fundamental importance in the current short-wavelength solar wind turbulence research. 20 What are the relative contributions to fluctuating energy dissipation from the quasi- linear scenario and the intermittency scenario in short-wavelength turbulence? Within the quasilinear scenario, what are the relative contributions of magnetosonic- whistler waves and kinetic Alfv´ en waves to short-wavelength turbulence? What are the 3D characteristics of whistler turbulence evolution and their effects on energy distribution and dissipation in the homogeneous solar wind plasmas? The research described in this dissertation is focused on understanding one of the two char- acteristic turbulence mentioned above, the whistler turbulence, in a collisionless, homoge- neous, magnetized plasma. The goal of this dissertation research is to advance the under- standing of the physics of whistler turbulence and its role in energy transport and dissipation in the solar wind. As described in section 1.1.1, in the high frequency short-wavelength regime (e.g., whistler mode), the MHD fluid description of turbulence fails. EMHD models, in which the ions are assumed to be stationary and electrons are represented as a fluid, are based on a low-frequency, long-wavelength approximation. Because of whistlers’ high frequency and short wavelength nature, EMHD models cannot adequately describe the plasma physics of this short wavelength regime. A kinetic approach is needed to fully reproduce the physics of dispersion and dissipation. PIC simulations, which represent both electrons and ions as superparticles, capture not only the wave-wave interactions represented by fluid descrip- tions such as EMHD, but also wave-particle interactions such as Landau and cyclotron damping which could be important dissipation mechanisms of cascading whistlers. The limitations of the PIC simulation method are not due to any analytic model approxima- tions, but rather are due to the constrained computing resources. PIC simulations in fully 3D system of whistler turbulence evolution are necessary. Pre- vious 2D simulations are intrinsically not gyrotropic, in contrast to what would be expected 21 in a homogeneous physical system. Such shortcoming may lead to different whistler cas- cade process from realistic situation in 3D systems. 3D kinetic PIC simulation of whistler turbulence evolution is an extremely challenging computational problem. The computa- tional requirement for such simulations is just within the reach of the capabilities of the most advanced massively parallel supercomputers which, for the first time, allow such a study to be performed. We extend the simulation approach in [79] to a fully 3D config- uration (Fig. 1.4(b)) and study whistler turbulence evolution evolving from an ensemble of preloaded whistler modes under different plasma parameters, such as initial fluctuation amplitude, electron beta e . The 2D simulations of [37, 79, 80] were carried out in a rela- tively small simulation box of dimensions 100 e 100 e ; the initial fluctuations imposed upon the system were of the order of those dimensions, so that only a forward cascade to shorter wavelengths could proceed. In order to determine the final state of whistler turbu- lence, it is necessary to consider larger simulation boxes to allow both a forward cascade to shorter wavelengths and an inverse cascade to longer wavelengths. We could then inves- tigate the competition of the forward vs. inverse cascade and wavevector anisotropy. We carried out the first fully 3D full-particle electromagnetic PIC simulations of whistler turbulence to accomplish the following specific research objectives: What are the 3D characteristics of whistler turbulence evolution and their effects on energy forward cascade, electron velocity distributions and dissipation in solar wind plasmas compared with those in 2D systems? How do 3D whistler turbulence characteristics, such as the energy cascade rates, the electron heating rates, the wavevector anisotropies, the fluctuation spectrum depen- dences, and the different mode couplings, vary under variable solar wind parameters, i.e., fluctuation amplitude and electron beta? 22 Is there any effective way to distinguish whistler modes and KAWs in short- wavelength solar wind turbulence and therefore guide future spacecraft measure- ments and observations? What are the relative contributions to fluctuating energy dissipation from the quasi- linear scenario and the intermittency scenario in whistler turbulence? What is the preferential evolution direction, energy cascade rate and wavevector anisotropy of whistler turbulence in a system large enough to permit both forward and inverse cascades? How does the forward cascade process compete against the inverse cascade process, and compare against previous forward cascade process? 1.4 Dissertation outline The contents of the dissertation are organized as follows. Chapter 2 introduces the computational approach used in this work, i.e., electromag- netic particle-in-cell modeling. Special attentions are give to various optimization and parallelization methods and numerical techniques implemented in this disserta- tion. Chapter 3 presents the first fully 3D PIC simulation results of whistler turbulence forward cascade in the solar wind, and compare the different characteristics between 3D and 2D systems. Chapter 4 analyzes various plasma turbulence properties such as wavevector anisotropies, spectra dependences, and electron heating rates under variable fluc- tuation amplitude solar wind condition. 23 Chapter 5 presents both linear theory results and PIC simulation results of whistler dispersion relations, various turbulence properties and energy cascade rates as func- tion of propagation angle and variable electron beta e solar wind condition. Besides, it introduces frequency mismatch concept based on linear theory results to test whether the nonlinear wave-wave coupling of whistler turbulence forward cascade may lead to wavevector anisotropies. Suggestions on future observational work to resolve the wave identity in short-wavelength turbulence is also given. Chapter 6 investigates the linear wave-particle damping and nonlinear dissipation of fluctuation energy of whistler turbulence under variable solar wind condition. Chapter 7 address the issue of the competition of forward vs. inverse cascade and its effects on previous results. Chapter 8 concludes this dissertation research and discusses future work. During the course of this dissertation research, the results of Chapter 3-5 have already been published in [19, 35, 20]. The results of Chapter 6 and 7 are also currently being prepared for submission. 24 CHAPTER 2: COMPUTATIONAL APPROACH In this chapter, we give an overview of the computational methods used in this research work, and present the newly developed computational and numerical features. In section 2.1, a linear dispersion solver is introduced, followed by a discussion of electromagnetic particle-in-cell modeling and its parallelization and optimization in section 2.2. In sec- tion 2.3, numerical techniques developed and used for the whistler turbulence simulations are presented. Finally, computational challenges and supercomputing facilities are briefly described in section 2.4. 2.1 Linear Vlasov dispersion solver The most accurate description of a plasma is based on the kinetic theory. Kinetic theory must be used when the fluid approximation is not sufficient to describe the dynamics of the plasma and the velocity distribution must be utilized. The Vlasov equation 2.1 forms the basis of collisionless plasma kinetic theory. It describes the evolution of the distribution function in phase spacef(x;v;t), wheref(x;v;t) is the average number of particles per unit six-dimensional phase space volume df dt = @f @t +vrf + q m (E +vB) @f @v = 0 (2.1) Our numerical dispersion solver [33] solves a linearized Vlasov dispersion equation for electromagnetic fluctuations at an arbitrary direction of propagation relative to the back- ground magnetic fieldB o in a homogeneous, collisionless plasma. Each plasma compo- nent is assumed to have a bi-Maxwellian velocity distribution with no relative flow speeds; 25 charge neutrality and zero current are also assumed. Specifically, the code solves Eq. 2.2 for various plasma fluctuations of the system, with the elements of tensorD(k;!) as given by Eq. 2.3 wherek = ^ yk y +^ zk z , and the conductivitiesS(k;!) as described in section 7.1.1 of [33]. detjD(k;!)j = 0 (2.2) D xx =! 2 k 2 c 2 +k 2 c 2 X j S xxj D xy =k 2 c 2 X j S xyj D xz =k 2 c 2 X j S xzj D yx =k 2 c 2 X j S yxj D yy =! 2 k 2 z c 2 +k 2 c 2 X j S yyj D yz =k y k z c 2 +k 2 c 2 X j S yzj D zx =k 2 c 2 X j S zxj D zy =k y k z c 2 +k 2 c 2 X j S yzj D zz =! 2 k 2 y c 2 +k 2 c 2 X j S zzj (2.3) Equation 2.2 is a transcendental equation for the complex variable ! given the real wavevectork; for each choice of plasma parameters, Eq. 2.2 has many solutions! =!(k). Although, with various approximations, one may obtain analytic expressions for the least damped or the fastest growing of these solutions, we only use the exact numerical solutions of Eq. 2.2 given by the dispersion solver in this dissertation. 26 2.2 Electromagnetic particle-in-cell modeling A major component of this research is based on the particle-in-cell technique, a tool in studying kinetic problems in plasma physics. A particle-in-cell simulation attempts to emulate plasma behavior by modeling plasma as millions or billions of charged test “super- particles” and following the evolution of the orbits of individual test particles which are interacting with each other in a discretized field [7]. The evolution of the plasma is deter- mined by the laws of interactions between the superparticles using the discretized version of Maxwell’s equations, Newton’s second law of motion, and the Lorentz’ force equa- tion, etc. The superparticles can be located anywhere within the simulation domain but the macroscopic field quantities are defined only on the discrete grid points, therefore, two interpolation (gather/scatter) steps are needed to link the particle orbits and the field com- ponents. While the particle model follows the motion of the plasma on the finest space and time scale, the scope of the physics that can be resolved in a simulation critically depends on the computational power. In this section, we first introduce the basic electromagnetic particle- in-cell (EMPIC) model, then describe the parallelization of our 3D-EMPIC code, followed by optimization techniques developed in our code which improve the runtime performance. 2.2.1 Electromagnetic finite-difference time-domain method An EMPIC code, as used in this specific research, simulates plasmas using only the fun- damental physics laws, i.e. Maxwell’s equations for the macroscopic field and Newton’s second law for individual particle trajectories [100] rE = 4e(n i n e ) (2.4) rB = 0 (2.5) 27 y B X Y Z z B x B y y J E , z z J E , x x J E , Figure 2.1: The Yee lattice shows the location of the field components on the staggered, finite-difference mesh, where grid points are at the corners of the cell [107]. @E @t =crBJ (2.6) @B @t =crE (2.7) d mv dt =F =q E +v B c ; dx dt =v (2.8) where relativistic effects are included such that = 1= p 1jvj 2 =c 2 . Eqs 2.4 through 2.8 are the equations to be solved in our 3D-EMPIC code. In our code, currents are deposited using a rigorous charge conservation scheme [96]. Note that this scheme updates the field purely from the local data, thus in general runs efficiently in parallel. This scheme requires the use of a fully staggered grid mesh system in which E and Jdt are defined at midpoints of cell-edges while the B components are defined at the midpoints of the cell-surfaces. The staggered grid mesh system, known in the computational electromagnetic community as the Yee lattice, is shown in Fig. 2.1 [107]. 28 The self-consistent electromagnetic field is then solved using a local finite-difference time- domain (FDTD) method to the full Maxwell’s equations: E n+1 E n =dt(crB n+1=2 )dtJ n+1=2 (2.9) B n+1=2 B n1=2 =dt(crE n ) (2.10) where the superscriptsn + 1=2 andn + 1 represent the time level. The trajectory of each particle is integrated using a standard time-centering leapfrog scheme for Eq. 2.8: u n+1=2 u n1=2 = q m E n + u n+1=2 +u n1=2 2 n B n c dt (2.11) x n+1 x n = u n+1=2 n+1=2 dt (2.12) whereu = v, E and B are interpolated from the grids to the particle positions. The basic algorithm and normalization of our code can be found in appendix A.2. This code had been successfully applied in several previous relevant simulations such as [39]. 2.2.2 Parallelism on distributed-memory architecture The parallelization implementation of this EMPIC code on distributed-memory parallel computers is described in details in [100, 101, 99]. It is shown that the code runs with a high parallel efficiency 95% for simulations using over 10 8 particles and 512 processors. The spatial domain can be partitioned into 1-, 2-, or 3-dimensional domain decomposition, and each processor is assigned a subdomain and all the particles and grid points in it. When a particle moves from one subdomain to another, it must be passed to the appropriate pro- cessors, which requires inter-processor communication. To ensure that the gather/scatter operations can be performed locally with no inter-processor communication, each proces- sor stores “guard cells”, i.e., neighboring grid points surrounding a processor’s subdomain 29 Particle Move Current Deposit Guard Cell Summation E advance Guard Cell Exchange for E First half of B advance Start loop Second half of B advance Guard Cell Exchange for B Partition Manager if needed Partition Manager if needed “Particle Push” Stage “Field Solve” Stage Guard Cell Exchange for B Particle Trade Figure 2.2: Iteration flow chart for the parallel 3D-EMPIC. Routines in the rounded boxes require no inter-processor communication, whereas routines in the rectangular boxes involve communication [100]. which belong to another processor’s subdomain, shown as dark grey area in Fig. 2.3 (a). Inter-processor communication is thus necessary to exchange guard cell information. Figure 2.2 shows the iteration flow chart of our parallel 3D EMPIC code. Subroutines in the rounded boxes are non-communicative computation operations. Subroutines in the rectangular boxes are communication operations. On a parallel computer, each processor executes computation routines independently using its own data arrays. The computations are linked together through global Message Passing Interface (MPI) communications. The global periodic boundary conditions are also imposed in the communication subroutines to avoid additional loops over grid points and particles. 30 Table 2.1: 3D-EMPIC parallel efficiency for some benchmark cases Weak scaling case mesh size up to 2048 3 time spent per step (sec) efficiency 1024 MPI processes 3.51 0.96 16384 MPI processes 7.05 0.48 Strong scaling case mesh size 512 3 time spent per step (sec) efficiency 256 MPI processes 3.46 0.98 256 MPI processes 8 OpenMP threads 6.24 0.54 2.2.3 Parallelism on hybrid shared- and distributed-memory architec- ture Recently, our code [32, 19, 35, 20] was modified, upgraded, and undergone extensive opti- mizations. Shared-memory parallel computing has been implemented besides the origi- nal distributed-memory parallelism to improve the multithreading parallel performance on cluster equipped with multi/many-core nodes. The new code, follows the idea in [61], implements a hybrid MPI + OpenMP parallelism for domain decomposition, such that OpenMP directives are used on shared-memory multicore nodes, and MPI are used for communication between these physical distributed nodes. All non-communicative subrou- tines including those most computation-intensive operations, such as particles move and current deposit are parallelized using multithreads to fully take advantage all available pro- cessor cores on cluster’s multi-core nodes. Inter-node communication latency increases as the total number of communication pro- cessors increases. Jobs using very large number of communication processors could easily fail due to intermittent network errors. For example, our largest simulation use 2048 3 grid cells with 16384 MPI processes on the Yellowstone supercomputer at the National Center for Atmospheric Research Supercomputing Center (NCAR) [104]. We have run for two times, only one time successfully finished, the other fails due to network errors. The results 31 0 P 1 P 2 P 8 P 7 P 6 P 5 P 4 P 3 P Guard cells layer MPI process (Primary domain decomposition) pseudo-linked-list cells ordered by space-filling curve Y X (a) Domain decomposition 0 1 2 3 4 5 6 7 8 9 10 2 4 6 8 10 12 14 Log 2 (total time) (seconds) Log 2 (MPI processes + OpenMP threads) Ideal 100% parallel efficiency MPI only Hybrid MPI+OpenMP (b) Figure 2.3: (a) The physical volume is spatially decomposed into processes, P n . Each process consists of a number of blocks of computational cells in a pseudo-linked-list order (see section 2.2.4 for definition). Cells inside each block are traversed concurrently by threads (denoted by red dots with arrows) to compute the particle motion and solve the electromagnetic field. (b) Comparison of scalability of the hybrid parallel implementation and the purely distributed parallel implementation. All cases are performed on dual quad- core AMD Opteron 2:3 GHz processor nodes. For hybrid cases, 8 threads are spawned mapping 8 physical cores of a cluster node. indicate that the communication cost between these 16384 MPI processes are significant increased, especially for those collective commutations. Our current code runs at a parallel efficiency of 96% for 1024 MPI processes, whereas the efficiency is reduced to 48% for 16384 processors for this weak scaling case, shown in Table 2.1. The hybrid parallel implementation is aimed to reduce these communication cost and random network failure, via reducing the total number of MPI processes, and utilizing these saved processor cores, which were originally independent MPI partitions, for mul- tithreading parallelization. The overall parallel efficiency could be higher than the purely distributed case if a large number of communication processors are involved. Besides, such a hybrid parallelism could be easily transplanted to various heterogeneous platforms and configured to explore each different system’s full computing power. Figure 2.3(a) shows the hybrid MPI + OpenMP parallel implementation of our current EMPIC code for a 2D domain decomposition. In this hierarchical decomposition, each 32 MPI process is assigned a spatial-decomposition subsystem, and OpenMP threads of each process handle blocks of “pseudo-linked-list” cells inside the subdomain in a predefined order, such as Morton Z curve order shown here. Figure 2.3(b) show the scalability of the hybrid parallel implementation for a 2D homogeneous plasma simulation of 1024 2 grid cells with 64 particles per cell per species. In hybrid cases, we spawned 8 threads per process mapping 8 physical cores of a node. The hybrid version’s parallel efficiency is roughly 50% of that of the purely distributed version. However, in Fig. 2.3(b), with half of the total amount of MPI processes, the hybrid version could reach same same speedup. With the number of MPI processes continues to increase, the distributed version suffers a decrease in parallel efficiency due to communication cost, as indicated by Table 2.1. In contrast, the hybrid version use at least a factor of 8 less processes, and thus have much less communication latency than the distributed version, and may have a better the overall runtime performance. We note that a GPU version of this 3D-EMPIC code have been implemented using Compute Unified Device Architecture (CUDA) programming model [18] which exploits even more advantage of multithreading shared-memory parallelism on graphics processing units. Nowadays, processors have more and more multi-cores on a single chip, such hybrid parallel implementation is essentially important to accommodate future high performance computing supercomputers. 2.2.4 Adaptive particle sorting An important approach to optimize a EMPIC code is to improve data and computa- tion locality. If one could always fully utilize processors’ cache memory (faster speed and smaller size) and reduce cache miss, the simulation speedup could be dramatically increased, especially for large-scale particle simulation like our case, of which the length of particle array on each node is on the order of 100 10 6 . Specifically, the key idea is 33 to preserve nearest-neighbor cells spatial proximity in memory, so cache reuse is maxi- mized and reloading from main memory to cache is minimized. This requires a dedicated order for both accessing the block of cells and storing particles in the particle array, since the most time-consuming parts, particle move (gather) and current deposit (scatter), require both particle and field grid information. Such order is achieved through space-filling curves [60], which mapping from the multi-dimensional space to one-dimensional list to preserve spatial proximity of consecutive list elements. Two specific space-filling curve, Morton Z Curve and Hilbert curve are been implemented, and we found no obvious performance difference, therefore, we choose Morton Z Curve order for all our 3D simulations. A particle sorting scheme using this spatial proximity preserving algorithm is devel- oped and implemented. With simulation going on, particles are randomized and could be anywhere in the domain. Two particles consecutively in the particle array (also continuous in memory) could be located far away in the simulation space domain, and thus, require dif- ferent portion of grid cells for gather/scatter steps and create huge amount of cache miss in results, degrading the computation speed. Thus, to explore computation locality, particles needs to be sorted in respect of their position for a certain period. The primitive way is to sort the particle array by particle cell index line by line, in the column-major order in For- tran. Lets denote this order as “scan-line”. In this work, Morton Z curve is used to create a one-dimensional cell list order for sorting the particle array in order to preserves spatial proximity and enhance computation locality. In gather/scatter steps, each thread/process loop over the particle array to compute force or current density, which is essentially loop- ing over the cells in the order defined by space-filling curves. We denote this as pseudo- linked-list cells, as it is opposed to the real linked-list cells which restrictedly follow the cell order and loop over the particles in each cell. We do not use the real inked-list cells in gather/scatter steps, because it requires updating the link list every step. This is very costly in our problem since our particle array is very large and irregular memory access pattern should be avoided if possible. Therefore, looping over the particle array is still the primary 34 order, and looping over particle array is not exactly following the space-filling curves’s cell order. The detailed sorting algorithm implements an out-of-place linear-timeO(n) counting sort [98], and only a non-full sort for particle’s cell index is necessary. Due to the overhead involved in sorting, we use an adaptive cost/benefit algorithm [74] to decide the time period between sorting and ensure sorting penalties are recovered. Sorting either too infrequently or frequently results in non-optimal execution time due to either increasing iteration time or excessive sorting overhead respectively. Sorting is only done when its performance impact outweighs the cost required to perform the sort, that is the following condition is satisfied: tSumnLastSorttSortedtSort (2.13) wheretSum is the total time elapsed since the last sort,nLastSort is the number of steps since the last sort, tSorted is the average time spent per step within three steps after the last sort, andtSort is the time spent for the last sort. This condition can be derived under the assumption that the time difference between the measured time per step and the time per step if all particles perfectly sorted ( tSorted) is increases linearly with increasing time step, such that y = mx, where m is the slope constant, x is the step, and y is the time difference for that step. We also assume tSort and tSorted are constant, so the optimization equation is: min( tSort + 0:5mx 2 +tSortedx x ) (2.14) The derivative of Eq. 2.14 gives x = q 2tSort m or tSort = 0:5mx 2 . In the code, this is satisfied whenever condition 2.13 is true. 35 0 100 200 300 400 500 600 700 800 900 1000 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 iteration step time per step (seconds) no particle sorting scan −line sorting Morton order sorting (a) overall view 0 10 20 30 40 50 60 70 80 90 100 1.65 1.7 1.75 1.8 1.85 1.9 1.95 2 2.05 2.1 2.15 2.15 iteration step time per step (seconds) (b) zoomed-in view Figure 2.4: Comparison of simulation time per step for cases with and without the sorting mechanism, and with different space-filling curve sorting orders. Figure 2.4 shows the comparison of the simulation time per step for a typical 512 3 homogenous plasma simulation with sorting in Morton curve orders (blue line), with sort- ing in scan-line order (green line), and without sorting (red line). In the runs with the sorting mechanism, the computing time remains a constant between each sorting period while in the run without the sorting mechanism, the computing time gradually increases at 36 Table 2.2: Parallel I/O efficiency on NCAR’s GPFS file system mesh size 512 3 , file size 3 GB number of I/O processes time spent in write (sec) efficiency 1 5.33 1.0 8 1.67 0.4 64 0.33 0.25 mesh size 2048 3 , file size 192 GB number of I/O processes time spent in write (sec) efficiency 1 341.33 1 512 12.43 0.05 each iteration step, as clearly shown in Fig. 2.4(b). Time per step for the no-sorting run reached about 1.75 times as much as the time per step for the sorting cases at iteration step 1000. The speedup for the total simulation time from particle sorting cases is already about 1.5 after iteration step 1000 and will increase more for longer simulation time period, for example, our simulation cases usually run for 90000 steps. The use of Morton space-filling curve order over the scan-line order for particle sorting delivers about an extra 5% more speedup for these 3D simulations. 2.2.5 Parallel I/O and fast Fourier transform Parallel I/O operations are implemented in our current 3D-EMPIC code to facilitate very large size data output and loading. We use MPI-IO subroutines provided by MPI-2 library to achieve parallel I/O. However, we note that these parallel I/O operations usually need dedicated parallel file systems to support hundreds if not thousands processes to concur- rently write to or read from file systems, and specific tuning is also desired. Our experience is that MPI parallel I/O operations on Lustre distributed file system usually perform faster and more accurately, whereas operations on Network File System (NFS) are slower and subject to error-prone when the number of writing processes are increased to hundreds. 37 Table 2.2 shows the parallel efficiency for some typical output file size in our simulation on NCAR’s General Parallel File System (GPFS). Although the parallel efficiency is not good, it does significantly reduce the time spent in writing simulation data to file systems. We note that the Yellowstone supercomputer and its associated GPFS file system are recently deployed, and there does have some random I/O failure or data lost happened. Parallel fast Fourier transform (FFT) is also implemented in the current code by calling the cluster FFT functions from Intel Math Kernel Library. These functions are designed to perform fast Fourier transforms on a cluster environment using MPI to communicate between different processes. We usually need three dimensional Fourier transform, how- ever, the cluster FFT functions only provide one dimensional parallelization, and the largest partition could not exceed the length of the last dimensions in Fortran, that is a 512 3 cell size simulation could be at most partitioned in 512 MPI processes for FFT. The time spend for FFT subroutines is about 5 7 times as much as the time spent for particle sorting. Luckily, we only need to call it periodically, and only affect the overall runtime perfor- mance slightly. 2.3 Numerical techniques To facilitate the whistler turbulence simulation and data diagnostics, we have utilized sev- eral numerical techniques for various purposes. Our simulations run on small perturbed fluctuations, therefore, a lower numerical noise is very important to separate the physics and numerical artifacts. Energy conserving interpolation scheme [13] for particle push stage have been implemented to reduce artificial numerical heating. A local and inexpen- sive method, Marder passes for enforcing Gauss’ law [53], has also been used periodically to clean arithmetic error induced violations of Gauss’ law. We also apply the current den- sity smoothing techniques originated from TRISTAN code [14] to reduce the unnatural noise created from finite number of superparticles. Any source of Maxwell’s equation, i.e., 38 current density in our code, has been smoothed every time step before it is used to solve electric fields in Eq. 2.6. In this section, we describe three other numerical techniques that have been used to ensure the success of our turbulence simulation. 2.3.1 Turbulent spectrum initialization To create an initially broadband incoherent fluctuation spectrum close to a turbulence, we load an ensemble of right-hand polarized whistler waves of same magnetic fluctuation amplitude but random phase at t = 0. Frequency and wavevector relationships of these waves are solutions from the linear dispersion solver in section 2.1. Each wave’s wave- length (and thus wavevector) is specifically chosen that it could be resolved accurately on the simulation’s computational grid. Although our dispersion solver could calculate the fluctuation amplitude between elec- tromagnetic field component of each wave, it does not give information about phase relation between components, that is the polarization property. Therefore, we follow the approxi- mate relationships derived from [95] and parallel propagating whistler mode as below, By = 0 Bx =i Bz = 0, if k ? k k < 0, Bz = (2.15) 39 fork z > 0 : Ex = 0 Ey =i Ez =i, ifa 2 < 0, Ez =i fork z < 0 : Ex = Ey =i Ez =i, ifa 2 < 0, Ez =i (2.16) where a 2 = kykz 2 e 1+k 2 y 2 e . We note that Eqs. 2.15 and 2.16 are based on the assumption that the wave propagates in the k y k z plane, where k z is parallel direction to B o and k x is the positive out of the plane direction. Therefore, for a 3D initial spectrum, we rotate the k y k z spectrum 90 o aboutk z axis tok x k z withk y being the positive out of the plane direction. We also rotate it 45 o and 135 o in the same fashion. The only difference is that we need to choose different wavelength modes so that the computational grid could resolve accurately. Current densityJ can be derived from Amp` ere’s law Eq. 2.6, and is imposed as a pertur- bation of the average velocity of the electron velocity distributionv e , such thatJ =qn e v e assuming protons are at rest due to their large inertia in this relatively high frequency regime. We neglect the electrostatic and particle density contributions to the initial fluc- tuations, because the calculations should adjust the initial fluctuations to self-consistent fluctuations within a few electron plasma periods [80]. 40 Figure 2.5: Cold plasma dispersion curves,!(k)=! p = cos 1=2 kx, showing compensa- tion and smoothing fora 1 = 0:5 and variousa 2 values. The dashed line is the cold plasma dispersion for no compensation and no attenuation [7]. 2.3.2 Spectral filtering To eliminate numerical noises and unreliable high mode number structures, a spectral filter in Fourier domain is applied to current densityJ for some cases. We use the smoothing functionSM(k) described in [7], SM(k) = exp(a 1 sin 2 kx 2 a 2 tan 4 kx 2 ) (2.17) where k is the mode number, x = 2 Lx , and a 1 and a 2 are two adjustable parameters. a 1 sin 2 kx 2 is used to cancel forkx. 1 theO(k 2 ) error in the dispersion relation due to spectral transform, and is called compensating or boosting. Seta 1 = 0:5 makes dispersion error to order (kx) 4 . a 2 tan 4 kx 2 is to attenuateJ(k) at short wavelengths, generally called smoothing. A comparison of this spectral filter effect for cold plasma dispersion 41 curves !(k)=! p = cos 1=2 kx is shown in Fig. 2.5. Similar spectral filters have been successfully applied in hybrid simulations and other PIC simulations [17, 97]. 2.3.3 Wavelet transform To analyze simulation results as functions of both frequency and time, we apply the con- tinuous wavelet transform (CWT), a technique widely used in observations [47, 52], on the magnetic field fluctuations. Like the Fourier transform, the continuous wavelet trans- form uses inner products to measure the similarity between a signalf(t) and an analyzing function. In the Fourier transform, the analyzing functions are complex exponentialse i!t , and the resulting transform is a function of a signal variable!. In the CWT, the analyzing function is a wavelet, . The CWT compares the signal f(t) to shifted and compressed or stretched versions of a wavelet. Stretching or compressing a function is collectively referred to as scaling and corresponds to the physical notion of scale. By comparing the signal to the wavelet at various scales and positions, and if the signal is real value, the CWT coefficient is a real-valued function of scale and position. For a scale parametera> 0, and positionb, the CWT is: C(a;b;f(t); (t)) = Z 1 1 f(t) 1 p a ( tb a )dt (2.18) where C is the CWT coefficient and denotes the complex conjugate. Note that not only do the values of scale and position affect the CWT coefficients, the choice of wavelet also affects the values of the coefficient. The coefficient C represents how closely correlated the wavelet is with this section of the signal. The larger the number C is in absolute value, the more the similarity. In general, the signal energy does not equal one and the CWT coef- ficients are not directly interpretable as correlation coefficients. In our case the signal is magnetic field as function of timeB(t), the position thus represents the timet. While there is a general relationship between scale and frequency, no precise relationship exists, and 42 it is better to talk about the pseudo-frequency corresponding to a scale. Using a function scal2frq provided by MATLAB [105], we could obtain approximate scale-frequency rela- tionships for specified wavelets and scales. Morlet wavelet is chosen as the mother wavelet for comparison, since it is composed of a complex exponential multiplied by a Gaussian window, and thus is close to small perturbed magnetic fluctuations. 2.4 Computational requirements and facilities 2.4.1 Computational challenge The simulation parameters required by a fully 3D simulation of whistler turbulence present an extremely challenging problem computationally. First, in order to resolve the temporal turbulence cascade process, full particle EMPIC simulations must be performed for long periods of time (thousands of electron cyclotron periods). With time step typically atdt 0:05! pe due to numerical stability constraints, a simulation would require more than tens of thousands simulation time steps. Second, and more importantly, the simulation domain must be sufficiently large to contain several wave lengths for the longest wave length modes in each direction while the mesh size must be sufficiently small to resolve the shortest wave length modes from the cascade process. The mesh resolution also needs to satisfy De , where De is the Debye length, to avoid artificial numerical heating. Finally, a sufficiently large number of particles per cell is required to minimize the numerical noise. The simulations of whistler turbulence forward cascade presented in [19, 35, 20] used a simulation domain ofLx =Ly =Lz = 51:2 e or 512 De . With a cell resolution of = De = 0:1 e , the simulation domain contains 512 3 cells. The number of superparticles used is typically 17:2 10 9 . The total memory requirement is 2048 GB. With a time step at dt = 0:05! pe , the simulation would typically require 90000 simulation steps. 43 On the Yellowstone, using 1024 subdomains on 1024 processor cores, a typical whistler turbulence forward cascade simulation requires about 50 wall clock hours to complete. The size of the simulation domain of [37] only allows for forward cascade of the initial modes to shorter wavelength modes. To further study the competing effects of the forward and inverse cascade for this case, the simulation domain must be further enlarged to allow possible cascade to longer wave length modes. In this case, we increase the simulation size by a factor of 64 to include both forward and inverse cascade of whistler turbulence, that is, a domain of 2048 3 cells and 5:5 10 11 superparticles. The total memory requirement for such a run is 32768 GB 32 TB. Using a domain decomposition of 16384 subdo- mains, this simulation currently runs at a computing speed of about 6 hours wall clock for every 3000 simulation steps. At the time this dissertation is written, as far as the author’s knowledge, this simulation is the world’s largest kinetic simulation of solar wind turbu- lence. Simulations at similar scales have also been attempted in [46], which used the VPIC code [13] for 3D simulations of a domain of 992 992 1984 cells. 2.4.2 State-of-the-art supercomputing facilities The primary modeling and simulation studies are carried out extensively at the Pleiades supercomputer of the NASA Advanced Supercomputing (NAS) Facility. Recently, our code have been deployed at the Yellowstone supercomputer of the NCAR-Wyoming Super- computing Center. Early phase of development and test were conducted at USC High Per- formance Computing and Communications Center (HPCC). Both the Pleiades and the Yel- lowstone supercomputers are 1:5 petaflops clusters, and currently ranked within the top 20 of the world’s fastest supercomputers [106]. USC cluster is also ranked as the nation’s 5th fastest academic supercomputer. 44 CHAPTER 3: THREE-DIMENSIONAL PARTICLE-IN-CELL SIMULATIONS OF WHISTLER TURBULENCE FORWARD CASCADE: COMPARISON WITH TWO-DIMENSIONAL SIMULATIONS In this chapter, we present the results of the first fully three-dimensional particle-in-cell simulations of whistler turbulence in a magnetized, homogeneous, collisionless plasma. With an introduction in section 3.1, initial plasmas and computational parameters and solar wind fluctuations setup are given in section 3.2. The detailed simulation results are dis- cussed and the comparison with previous two-dimensional simulations are presented in section 3.3. Section 3.4 summarizes the new findings and concludes this chapter. 3.1 Introduction As introduced in chapter 1, there is currently substantial debate about the properties of the high frequency short wavelength turbulence that is far above the first proton-scale spec- tral break and extends down to and through the second electron-scale spectral break [88]. While new 3D gyrokinetic simulations of kinetic Alfv´ en turbulence have been made and the results show some characteristic properties observed in the solar wind, kinetic simulations of whistler turbulence are relatively few [37, 79, 80] and all in 2D, despite substantial fluid model simulations have been made [10, 25, 26, 83]. These 2D PIC whistler turbulence simulations are intrinsically not gyrotropic, in contrast to what would be expected in a 45 homogeneous physical system. A fully 3D kinetic model is apparently needed to throughly understand whistler fluctuations and to be able to compare results with observational data in the solar wind turbulence. Here we describe the first PIC simulations to examine the forward cascade of whistler turbulence in a fully 3D plasma model using particle-in-cell code 3D-EMPIC described in chapter 2. We choose our initial spectrum of whistler waves (see the next section 3.2) to have the longest wavelengths that fit into the simulation box; therefore by construction the cascade processes here are necessarily in the forward direc- tion. 3.2 Initial conditions and simulation setup We consider hereafter an electron-proton plasma where subscripte denotes electrons and p denotes protons. We denote the jth species plasma frequency as ! j q 4n j e 2 j =m j , the jth species cyclotron frequency as j e j B o =m j c, the jth species thermal speed as v tj p k B T kj =m j , initial electron thermal speed as v teo , and j 8n j k B T kj =B 2 o . We define , the angle of mode propagation, by k B o = kB o cos(). The dimen- sionless parameter characterizing the initial magnetic fluctuation energy is denoted by o k jBj 2 t=0 =B 2 o . Here “three-dimensional” means that the simulation includes variations in both three spatial dimensions, as well as three velocity space response of each ion and electron superparticle. The plasma is homogeneous with periodic boundary conditions in all three directions. The uniform background magnetic field isB o = ^ zB o so that the sub- scriptsz andk represent the same direction. Thus wavevectork = ^ xk x + ^ yk y + ^ zk k and k ? = p k 2 x +k 2 y . The total fluctuating magnetic energy density is obtained by summing over all simulation wavevectors. The 2D reduced magnetic fluctuation energy spectra by summation over the third Cartesian wavevector component, e.g., jBj 2 =8jB(k y ;k z )j 2 kx jB(k)j 2 (3.1) 46 In contrast, the 1Dk ? reduced magnetic fluctuation energy spectrum, which we here denote asE(k ? ), must be summed over bothk k and, the azimuthal angle of the perpendicular wavevector. Thus, in the continuum limit, jBj 2 =8 Z d 3 kjB(k)j 2 =8 = Z dk ? k ? E(k ? ) (3.2) In a 2D k k k ? system, the above equation is reduced to simply as the sum of the 2D spectrum over the Cartesian parallel wavenumber, jBj 2 =8jB(k ? )j 2 k k jB(k)j 2 (3.3) The average fluctuation wavevector anisotropy angle B is defined through the wavevector anisotropy factor [86]: tan 2 B k k 2 ? jB(k)j 2 k k 2 k jB(k)j 2 (3.4) An isotropic spectrum corresponds to tan 2 B = 1:0. Our concern here is the anisotropy of the forward cascaded fluctuations, so this quantity is evaluated by summations over wavevectors which satisfy 0:65jk e j 3:0. The system has a simulation domain size of 512 512 512 cells, with a grid spacing of = De = 0:1 e , it is equal to a spatial dimensions of L x = L y = L z = 51:2 e . For the 2D run shown later in this section, x direction is set to contain only 1 cell. The number of superparticle pairs per cell is 64, which gives roughly 17 billions of total super- particles in the 3D system. Simulation time step is chosen to bet = 0:05! 1 e to enforce Courant-Friedrichs-Lewy (CFL) condition. The total energy variation is less than 0:4% throughout the simulations. For these parameters, the fundamental mode in the simulations has wavenumberk e = 0:1227, and short wavelength fluctuations should be resolved up to the component wavenumber of k k e = k ? e = 4. The initial physical dimensionless parameters are e = 0:1, v te =c = 0:1, ! 2 e = 2 e = 5:0, m p =m e = 1836, T e =T p = 1:0 and 47 Table 3.1: Plasma parameters in the simulations Simulations Typical solar wind at 1AU e 0.1 1 v te =c 0.1 0.005 ! e = e 2.24 140 T e =T p 1 1 m p =m e 1836 1836 o 0.1 0:01 1:0 B o (nT ) 641 10 N e (cm 3 ) 20 20 T e (eV ) 5118 10 V te (km=s) 3:0 10 4 1500 V A (km=s) 3131 50 V sw (km=s) N/A 300 600 ! e (s 1 ) 2:5 10 5 2:5 10 5 e (s 1 ) 1:1 10 5 1800 e (km) 1.2 1.2 De (km) 0.12 0.006 o = 0:1. Plasma parameters used for the simulations are summarized in Table 3.1, and computational parameters are listed in Table 3.2. We note that variable e condition is addressed in chapter 5. Although due to explicit PIC algorithm constraints discussed in section 2.2.1, we have to choose an artificially increased electron thermal speed ratiov te =c to make simulations feasible to perform, the relativistic effects in our simulations are still negligible and similar to typical solar wind condition. Increasingv te =c results in a decrease of electron plasma-frequency-to-cyclotron-frequency ratio! e = e , however, linear disper- sion calculations shown in Fig. 5.3 imply that as long as the value! e = e & 1, properties of whistler modes remains unchanged. Three-dimensional spectrum of right-hand polarized whistler waves of same amplitude are imposed att = 0. Initial wavenumbers arek k=? e =0:1227,0:2454, and0:3682, andk ? e = 0. There is no solution under magnetosonic-whistler regime in purely perpen- dicular directionk k e = 0. The method of loading such an initial spectrum is described 48 Table 3.2: Computational parameters in the simulations Simulations = De ; = e 1.0, 0.1 L x;y;z = De ;L x;y;z = e 512, 51.2 ^ v te 1.0 t! e 0.05 superparticle pairs per cell 64 total superparticles 1:72 10 10 in section 2.3.1. The frequencies and amplitude of field components are derived from the linear dispersion equation for magnetosonic-whistler fluctuations in a collisionless plasma, but the subsequent evolution of the fields and particles are computed using the fully non- linear PIC simulation code. 3.3 Particle-in-cell simulations and discussions In this section, we present three simulation runs. Run 1 is a two-dimensional simulation with an initial spectrum similar to the corresponding runs of [37, 79]. It serves as both a comparison case to 3D runs and a validation case to previous 2D studies. Total 42 modes consisting of multiplying 7 parallel modes and 6 perpendicular modes are present in the k y -k z plane of the initial whistler spectrum. The initial spectra of Run 1 is of same pattern as the upper right panel of Fig. 3.1. Runs 2 and 3 are three-dimensional simulations, with the only difference in the initial spectra. Run 2 has same 42 modes in thek y -k z plane as those in Run1, and another extra 36 modes in the k x -k z plane, omitting double counting another 6 modes atk ? = 0. Run 3 is closest to the actual gyrotropic turbulence system. It has the same 78 modes as Run 2, but additional 72 modes corresponding to a rotation of the Run 2 modes through an angle of 45 o about thek z axis. Although we could make the initial spectra even more gyrotropic, simulation results (Fig. 3.4) imply that Run 3 spectra 49 Figure 3.1: Upper panels: Reduced magnetic fluctuation energy wavevector spectra at t = 0: (left panel) Run 2 in the plane perpendicular toB o , (center panel) Run 3 in the plane perpendicular toB o , and (right panel) Run 3 in thek y -k z plane containing B o . Lower panel: 3D magnetic fluctuation energy spectra of Run 3 att = 0. 50 are good enough to drive wave-wave interactions to keep wave cascade process going. Figure 3.1 illustrates the distribution of the initial modes for the two 3D cases. We note that although different runs have different number (N) of modes imposed in the system, the initial total fluctuating magnetic field energy density is still the same for all runs, i.e., o = P N n=1 jB n j 2 t=0 =B 2 o = 0:1. To validate the correctness of the initial whistler mode loading process, we plot the angle of perpendicular magnetic fieldB ? and electric fieldE ? as a function of time. Note that a single whistler wave is right hand circular traveling wave. When both positive and negative propagating whistler waves are superimposed, the ensemble become right hand circular standing wave as demonstrated in Fig. 3.2. Figure 3.3 shows 3D and 2D reduced magnetic energy spectra from the three simula- tions at a late-timej e jt = 447. As in the 2D Run 1, the imposition of an initial spectrum of relatively long-wavelength whistlers in the 3D runs leads to a forward cascade to shorter wavelengths and the development of a broadband, turbulent spectrum. The late-time spec- trum in thek x -k y plane perpendicular toB o is approximately gyrotropic, which is a con- sistency check on our homogeneous plasma simulation with initial conditions which are symmetric betweenx andy. The late-time reduced spectra in thek y -k z plane in all three runs, and the 3D spectra of Run 3 in lower left panel show that, for the first time, in a fully 3D system, the forward whistler turbulence cascade is anisotropic, preferentially transfer- ring energy to fluctuations with wavevectors quasi-perpendicular, rather than quasi-parallel, toB o . The most important new finding of 3D whistler simulations is the spectral break at the inverse electron inertial length scale. Figure 3.4 show thek ? reduced spectra from all of the three runs at a late-timej e jt = 447. The spectra from 2D Run 1 (red line) is proportional tok 4:4 ? over 0:25< jk ? e j< 3:0, similar to thek 4:5 ? late-time results from previous 2D PIC simulations [37, 79, 80]. No spectral break is reported from any of these previous works, whereas observations of solar wind turbulent spectra with breaks near the inverse 51 B ⊥ E ⊥ Run 1 Run 2 Run 3 Figure 3.2: Perpendicular magneticB ? and electricE ? field polarization as a function of time from all three runs. 52 Figure 3.3: Reduced magnetic fluctuation energy wavevector spectra atj e jt = 447 from Run 1 (upper left panel), Run 2 (center panels) and Run 3 (right panels). The upper panels represents spectra in the k y -k z plane, whereas the lower panels represents spectra in the plane perpendicular toB o . Lower left panel is the corresponding 3D magnetic fluctuation energy spectra of Run 3. 53 10 −1 10 0 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 E(k )/B 2 o |k e | Run 1 Run 2 Run 3 k −3.1 k −4.4 k −4.8 Figure 3.4: k ? reduced magnetic fluctuation energy spectra atj e jt = 447 from the 2D Run 1 (red line), and the 3D Run 2 (green line) and Run 3 (blue line) simulations. The blue dashed lines represent power law functions as labeled for comparison against the simulation results. electron inertial length were reported by [76, 4] and demonstrated from gyrokinetic simu- lations of kinetic Alfv´ en turbulence [45]. In contrast to 2D whistler simulations, 3D runs do show an obvious spectral break at electron inertial scale, with power law dependence at both side of the break, as illustrated in Fig. 3.4. The reduced spectra at relatively long wavelengths (0.36 < jk ? e j< 1.0) suggest ak 3:1 ? dependence, whereas at shorter wave- lengths (1.0 < jk ? e j< 4.0) the spectra become steeper with ak 4:8 ? dependence. The upturn in the spectra atjk ? e j> 4.0 corresponds to the noise level of the simulation. The measured power law spectral index =3:1, which is very close to reported spectral 54 indices such as =2:8 from recent observations [76]. The suggestion of two distinct power-law regimes of the turbulence is also similar to the prediction of [57] who used an EMHD model with isotropic fluctuations to derive ak 11=3 magnetic fluctuation spectrum at 1<jk e j. Figure 3.4 is the first time that the second spectral break at electron scale is reported for whistler turbulence, which reinforce the hypothesis that whistlers dominates and heat electrons at very short-wavelength solar wind turbulence. We note that the results for Run 2 and Run 3 are very similar, indicating that the late-time spectra are relatively independent of the detailed choice of initial fluctuations. Alexandrova et al. [4] reported observations of solar wind turbulence with an expo- nential decrease at scales shorter than the electron gyroradius e . These measurements are not necessarily in contradiction with our short-wavelength power-law result because our simulations have been run at e = 0:1 where Landau damping of whistler fluctuations is relatively weak, whereas the Alexandrova et al. [2009] measurements were almost all at larger e values where Landau damping is stronger and exponential decreases of spectra are more likely. We will address thek ? reduced spectra e dependance in chapter 5. Figure 3.5 compares the spectral wavevector anisotropies tan 2 B of the three runs a function of time. The results of the 2D Run 1 are very similar to that of Run I from [79], with tan 2 B 4 atj e jt' 400, again proving the correctness of our code. Both 3D runs exhibit stronger anisotropies, with Run 2 reaching tan 2 B ' 6 and Run 3 attaining tan 2 B > 10 at late times. Our interpretation of these results is that the much larger number of quasi-perpendicular modes available in the three-dimensional cases allows the nonlinear wave-wave interactions and perpendicular cascades to proceed much more efficiently than in two dimensions. The successively increasing number of modes with k ? components among these three runs implies more channels for perpendicular cascading and therefore a successively faster rate of such energy transfer for the fluctuations. Had we imposed more dense modes in wavevector space, the trend of these three runs indicates that such a 55 0 50 100 150 200 250 300 350 400 450 1 3 5 7 9 11 13 tan 2 B | e |t Run 1 Run 2 Run 3 Figure 3.5: Time histories of the spectral wavevector anisotropy factor tan 2 B from the 2D Run 1 (red line) and the 3D Run 2 (green line) and Run 3 (blue line) simulations. Evaluations here are taken on spectra over 0:65jk e j 3:0. condition should correspond to a highly anisotropic late-time condition with tan 2 B >> 1, and a spectral index value possibly closer to -2.8. Figure 3.6(a) shows the total fluctuating magnetic field energy density as a function of time. Results of all three runs are qualitatively similar, demonstrating a gradual decrease of the total fluctuating magnetic field energy density. However, the 3D simulations exhibit faster rates of fluctuation energy dissipation than the 2D Run 1. We provide a possible explanation for this difference at the end of this section. Figure 3.6(b) shows that electrons have greater parallel kinetic energy (solid lines) gains in the 3D runs than in the 2D Run 1. The perpendicular energy (dashed lines) gain of electrons is much weaker in all three runs, and there is no significant difference in this energy gain among the three cases. Figure 3.7 illustrates the reduced electron velocity distributions f e (v k ) and f e (v ? ) at t = 0, atj e jt = 447, and their differences. This demonstrates clearly that whistler turbu- lence preferentially heats electrons in directions parallel and anti-parallel toB o . Consistent 56 0 50 100 150 200 250 300 350 400 450 1 1.2 1.4 1.6 1.8 2 2.2 | e |t T e /T e,t=0 0 50 100 150 200 250 300 350 400 450 0.2 0.4 0.6 0.8 1 | B| 2 /| B| t=0 2 Run 1 Run 2 Run 3 (a) (b) Figure 3.6: Time histories of the (a) total fluctuating magnetic energy density and (b) the parallel (solid lines) and perpendicular (dashed lines) components of the total electron kinetic energy a function of time, for the 2D Run 1 (red line) and the 3D Run 2 (green line) and Run 3 (blue line), respectively. with Fig. 3.6(b), electrons are heated much more in parallel direction in the 3D runs than in corresponding 2D run, while there is minimum amount of difference in perpendicular heating. 57 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 v ! / √ 2v te o f e (v || ) Run 1 −4 −3 −2 −1 0 1 2 3 4 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 v ⊥ / √ 2v te o f e (v ) −4 −3 −2 −1 0 1 2 3 4 v ⊥ / √ 2v te o v ! / √ 2v te o Run 3 f e, e t = 0 f e, e t = 447 (f e, e t=447 − f e, e t=0 )/f(0) e, e t=0 Figure 3.7: Reduced electron velocity distributions as functions ofv k (top panels) andv ? (bottom panels) for Run 1 (left panels) and Run 3 (right panels). The grey lines represent initial velocity distributions at t = 0, values, the red and blue solid lines correspond to Run 1 and Run 3 atj e jt = 447, and the dashed lines represent the difference of the two distributions normalized by the initial value of that distribution. 58 Note that the electron Landau damping is associated with obliquely and quasi- perpendicular whistler modes and scatters electron parallel velocity as described in section 1.1.2, and 3D runs have much larger number of quasi-perpendicular modes available for facilitating such dissipation. Therefore, the fact that electron are heated much more in par- allel direction of 3D turbulence than 2D turbulence in Fig. 3.6(b) and 3.7 imply that, at least for the e = 0.10 cases considered here, it is the Landau resonance due to obliquely prop- agating whistlers which provides the stronger mechanism for damping of the turbulence. The more rapid transfer of whistler fluctuation energy to strongly oblique propagation illus- trated in Fig. 3.5 indicates that 3D whistler turbulence is more efficient at dissipating such turbulence and provides a plausible explanation for the more rapid decrease of the fluctuat- ing energy illustrated in Fig. 3.6(a). 3.4 Conclusions The first fully 3D PIC simulation of whistler turbulence in a homogeneous collisionless plasma with a uniform background magnetic fieldB o has been carried out. We impose an initial relatively isotropic spectrum of long-wavelength whistlers upon the system, with an initial e = 0:1, and compute the free decay of these modes to a broadband shorter wave- length regime of turbulence. The 3D system exhibits a forward cascade to shorter wave- lengths and yields transfer of fluctuation energy to a broadband, turbulent spectrum with wavevectors preferentially quasi-perpendicular to B o . Both the forward cascade and its consequent wavevector anisotropy are qualitatively similar to results obtained from previ- ous 2D PIC simulations of whistler turbulence. Quantitatively, however, in the 3D runs, tur- bulence develops a stronger anisotropic cascade more rapidly to a larger value, and exhibits faster rates of fluctuation energy dissipation than that in 2D simulations. Furthermore, the more modes that are used as initial conditions, the more rapid the cascade and the more strongly anisotropic the turbulence becomes. 59 The biggest finding of 3D whistler simulations is the spectral break at the inverse elec- tron inertial length scale. The reduced magnetic fluctuation energy spectrum of our 2D simulation shows ak 4:4 ? dependence atjk ? e j< 1, similar to the steep power-law behav- ior exhibited in earlier such simulations [79, 80]. In contrast, our 3D simulations show a less steep reduced spectrum of k 3:1 ? atjk ? e j< 1, a spectral break nearjk ? e j' 1, and a steeper spectrum withk 4:8 ? at shorter wavelengths, both qualitatively and quantita- tively similar to recent observations of short-wavelength turbulence in the solar wind. The much larger volume of perpendicular wavevector space in 3D appears to facilitate the trans- fer of fluctuation energy toward perpendicular directions, and may be responsible for the observed steepening of power spectra at electron scales. 60 CHAPTER 4: EFFECTS OF INITIAL FLUCTUATION AMPLITUDE ON WHISTLER TURBULENCE FORWARD CASCADE In this chapter, we describe the first ensemble of three-dimensional particle-in-cell plasma simulations of whistler turbulence representing a wide range of initial fluctuation ampli- tudes in a collisionless, homogeneous, magnetized plasma on which an initial spectrum of relatively long-wavelength whistler fluctuations are imposed. Introduction and initial con- dition of solar wind fluctuations are given in section 4.1 and 4.2, respectively. The detailed simulation results are presented in section 4.3. Discussion and conclusions are given in section 4.4. 4.1 Introduction Chapter 3 has already described the first fully 3D PIC simulation of whistler turbulence, using a single set of initial plasma parameters in a collisionless, homogeneous, magne- tized plasma model. The results not only reproduce the forward cascade and development of quasi-perpendicular wavevector anisotropy as shown in previous work, but also show much faster cascading rates, stronger dissipation rates and higher anisotropies than any of previous work. More importantly, for the first time, a magnetic spectral break near elec- tron scales is discovered in whistler turbulence, and the measured power-law dependence is close to the reported values from recent observations of short-wavelength solar wind turbulence [76, 4]. 61 We note that, previous 2D PIC simulations of whistler turbulence, unlike EMHD or some Hall MHD models [30, 26, 83, 31], did not show a “universal” spectrum with a fixed power-law dependence onk ? , but rather displayed reduced magnetic spectra which became less steep and more anisotropic as the initial amplitudes increased. Also, no spectral break near electron scales is reported from any of these works. To be able to reach more con- clusions and answer the question of whether a “universal” spectrum exists as opposed to the results from 2D PIC simulations, we here present the first ensemble of 3D whistler turbulence simulations in a wide range of initial fluctuation amplitudes in a collisionless, homogeneous, magnetized solar wind plasmas. 4.2 Initial conditions The plasma is homogeneous with periodic boundary conditions in all three directions. At t = 0 we impose a three-dimensional spectrum of relatively long wavelength fluctuations of same amplitude which each satisfy the linear dispersion properties of whistler waves. The initial spectrum is approximately isotropic in the sense that tan 2 ( B )' 1:4, which is relatively close to the value of unity which would obtain for a strictly isotropic spec- trum. The initial spectrum is the same as that of Run 3 in chapter 3. Thus there are 42 modes in the k y -k z plane corresponding to k y e = 0:1227n y (n y = 0;1;2;3) andk z e = 0:1227n z (n z =1;2;3); 36 modes in thek x -k z plane corresponding to k x e = 0:1227n x (n x =1;2;3) and k z e = 0:1227n z (n z =1;2;3); and an additional 72 modes corresponding to rotations of thek x -k z modes through 45 o and 135 o about the k z axis. The single parameter which we vary over this ensemble is the initial energy density of the whistler magnetic fluctuations: o = 0:02; 0:05; 0:1; 0:2; 0:5 and 1.0. The other plasma and computational parameters are the same as those in chapter 3, as listed in Table 3.1 and 3.2, respectively. 62 4.3 Particle-in-cell simulations and discussions We have carried out an ensemble of 3D PIC simulations of the free decay of whistler tur- bulence by forward cascade in a homogeneous plasmas. In each case, the initial ensemble of relatively long-wavelength whistler fluctuations undergoes a forward cascade to shorter wavelengths which is gyrotropic in the perpendicular plane as illustrated in Fig. 3.3, but bears the usual quasi-perpendicular wavevector anisotropy as shown in Fig. 4.3. Figure 4.1(a) shows the normalized total fluctuating magnetic field energy density and Figure 1(b) illustrates the normalized parallel (solid lines) and perpendicular (dashed lines) electron kinetic energies as a function of simulation time for all runs of the ensemble. The fluctuating magnetic field energy decreases with time; this decrease corresponds to an increase in the electron thermal energy with the parallel electron temperatures gaining considerably more energy than the perpendicular electron temperatures. Comparison of the o = 0.1 runs shows that, at late times, T ke =T ?e ' 1.3 in 2D whereas T ke =T ?e ' 1.6 in 3D. We interpret the greater dissipation efficiency of the 3D turbulence as due simply to the greater number of obliquely propagating modes. Figure 4.1(c) and (d) illustrate respectively the normalized magnetic energy damp- ing rate djBj 2 =dt=jBj 2 and the normalized electron parallel temperature heating rate dT ke =dt=jBj 2 for each run. Both the heating rate at all values of o as well as the damping rate at smaller values of o are weak functions of the initial fluctuation energy, consistent with the hypothesis of this dissipation being due primarily to linear damping, at least for o 0:2 cases. Further evidence in support of this hypothesis follows from the damping rate of whistler fluctuations at quasi-perpendicular propagation in a e = 0:1 plasma being of the order of = e =0:002 [79], consistent with the late-time field dissipation rate illustrated in Fig. 4.1(c). The dashed line in Fig. 4.1(c) isdjBj 2 =dt=jBj 2 =0:004, which corresponds to = e =0:002. Figure 5.1 shows that, at kc=! e < 1, linear the- ory predicts that the Landau resonance at oblique propagation is a stronger mechanism 63 0 0.2 0.4 0.6 0.8 1 | B| 2 /| B| t=0 2 0 50 100 150 200 250 300 350 400 450 0 2 4 6 T e /T e,t=0 | e |t −20 −15 −10 −5 0 5 x 10 −3 d| B| 2 /dt/| B| 2 0 50 100 150 200 250 300 350 400 450 −5 0 5 10 15 x 10 −3 dT ||e /dt/| B| 2 | e |t o = 0.02 o = 0.05 o = 0.1 o = 0.2 o = 0.5 o = 1.0 (a) (b) (d) (c) Figure 4.1: Simulation time histories of (a) the total fluctuating magnetic energy densityjB(t)j 2 =jBj 2 t=0 , (b) the parallel (solid lines) and perpendicular (dashed lines) components of electron temperatureT ke =T keo andT ?e =T ?eo whereT eo T e;t=0 , (c) the normalized magnetic energy damping ratedjBj 2 =dt=jBj 2 and (d) the normalized electron parallel temperature heating rate dT ke =dt=jBj 2 . Here the light blue lines represent the o = 0:02 run, red lines represent the o = 0:05 run, dark blue lines represent the o = 0:1 run, green lines represent the o = 0:2 run, black lines represent the o = 0:5 run, and magenta lines represent the o = 1:0 run. The dashed line in (c) isdjBj 2 =dt=jBj 2 =0:004. 64 1 5 10 15 tan 2 B (a) wavevector anisotropy factor 0 50 100 150 200 250 300 350 400 450 0.4 0.8 1.2 1.6 T ||e / T e | e |t (b) electron temperature anisotropy o = 0.02 o = 0.05 o = 0.1 o = 0.2 o = 0.5 o = 1.0 Figure 4.2: Simulation time histories of (a) the average wavevector anisotropy factor tan 2 B , and (b) the average electron temperature anisotropyT ke =T ?e from the various o runs. The colors correspond to values of o are as labeled. for whistler damping than the cyclotron resonance. Thus the Landau resonance could be expected to contribute to most of the linear dissipation shown in low o cases. ThejB(t)j 2 damping rate at early times of relatively large o = 0:5 and 1.0 increases with the initial fluctuating field energy, indicating that fully nonlinear processes play a major role here. Figure 4.2(a) illustrates the wavevector anisotropy factor tan 2 B as a function of sim- ulation time for each run. The 3D turbulence ranging over 0:02 o 0:2 becomes more anisotropic with increasing initial fluctuation energy o , similar to the 2D simulations 65 of [79]. But quantitatively, for each value of o in this range, the 3D simulation is much more strongly anisotropic than the corresponding 2D simulation. Furthermore, the large amplitude o = 0:5 and 1.0, 3D runs show a decrease in the wavevector anisotropy with increasing initial fluctuation amplitude. Similar trend could be observed in Fig. 4.2(b). The electron temperature anisotropy continuously increase as initial fluctuation increases until o reaches 0.5 where the dominating dissipation mechanism change to fully nonlin- ear processes. In quasi-linear regime, e.g., o 6 0:2, the preferential electron temperature anisotropyT ke =T ?e indicates that it is the electron Landau damping plays a major role in fluctuation energy dissipation. Figure 4.2, as well as Fig. 4.1(c), further suggest two distinct regimes in terms of the response to changes in o . The monotonic increase of the wavevector anisotropy with increasing o at 0:02 o 0:2 suggests linear (Landau) damping as the dissipation mechanism, whereas the decreasing anisotropies with further increases in o indicates fully nonlinear plasma processes are dominating in this regime. In each simulation the wavevector anisotropy rapidly increases at early times, as the for- ward cascade pumps magnetic fluctuation energy into the short wavelength part of the spec- trum where there is very little initial energy. At roughlyj e jt' 100, this process has estab- lished a quasi-steady spectrum with the usual quasi-perpendicular wavevector anisotropy as illustrated in the early-time panels of Fig. 4.3(a-c). At later times, electron dissipa- tion begins to take energy out of the short-wavelength part of the spectrum, particularly at oblique propagation angles, giving rise to the four-pointed star shape of the spectrum as shown in the late-time panel of Fig. 4.3(d). So we interpret the time development of the wavevector anisotropy as due to two different processes: at early times, the wave-wave processes of the forward cascade drive the anisotropy, whereas at late times wave-particle dissipation removes fluctuation energy so as to further enhance this anisotropy. A time series with detailed early times of 3D magnetic energy spectra from the simulations of o = 0:05; 0:1; 0:2 is illustrated in Fig. 4.4. 66 Figure 4.3: k z k y cuts from 3D wavevector iso-surfaces of magnetic fluctuation energy density at four times as labeled during the simulation of o = 0.2. 67 ε o = 0.2 ε o = 0.1 ε o = 0.05 Figure 4.4: Time series of 3D wavevector iso-surfaces of magnetic fluctuation energy density as labeled during the simulations of o = 0:05; 0:1 and 0.2. 68 10 −1 10 0 10 −8 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 |k e | E(k )/B 2 o o = 0.02 o = 0.05 o = 0.1 o = 0.2 o = 0.5 o = 1.0 noise level noise level k −2.8 k −5.0 Figure 4.5: k ? reduced magnetic fluctuation energy spectra atj e jt = 134 from the various o runs. The dashed lines represent power law functions as labeled for comparison against the simulation results. The colors correspond to values of o are as labeled. Figure 4.5 plots thek ? reduced magnetic fluctuation energy spectrumE(k ? ) from each run at e t = 134. In contrast to the previous 2D whistler simulations, there are clearly distinct breaks in each reduced spectrum atk ? e ' 1, with considerably steeper slopes at still shorter wavelengths, especially for o 0:1. Increasing initial fluctuation amplitudes over 0.02 o 0.5 yields a monotonic increase in the energy of the cascaded fluctuations, as well as to a decrease in the slope of the spectrum at bothk ? e > 1 andk ? e < 1. However, while o continue to increase to much larger values, the spectral indices at bothk ? e < 1 andk ? e > 1 seem asymptotically reach constant values, e.g.,' -2.8 at k ? e . 1. More importantly, the value = -2.8 is exactly the same as the spectral index 69 measured at recent solar wind observations [76, 4], and the value is only achieved at very large amplitude, i.e., o = 0:5 and 1.0, where fully nonlinear plasma processes begin to dominate the whistler turbulence cascades as discussed previously in Fig. 4.1, 4.2. This may indicate that in the actual solar wind condition, nonlinear dissipation plays a substan- tial role in shaping the magnetic energy spectra. Linear damping mechanism is partially contribute to the steady state solar wind spectra. To support this argument, we note that, Fig. 7(a) of [79] and arguments therein confirm that, linear Landau damping of whistlers at sufficiently oblique propagation scales as k 2 to and beyond k e ' 1. If the wave amplitude is sufficiently small that the turbulence could be described as an ensemble of normal plasma modes with properties derived from linear theory, such as the magnetic fluctuation B / e i(kr!rti t) , the linear Landau damping k 2 will indicate the magnetic energy spectra scaled asE(k)/ e 2 t / e k 2 t . Therefore linear damping alone cannot explain this short-wavelength spectral break, nor the power law spectral index -2.8 and -4.5 at both side of spectral break. Nonlinear wave-particle interactions and/or dissipation associated with reconnections at intermittent current sheet structures are the more likely source of this steepening with nonlinear wave-wave interactions also playing a possible role. Whatever the cause, it is clear that our 3D simulations of whistler turbulence have captured the further steepening of very short-wavelength spectra as reported by recent observations, and in our view, it is a potential hallmark of whistler turbulence at such short wavelengths. Our simulations do not show a strongly variable power law spectra index, whereas the spectral indices are preferable to be interpreted as asymptotically reaching a fixed value similar to those of the solar wind observations, at least for e = 0:1 cases considered in this chapter. The results also suggest that a large initial amplitude is the spectral index limit in which a comparison may be drawn between EMHD computational models [83] and PIC simulations. To make this comparison more complete, whistler turbulence simulations at considerably smaller values of e to reduce Landau damping are carried out in chapter 5. 70 −6 −4 −2 0 2 4 6 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 v ! / √ 2v te o f e (v || ) −6 −4 −2 0 2 4 6 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 v ⊥ / √ 2v te o f e (v ) Figure 4.6: Reduced electron velocity distributions as functions of v k (top panel) and v ? (bottom panel) atj e jt = 447 from the various o runs. The colors of the solid lines correspond to values of o as stated in the caption to Fig. 4.1, with the grey line representing the initial distributions att = 0. The dashed lines represent the difference between the late- time and the initial distribution for the case o = 0:1 normalized by the initial distribution. Figure 4.6 shows the reduced electron velocity distributions as functions f e (v k ) and f e (v ? ) atj e jt = 0 and 447 for each run. As one would expect, increased initial fluc- tuation energy (i.e., larger o ) corresponds to greater electron heating. Electron heating is stronger in directions parallel to B o than in perpendicular directions for low o cases 71 where linear damping is a primary dissipation mechanism. This is again indicative of elec- tron Landau damping being a primary linear dissipative mechanism. However, preferential parallel electron heating may also be due to parallel electric fieldE k associated with recon- nections at current sheets. This is supported by the fact that higher o cases where nonlinear dissipation mechanism dominates still showT ?e =T ke > 1, as illustrated in Fig. 4.2(b). Note that an increasing o leads to a flattening of the subthermal part of both velocity component distributions, as well as to an enhancement of the suprathermal tails of those distributions. This may be related to the work of [75] who used analytic theory to predict that kinetic Alfv´ en waves should cause a flattening off e (v k ), and suggested that a similar flattening as well as a suprathermal tail formation should arise as a consequence of electron scattering by magnetosonic-whistler turbulence. The flattening of f e (v k ) and f e (v ? ) is an irreversible change in the electron velocity distribution for relatively large values of o because Fig. 4.1(a) shows the fluctuating magnetic fields are strongly dissipated at late simulation times. Thus there can be no recovery to the large-amplitude initial conditions. 4.4 Conclusions We have here described the first ensemble of variable initial fluctuation amplitude o of 3D particle-in-cell plasma simulations of whistler turbulence. The computational model rep- resents a collisionless, homogeneous, magnetized plasma on which an initial spectrum of relatively long-wavelength whistler fluctuations are imposed. The simulations follow the temporal evolution of the system as it decays into a broadband, anisotropic, turbulent spec- trum at shorter wavelengths via a forward cascade. At relatively weak o 0:2, increasing o leads to more anisotropic magnetic wavevector spectra, and stronger heating of electrons in directions parallel to the background magnetic field. The anisotropies in 3D are stronger than in comparable 2D runs for all these o values; we attribute this difference to the larger 72 volume of perpendicular wavevector space in the 3D runs, which facilitates the transfer of fluctuation energy toward the perpendicular. The ensemble confirms that there is a characteristic spectral break near inverse electron inertial length in 3D whistler turbulence simulations, with both sides of the break show- ing a power law dependance, independent of the initial fluctuation amplitudes. With o keeps increasing, the power law index is asymptotically approaching -2.8, a value that is reported from recent solar wind observations [76, 4]. Our simulations results do not show a strongly variable power law spectra index, but rather an asymptotic “universal” value of -2.8 very close to the solar wind measurements, at least for e = 0:1 cases considered here. More importantly, our simulation results suggest two regimes of distinct dominating dissipation in response to various initial fluctuation amplitude. At lower o 0:2, linear Landau damping is the primary dissipation mechanism, whereas at higher o 0:5, non- linear processes is more likely a dominant dissipation mechanism. Note that the spectra index’s magic number -2.8 is only approached at higher o 0:5. This may indicate that in the actual solar wind condition, nonlinear dissipation plays a substantial role in shap- ing the magnetic energy spectra. However, our simulations presented in this chapter are at e = 0:1 while typical solar wind has e 1:0. To be able to reach definitive conclusions, further investigate on e dependance are carried out in chapter 5. 73 CHAPTER 5: EFFECTS OF ELECTRON BETA ON WHISTLER TURBULENCE FORWARD CASCADE In this chapter, we describe an ensemble of three-dimensional particle-in-cell simulations of whistler turbulence at variable initial e on a collisionless, homogeneous, magnetized plasma model. The value of e covers a wide range, representing conditions from real- istic solar wind e = 1:0 to ideal electron magnetohydrodynamics e = 0:01. With an introduction in section 5.1, initial plasma parameters for variable e are given in section 5.2. In section 5.3, we present whistler mode dispersion relations based on linear theory calculations. The detailed simulation results are discussed and the comparison with dif- ferent initial e and fluctuating amplitudes are presented in section 5.4. Suggestions on future observational work to resolve the wave identity in short-wavelength turbulence is also given in section 5.5. Section 5.6 introduces the frequency mismatch concept based on linear theory to test what conditions of the nonlinear wave-wave coupling favors a forward cascade and wavevector anisotropy of whistler turbulence. Finally, section 5.7 summarizes the new findings. 5.1 Introduction Fully 3D PIC simulations of whistler turbulence with variation of initial fluctuating ampli- tudes have been described in chapters 3 and 4. In these simulations a spectrum of rela- tively long wavelength whistlers (0:1< k k;? e < 0:3) is imposed as an initial condition, and the fluctuations are allowed to freely decay in time. These simulations exhibit the 74 development of a broadband spectrum of turbulence via a forward cascade, in which rela- tively long wavelength fluctuations transfer their energy via nonlinear processes to shorter wavelength modes. Furthermore, the short wavelength fluctuations preferentially develop quasi-perpendicular propagation, that is, withk ? k k . However, these 3D PIC computations were carried out only with initial e = 0:1. Here we describe 3D simulations of the forward cascade of whistler turbulence for three separate cases: the initial conditions of e = 0.01, 0.10, and 1.0 with the 1.0 case representing the condition most typical of solar wind plasmas, and the 0.01 case representing the condition approximately close to ideal electron magnetohydrodynamics. This parametric variation embraces a very broad range of plasma physics effects; here we describe consequent vari- ations in the forward cascade rate, the wavevector anisotropy, the electron heating rate and the electron anisotropy in both the time-space domain and the frequency-wavevector domain. We also repeat fluctuating amplitude variation as those in chapter 4 for lower and higher e , and compare characteristick ? reduced magnetic energy spectra across both e and initial amplitudes. It provides fresh information to help determine whether whistler fluctuations contribute substantially to electron-scale turbulence in the solar wind. We note that the 2D simulations of [78] have been conducted for different values of initial e . Here we choose an alternative approach to compare different initial e . 5.2 Initial conditions For all our simulations, the plasma is collisionless and homogeneous with periodic bound- ary conditions. The computational parameters are the same as in our previous 3D simula- tions, as listed in Table 3.2. The wavevectors of the initial spectrum of whistler fluctuations imposed upon the system are also the same as those described in chapters 3 and 4. Note that our initial spectra wavevector are normalized to electron inertial length e . Therefore, these wavevectors normalized in terms of the electron gyroradius e vary with e because 75 Table 5.1: Different plasma parameters in three variable e simulations e 0.01 0.1 1 ! e = e 0.71 2.24 7.07 o 0.02 0.2 2.0 e 2.0 2.0 2.0 k e =k e p e . The initial plasma conditions are mostly the same as those used in chapters 3 and 4: m p =m e = 1836, T e = T p and v 2 te =c 2 = 0.01. To study the consequences of e variations, we vary the magnitude ofB o to obtain three cases: ! 2 e = 2 e = 50 implying e = 1.0,! 2 e = 2 e = 5 implying e = 0.1, and! 2 e = 2 e = 0.5 implying e = 0.01. We note that the variation of e corresponds to a decrease of! e = e , however, linear dispersion calculations shown in Fig. 5.1 imply that as long as the value! e = e & 1, properties of whistler modes remain unchanged. These different values are highlighted in Table 5.1, and other common normalized parameters can be found in Table 3.1. As in our previous PIC simulations of whistler turbulence, each initial mode is given the same amplitude, and e is increased through a decrease inB o . However, in the 2D sim- ulations of [78] the dimensionless parameter characterizing the initial magnetic fluctuation energy was o , defined by k jBj 2 t=0 =8 = o B 2 o =8. At constant o , an increase in e due to a reduction ofB o implies a decrease in the initial dimensional magnetic fluctuation energy; as a consequence the dimensional electron heating rate necessarily decreases, as shown by their 2D simulations. Here we take a different approach, which we believe is a better way to compare results of turbulent heating due to e variations. We choose the initial dimensional magnetic fluctuation energy density to be fixed asB o and e are varied. Thus, we define e by k jBj 2 t=0 =8 = e n e k B T e;t=0 . For the simulations described here e = 2.0, so if e = 0.01, then o = 0.02; if e = 0.1, then o = 0.2; and if e = 1.0, then o = 2.0. 76 Figure 5.7 shows that the heated electron velocity distributions become strongly non- Maxwellian at very low values of initial e , so instead of using the electron temperature as a diagnostic of electron heating we use the electron kinetic energy density which we define as K e m e 2 Z d 3 vv 2 f e (v) (5.1) with parallel and perpendicular components K ke m e 2 Z d 3 vv 2 k f e (v) (5.2) and K ?e m e 2 Z d 3 vv 2 ? f e (v) (5.3) 5.3 Linear analysis Although turbulence could not be fully described as an ensemble of weakly interacting incoherent normal modes, the study of electromagnetic waves in whistler modes regime could still provide predictions and explanations on the characteristic cascades and energy dissipations of whistler turbulence. As we described earlier in section 1.1.1, the whistler wave dispersion relation is very complex, such as Eq. 1.1. One usually has to use numer- ical approaches to solve the full kinetic dispersion equation for whistler mode. Some recent studies have utilized the cold plasma approximation to yield analytic expressions for whistler dispersion in turbulence applications [59, 58, 27]. However, thermal plasma effects produce substantial differences in whistler dispersion properties for the three differ- ent values of e which we consider here. Therefore, we here use numerical solutions of the full linear kinetic dispersion equation for Maxwellian velocity distributions, as described in section 2.1, to obtain whistler dispersion properties in the following analysis. 77 0 1 2 3 4 5 6 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 /| e | /| e | e = 0.01 0 0.2 0.4 0.6 0.8 −0.05 −0.04 −0.03 −0.02 −0.01 0 /| e | e = 0.01 0 0.5 1 1.5 2 2.5 3 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 /| e | /| e | e = 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 −0.2 −0.15 −0.1 −0.05 0 /| e | e = 0.1 0 0.25 0.5 0.75 1 1.25 1.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 /| e | /| e | e = 1.0 |kλ e | 0 0.1 0.2 0.3 0.4 −0.25 −0.2 −0.15 −0.1 −0.05 0 /| e | e = 1.0 ω/|Ω e | = 0 o = 55 o = 80 o or 70 o Figure 5.1: Linear theory results for whistler dispersion as a function of wavenumber jk e j for three propagation angles at three e as labeled. Green, blue, red lines denote = 0 o ; 55 o , and 80 o or 70 o for e = 1:0 cases. Left column shows dispersion solutions for real frequency! and damping rate . Solid, dotted, dash dotted, dashed lines represent solutions for! 2 e = 2 e = 2500, 50, 5 and 0.5. Right column shows dispersion solutions for vs.!. The values! 2 e = 2 e correspond to those listed in Table 5.1. 78 Figure 5.1 left column illustrates whistler dispersion for parallel ( = 0 o ), oblique ( = 55 o ) and quasi-perpendicular ( = 80 o ) propagation angles at e = 0.01, 0.10 and 1.0. The green lines denote = 0 o , the blue lines denote = 55 o , and the red lines denote = 80 o or 70 o for e = 1:0 cases. The positive frequency range indicates the dispersion equation solution for!(k), and the negative portion indicate the damping rates . As we described earlier, the value of! e = e in our simulation is much smaller than the value in the solar wind. However, our numerical solutions indicate that the characteristics of whistler dispersion is insensitive to ! e = e as long as ! e = e & 1. Our choices of ! 2 e = 2 e = for the dispersion solver are 0.5, 5, 50, and 2500. The numerical solutions of the dispersion relations based on these choices are superimposed on Fig. 5.1. The solid lines represent ! 2 e = 2 e = 2500, the dotted lines represent ! 2 e = 2 e = 50, the dash dotted lines represent ! 2 e = 2 e = 5, and the dash lines represent! 2 e = 2 e = 0:5. As all of the three panels in figure show, for each different e and choice, the dispersion relations for! 2 e = 2 e & 1 are almost the same. When ! 2 e = 2 e further decreases to 0.5 (dashed lines), the dispersions for !(k) diverge from the higher ! 2 e = 2 e solutions, especially for the oblique propagating waves. The ! tends to be smaller. Nevertheless, the damping rates only experience a smaller amount of change for e = 0:01 cases. In our simulations, only e = 0:01 cases have a ! e = e = 0:71 1; so we expect very small amount of difference to ! e = e > 1 cases for e = 0:01 simulations. For e = 0:1; 1:0 simulations, the characteristics of whistlers should remain unchanged compared to! e = e 1 solar wind conditions. The left column clearly shows that the damping rate of high e = 1:0 cases are much more stronger than lower e cases. For example, high e = 1:0 case reaches a relatively large damping rate =j e j =0:06 whenjk e j< 0:8, while lower e cases need to extend wavevector to at leastjk e j> 1:2. There are two types of wave-particle damping involved here. At oblique propagationkB o 6= 0, electron Landau damping as described in section 1.1.2 is the major wave damping mechanism. Note that the damping rate for whistlers at low e = 0:01 regime is very small so that Landau damping could almost be neglected at 79 jk e j < 2:0. AtkB o = 0, there is no Landau damping for electromagnetic waves, so the increase in damping at successively higher e is due to the whistler fluctuations moving into cyclotron resonance with suprathermal and thermal electrons. The electron cyclotron damping is essentially zero atjk e j . 1:3; 0:8; 0:32 for e = 0:01; 0:1; 1:0, respectively, but has a sudden onset after that. In contrast, electron Landau damping, which arises only at > 0 o , has no clear onset, but increases gradually with increasing wavenumber across the full range of wavenumbers. Figure 5.1 right column shows the damping rate as a function of real frequency !. We only keep dispersion solutions that correspond to the simulations cases, i.e., the left figure represents e = 0.01 and! 2 e = 2 e = 0.5, the middle figure represents e = 0.1 and ! 2 e = 2 e = 5, and the right figure represents e = 1.0 and ! 2 e = 2 e = 50. Clearly, quasi- perpendicular and oblique propagating modes, for the same!, always have much bigger damping rate than parallel propagating modes for e = 0:01; 0:1 cases and for e = 1:0 at!=j e j. 0:2. The turning point of e = 1:0 case marks the rapid increase of its electron cyclotron damping in parallel direction. Figure 5.1 also indicates that quasi-perpendicular and very oblique modes are always in lower frequencies regime compared to quasi-parallel modes. Figure 5.2 illustrates the 2D dispersion relation contours as functions ofjk k j andjk ? j for whistler modes at three e values. The left column shows the real frequency!, while the right column shows the damping rate . The parameters chosen correspond to those used in the simulations. Clearly, at e = 0:01, the damping is so small that the numerical solutions for! could be found everywhere in wavevector domain that we are interested in, i.e.,jk e j 3:0. In contrast, e = 1:0 case have much larger damping at smallerjk e j. For example, = e = 0:3 atjk k e j = 1 for e = 1:0 case, whereas the value for lower e cases are less than 0.02. Because the damping of wave amplitude is too big for e = 1:0 case, our dispersion solver could not find any valid solutions in whistler mode range beyond jk k e j > 1:2 andjk ? e j > 2, resulting in a much smaller dispersion contour comparing 80 Figure 5.2: 2D contours for whistler dispersion as functions of wavenumberjk k j andjk ? j at e = 0.01, 0.1 and 1.0 as labeled. The left column’s color represents the real frequency !, and the right column’s color represents the damping rate absolute valuej j. 81 with lower e cases. We will frequently refer to Figs. 5.1 and 5.2 when we discuss the simulations results in the next section. 5.4 Particle-in-cell simulations and discussions Simulation results and discussions are organized into categories of integrated quantities (e.g., fluctuation energy history), electron properties, wavevector spectral properties and frequency spectral properties. Most of these results are from three characteristic runs rep- resenting three different e as highlighted in Table 5.1. Unless otherwise noted, we refer to these three runs in the discussion. To analyze certain properties of whistler turbulence, we conducted further runs, such as varying initial fluctuating amplitude, initial spectrum shape, etc. We will note these special runs when it is necessary. Figure 5.3(a) shows the total energy density (dashed lines) and the total magnetic fluc- tuation energy density (solid lines) as a function of time. It demonstrates that the total energy is well conserved in each of the simulations, that the total fluctuating field energy is dissipated with time, and that the rate of dissipation increases dramatically with increas- ing initial e . The total energy variation is less than 0.4% throughout each of the three simulations. Figure 5.3(b) and (c) illustrate the cascaded magnetic fluctuation energy (which corre- sponds to wavenumbers satisfying 0:65jk e j 3:0) as a function of time in logarithmic and linear scales, respectively. They show that there are two distinct stages for each sim- ulation. During early times wave-wave interactions carry magnetic fluctuation energy to shorter wavelengths, where there is initially only a low level of thermal fluctuations present. The rate of forward cascade in wavenumber is a sensitive function of e , increasing rapidly with this parameter, as more clearly shown in Fig. 5.3(c). The early-time stage ends at the maximum value of the cascaded energy; at this time the overall damping rate becomes 82 0 500 1000 1500 2000 2500 3000 3500 4000 10 −1 10 0 | B| 2 /8 K e,t=0 (a) magnetic fluctuation energy density 10 0 10 1 10 2 10 3 10 −4 10 −3 10 −2 10 −1 | B| 2 /8 K e,t=0 (b) cascaded magnetic fluctuation energy density in log scale 0 200 400 600 800 1000 0.02 0.04 0.06 0.08 0.1 | B| 2 /8 K e,t=0 | e |t (c) cascaded magnetic fluctuation energy density in linear scale e = 0.01 e = 0.1 e = 1.0 Figure 5.3: Simulation results as functions of time for three runs with initial values of e as labeled. (a) Total energy density (dashed lines) and total magnetic fluctuation energy density (solid lines). (b) Cascaded magnetic fluctuation energy density (summed over 0:65jk e j 3:0) normalized to initial electron kinetic energy density. (c) An enlarged (b) in linear scale of both axes. 83 faster than the overall cascade rate so that dissipation dominates the late-time stage of the simulations, as indicated by the temporal decrease in the cascaded magnetic energy density. To investigate the damping rate transition, we perform extra simulations of variable initial amplitude e at e = 0.01 and 1.0, with the combination of variable e at e = 0:1 cases from chapter 4. Figure 5.4 illustrates the normalized magnetic energy damping rate djBj 2 =dt=jBj 2 for e = 1:0; 2:0; 5:0 at e = 0.01, 0.1 and 1.0. Consistent with Fig. 4.1(c), for all three e sets, the damping rate at smaller values of e are weak functions of the initial fluctuation energy, indicating this dissipation is primarily due to linear damping for smaller amplitude cases, such as e 2:0. From linear theory results in Fig. 5.1, the damping rate of whistler fluctuations at quasi-perpendicular propagation are of the order of = e = -0.0004, -0.002, -0.01 for e = 0.01, 0.1 and 1.0 plasmas, respectively. The cor- responding values are shown in dashed lines in each panel. The late-time field dissipation rate from each e simulation sets are consistent with these damping rate from the linear theory results, again implying linear dissipation as primary dissipation for e 2:0 of all e cases. ThejB(t)j 2 damping rates of all e cases start to increase at early times as e continues to increase, indicating that nonlinear processes gradually take over the major role in energy dissipation. Figure 5.5 shows how the wavevector anisotropy of the cascaded spectrum varies in time. These figures imply two important properties of the wavevector anisotropy. First, increasing e from 0.10 to 1.0 makes the late-time wavevector distribution more nearly isotropic; Fig. 3 from the 2D PIC simulations of [78] shows a similar trend. Second, increasing e from 0.01 to 0.10 increases the early-time rate of wavevector anisotropy development. Comparison of Fig. 5.3 and 5.5 shows the clear correlation between the electron scattering (i.e. dissipation) time and the wavevector anisotropy development time. In other words, both wave-wave and wave-particle interactions contribute to the develop- ment of the wavevector anisotropy here. 84 0 200 400 600 800 1000 1200 1400 −6 −4 −2 0 2 4 6 x 10 −3 d| B| 2 /dt/| B| 2 e = 0.01 0 50 100 150 200 250 300 350 400 450 −12 −10 −8 −6 −4 −2 0 x 10 −3 d| B| 2 /dt/| B| 2 e = 0.1 0 10 20 30 40 50 60 70 80 −0.05 −0.04 −0.03 −0.02 −0.01 0 d| B| 2 /dt/| B| 2 | e |t e = 1.0 e = 1.0 e = 2.0 e = 5.0 Figure 5.4: Time histories of the normalized magnetic energy damping rate djBj 2 =dt=jBj 2 for the three simulation sets with initial values of e and e as labeled. The dashed lines in upper, middle and lower panels are djBj 2 =dt=jBj 2 = -0.00002, - 0.004, -0.02, respectively. 85 0 0.2 0.4 0.6 0.8 1 0 500 1000 1500 2000 2500 3000 3500 0 10 20 30 tan 2 B | e |t e = 0.01 e = 0.1 e = 1.0 Figure 5.5: Wavevector anisotropy factor of the cascaded magnetic fluctuation energy as a function of time during the three simulations with initial e = 0.01, 0.10, and 1.0. 1 1.2 1.4 1.6 1.8 2 K e /K e,t=0 (a) electron kinetic energy density 0 500 1000 1500 2000 2500 3000 1 2 3 4 (K ||e /K e ) / (K ||e /K e ) t=0 | e |t (b) electron kinetic energy anisotropy e = 0.01 e = 0.1 e = 1.0 Figure 5.6: (a) Normalized electron kinetic energy density and (b) normalized electron kinetic energy anisotropy as functions of time during the three simulations with initial e = 0.01, 0.10, and 1.0. 86 Figure 5.6(a) illustrates the electron kinetic energy density as a function of time for each run; as with the field dissipation, the electron heating rate increases with e , although the asymptotic late-time electron kinetic energy density is approximately the same in all three cases. Figure 5.6(b) demonstrates the electron kinetic energy anisotropy as a function of time for the three runs; here the late-time anisotropy decreases with increasing e , the same general trend as exhibited by the late-time wavevector anisotropy in Fig. 5.5. It implies again, both wave-wave and wave-particle mechanism contribute to the growth of the wavevector anisotropy. Figure 5.7 illustrates the late-time reduced electron velocity distributions f e (v k ) and f e (v x ). As in Fig. 5.6(b),T k > T ? for each value of e , implying that the primary energy transfer channel of linear dissipation from whistler fluctuations to the electrons is via Lan- dau damping. Our interpretation of Fig. 5.7 is based on our solutions of the full kinetic linear disper- sion equation, as illustrated in the linear theory plots of section 5.3 above. At e = 0.01 the Landau resonance factor e != p 2jk k jv te is much greater than unity over 0.3jk e j and whistler waves resonate only with the few particles on the tail of the electron velocity distribution, e.g., the nonresonant phase speed of Fig. 2.1 of [33]. Thus, our e = 0.01 simulation shows that it is these fast electrons which exhibit the strongest heating, forming a broad tail on both the v ? and especially on the v k reduced distributions. The whistler anisotropy instability driven by an electron temperature anisotropyT ? >T k yields similar enhanced tails onf e (v k ) at very low e [36]. As e increases from 0.01 to 1.0, our solutions of the dispersion equation show a relatively weak change in the real frequency, but there is an order of magnitude increase in the electron thermal speed, so that e typically becomes of order of or less than unity, resonating with the large number of electrons in the thermal and subthermal parts off e (v k ), e.g., the resonant phase speed of Fig. 2.1 of [33]. Thus our e = 0.10 and 1.0 simulations show electron heating concentrated in the thermal part of the electron velocity distribution, so that especially at e = 0.10 wave-particle interactions 87 0 0.1 0.2 0.3 0.4 v ! / √ 2v te o f e (v || ) 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 v ! / √ 2v te o 0 0.2 0.4 0.6 0.8 1 v || /v te,t=0 0 0.2 0.4 0.6 0.8 1 v || /v te,t=0 −6 −4 −2 0 2 4 6 0 0.1 0.2 0.3 0.4 v ⊥ / √ 2v te o f e (v ) −6 −4 −2 0 2 4 6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 v ⊥ / √ 2v te o 0 0 0.2 0.4 0.6 0.8 1 v /v te,t=0 0 0 0.2 0.4 0.6 0.8 1 v /v te,t=0 t=0 e =0.01, | e |t=3000 e =0.1, | e |t=1000 e =1.0, | e |t=300 Figure 5.7: Reduced electron velocity distributions as functions ofv k (top panels) and ofv ? (bottom panels) at late times during the three simulations with initial e = 0.01, 0.10, and 1.0. The distributions are shown in linear scales in the left panels, and logarithmic scales in the right panels. 88 lead to a flattening of the subthermal part off e (v k ), as pointed out by [58]. Note, however, our results differ from those of [58] in two potentially important ways. First, Fig. 5.7 shows that the flattening off e (v k ) is relatively strong only at relatively large values of e . Second, although Fig. 1 of [58] shows no significant electron heating atv v te , our simulations demonstrate substantial electron energy gain at thermal and suprathermal electron energies. Figure 5.8 illustrates k k k ? reduced magnetic energy spectra from all three runs at several representative times. For all three runs the initial ensemble of relatively long- wavelength whistlers drives a forward cascade of fluctuating field energy toward shorter wavelengths which, as shown in Fig. 3.3, has a gyrotropic wavevector distribution. The three left-hand panels of this figure show relatively early-time spectra with successively decreasing wavevector anisotropies in the sense ofjk ? j > jk k j with increasing e . In section 5.6 we argue that wave-wave coupling of three whistlers should yield a relatively isotropic spectrum at e = 1:0, whereas wave-wave interactions at smaller values of e should yield more anisotropic spectra in the sense ofjk ? jjk k j, consistent with these three panels. The other nine panels of Fig. 5.8 show that the wavevector anisotropies continue to increase during those times when dissipation dominates. Figure 5.9 plotsk ? reduced magnetic energy spectrumE(k ? ) from all three e sets of e = 2:0 and 5.0 at their respective early-time stage. Consistent with results in chapter 4, all spectra show a characteristic distinct break atjk ? e j' 1, with considerably steeper slopes at still shorter wavelengths, especially for o 0:1. Similar to Fig. 4.5, as e increases to much larger values where nonlinear processes dominate the dissipation, the spectral indices for power law dependance atjk ? e j < 1 andjk ? e j > 1 seem to reach asymptotically constant values. However, such constant value decreases as e increases, reaching'2;2:8;3:7 for e = 0:01; 0:1; 1:0, respectively, atjk ? e j. 1. How do the simulation energy spectra compare with short-wavelength solar wind obser- vations? Fittings using both the two-power-law break model and the exponential model, as introduced in section 1.2.2, are superimposed on Fig. 5.9 indicated by blue dashed lines 89 Figure 5.8: k k k ? reduced magnetic fluctuation energy spectra, log[jB(k)j 2 =8n e k B T e;t=0 ], at four times for the three simu- lations with initial values of e as labeled. 90 10 −1 10 0 10 1 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 E(k )/K e,t=0 |k e | e =0.01 10 −1 10 0 10 1 |k e | e =0.1 10 −1 10 0 10 1 |k e | e =1.0 e = 0.01, e = 2.0 e = 0.01, e = 5.0 e = 0.1, e = 2.0 e = 0.1, e = 5.0 e = 1.0, e = 2.0 e = 1.0, e = 5.0 k −4.3 k −1.5 e −k e k −5.0 k −5.5 k −2.6 e −k e k −3.4 e −k e k −3.7 k −2.8 k −2 Figure 5.9: k ? reduced magnetic fluctuation energy spectra at corresponding early time from the three simulation sets with initial values of e and e as labeled. The blue and red dashed lines represent, respectively, the break model’s power law functions and the exponential model’s exponential rollover functions, as labeled for comparison against the simulation results. 91 and red dashed lines. Note that we use e instead of e for normalization in the exponen- tial model/ k ? exp(jk ? e j) for the reason described later in section 5.5. In e = 0:01 panel, the reduced spectra is able to span over almost two magnitude fromjk ? e j 0:1 to jk ? e j 10, and the results clearly favor the two-power-law break model than the expo- nential model. In the low e regime, wave-particle damping is minimized, and the PIC simulations approach EMHD models. The spectral index2:0 atjk ? e j < 1 for the break model, which is in acceptable agreement with typical EMHD’s spectra scaling -7/3 of whistler turbulence [83]. In e = 0:1 panel, the reduced magnetic energy spectrum is fitted with the spectral index = -2.8 atjk ? e j< 1 for the break model, or with = -2.6 for the exponential model. These values are in good agreement with recent observations’ value, i.e., =2:8 for the power law dependance atjk ? e j< 1 [76] and =8=3 for the exponential rollover [3]. In e = 1:0 panel, where the value of e is typical solar wind value, the reduced spectra show steeper power law dependance at atjk ? e j< 1 (/k 3:7 ? ) or steeper exponential rollover (/ k 3:4 ? exp(jk ? e j)) dependance than the reported val- ues from the corresponding range of the solar wind observations. Although the e = 1:0 results’ wavevector spectral indices could not match the solar wind observations frequency spectral indices closely, our reduced energy spectra definitely demonstrate the characteris- tic spectral break near inverse electron inertial length scale, with asymptotically constant values of spectral indices in the large fluctuation amplitude solar wind where nonlinear pro- cesses dominate. These properties are are very similar to the observations, and are potential hallmarks of whistler turbulence at such short wavelength turbulence. Figure 5.10 plotsk ? reduced power spectrum of electron density fluctuations PSD(k ? ) from all three e sets of e = 2:0 and 5.0 at their respective early-time stage. The power law fittings (blue dashed lines) and the exponential rollover fitting (red dashed line) are super- imposed on the same figure. Like previous magnetic fluctuation energy spectra, electron density spectra reach asymptotic levels as initial fluctuation amplitude e increases to large values, and these density spectra also exist distinct breaks aroundjk ? e j' 1. 92 10 −1 10 0 10 1 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 PSD(k )/K e,t=0 |k e | e =0.01 10 −1 10 0 10 1 |k e | e =0.1 10 −1 10 0 10 1 |k e | e =1.0 e = 0.01, e = 2.0 e = 0.01, e = 5.0 e = 0.1, e = 2.0 e = 0.1, e = 5.0 e = 1.0, e = 2.0 e = 1.0, e = 5.0 k −2.8 k −4.0 k −4.0 k −4.6 k −2.8 e −k e k −2.9 k −7 / 3 Figure 5.10: k ? reduced power spectrum of electron density fluctuations PSD(k ? ) at corresponding early time from the three simulation sets with initial values of e and e as labeled. The blue and red dashed lines represent, respectively, the break model’s power law functions and the exponential model’s exponential rollover functions, as labeled for comparison against the simulation results. 93 Comparing with reduced magnetic energy spectra, the spectral indices of all three e cases for both power law and exponential rollover dependance are relatively close to each other. In the panels of e = 0.1 and 1.0, the spectral indices atjk ? e j < 1 is around -2.8 -2.9, which is close to the reported -2.75 0.06 from a very recent solar wind electron density observation [23]. We also apply [3]’s exponential model to fit e = 1:0’s density power spectrum, and it matches the model model very well. We note that [3]’s model was derived from the observation of magnetic fluctuations spectra at e 1:0. In the e = 0:01 panel, the spectral index atjk ? e j< 1 is about -7/3, which is in agreement with the predication of -7/3 for a pure whistler cascade based on a compressible Hall MHD model [2]. The overall steeper spectra at higher e cases are possibly due to the fluctuation energy being damped out by heavy wave-particle interactions. This has been illustrated both from linear theory predication in Fig. 5.1 and electron velocity distribution from the simulations 5.7. Our electron density spectra share a close similarity, both qualitatively and quantitatively, to the solar wind observation and MHD predications. These results further support the hypothesis that whistler fluctuations populate at high frequency solar wind turbulence near inverse electron inertial length scales. Note that Figs. 5.9 and 5.10 are wavevector spectra while observations [76, 69, 3, 23] usually provide frequency spectra. Observations could only directly measure the spectrum as function of spacecraft frame frequency, and with certain assumptions of dispersion rela- tions, such as Taylor frozen-in-flow hypothesis [92], may convert the observed frequency spectrum into a wavevector spectrum. Taylor’s hypothesis states that! sc =! +kV sw kV sw , where! sc is the the frequency in the spacecraft frame. This hypothesis requires ! <kV sw and may not be applied to fluctuations withv parallel to the solar wind flow faster thanV sw . 94 Let us check whether the whistler wave at typical solar wind condition satisfy the restriction or not. Assuming the background magnetic fieldB o is perpendicular toV sw , then the whistler perpendicular phase speed is v ? = ! r k ? = k k c=! e 1 + (1 + p e )(k 2 c 2 =! 2 e ) j e j ! e k k ? c (5.4) The typical solar wind parameters where e = 1:0 are listed in Table 3.1. From Fig. 5.5, the wavevector anisotropy factor tan 2 B at e = 1:0 is roughly 5, thus we assume jk ? =k k j = p 5. Forjk ? e j = 0:5; 1:0; 1:5, whistler perpendicular phase speedv ? = 328 km/s, 309 km/s, 246 km/s, which are all in same magnitude ofV sw = 300 km/s 600 km/s. Therefore, the Taylor’s hypothesis is not quiet satisfied for e = 1:0 oblique propagating whistlers, and there is no guarantee that such a frequency-wavevector conversion is able to obtain an accurate wavevector spectrum existing in the short-wavelength solar wind turbulence. Since comparisons of the specific spectral indices values from our simulations wavevector spectrum and observations frequency spectrum may not be appropriate, we now turn into frequency domain. Figure 5.11 shows the magnetic energy frequency spectra computed over the full tem- poral extent of each simulation for all three e cases. The spectra is averaged on 64 64 series of data points which are located on ay-z plane atx =L x =2. The simulation results, again, show the preference for the power law fitting than the exponential rollover fitting. All three e cases show a distinct power law dependance change around!=j e j 0:15. The reported electron-scale breaks from recent observations [76, 69] are 40 50 Hz. Using the typical solar wind value e = 1800 from Table 3.1, we obtain the observational spectral break range at!=j e j = 0:14 0:17. This is in considerable agreement with our simulations’ value of frequency spectral break. More importantly, in contrast to wavevector spectra in Fig. 5.9, the spectral indices of all three e cases are very similar, which is also implied from solar wind observations 95 10 −7 10 −5 10 −3 10 −1 | B| 2 /8 K e,t=0 e = 0.01 10 −7 10 −5 10 −3 10 −1 | B| 2 /8 K e,t=0 e = 0.1 10 −2 10 −1 10 0 10 −7 10 −5 10 −3 10 −1 | B| 2 /8 K e,t=0 / | e | e = 1.0 −8/3 −3.5 −8/3 −3.5 −2.48 −3.2 Figure 5.11: Magnetic energy frequency spectra with initial values of e = 2:0 and e as labeled. The solid blue, green and red lines represent the spectra at e = 0.01, 0.1 and 1.0., the dashed blue lines represent the power law functions as labeled for comparison against the simulation results, and the two vertical black dashed lines in each panel bound the fre- quency range of the initially imposed whistler fluctuations. The spectrum is averaged over 4096 series of data points computed over the full temporal evolution of each simulation. [3]. Besides, these values are much closer to those from solar wind observations, not only at frequencies lower than the break, but also at frequencies higher than the break, such that =8=3 before the break and =3:5 after the break. A typical observational 96 frequency spectrum is illustrated early in Fig. 1.1 from [76]. e = 1:0 case shows a slightly less steep spectra than those from the other two e cases and the observations. We interpret it as due to more numerical noise in the e = 1:0’s spectrum. The high e case is heavily damped as indicated in Fig. 5.3, and therefore, although the frequency spectra are computed over the full temporal evolution using same number of data points, more of the data points of the e = 1:0 case at late times are just thermal noise as compared to lower e cases. As far as we know, at the time of this dissertation being written, nobody else has reported any simulations results of frequency spectral break and spectral index in short- wavelength turbulence regime, including whistlers, kinetic Alfv´ en wave and current sheet simulations. The close similarity of frequency spectra from our simulations and solar wind observations significantly reinforces the interpretation that whistler fluctuations could be the substantial constituent of solar wind turbulence at higher frequencies and short wave- lengths. Whistler turbulence not only undergoes a strong cascade toward shorter wavelengths and anisotropic wavevectors, it is also subject to the strong dispersion of the whistler mode. The combination of these three properties can lead to the following: if wave energy is injected at intermediate frequencies, some is transferred to higher frequencies, but some moves toward lower frequencies near the lower hybrid frequency lh at quasi-perpendicular propagation [32]. This is also illustrated in Fig. 5.11. The initial frequency spectra for each run are approximately constant between the two dashed lines and are at the thermal noise level for fluctuations at frequencies both greater than those at the right-hand dashed lines and smaller than those at the left-hand dashed lines. Thus, this figure indicates a transfer of wave energy to both higher and lower frequencies. In particular the e = 0.01 and 0.10 spectra show decided peaks near the lower hybrid frequency!' 0:024j e j. It implies us the importance to perform corresponding data analyses not only in time domain, but also in frequency domain. 97 Figure 5.12: Continuous wavelet transform of magnetic fluctuation energy as functions of both time t and real frequency ! from a single data series in the simulation at three dif- ferent values of e as labeled. The 2D color contours represent the value of the coefficient log[CWT(B x ) 2 +CWT(B y ) 2 +CWT(B z ) 2 ], where CWT denotes continuous wavelet trans- form function. Here the mother wavelet chosen is Morlet wavelet for all three cases. The two parallel white dashed lines in each panel indicate the frequency range of the initially imposed whistler fluctuations. 98 To analyze the frequency cascade to both higher and lower frequencies, continuous wavelet transform, described in section 2.3.2, is used to illustrate the turbulence evolution process in both temporal and frequency domain [52]. Figure 5.12 shows the 2D coefficient contour of magnetic fluctuation energy log[CWT(B x ) 2 +CWT(B y ) 2 +CWT(B z ) 2 ] as func- tions of both timet and real frequency! at three e runs. Here CWT denotes a continuous wavelet transform. A higher value in the contour indicates a higher correlation between the wavelet and that section of the input signal at that frequency range. Although we only use a single data series to perform the wavelet transform due to expensive computational time, it clearly shows that frequency cascade to higher bands happens primarily at the early-time stage, and the cascade to lower frequencies continues at the late-time stage. It confirms that an enhanced peak near the lower hybrid frequency lh appears in the late time spectra of e = 0.01 and 0.10 cases, as also shown in Fig. 5.11. These enhanced peaks near lh indicate that nonlinear Landau damping (induced scattering) might play a role here at late times [32]. Figure 5.13, illustrates! versus wavenumber whistler turbulence dispersion relations, for the first time, at three different values of e from 3D PIC simulations. These dispersion plots are computed over the same time intervals as those used in Fig. 5.3(a). Here a specific perpendicular wave number is chosen asjk ? e j = 0:61. We choose this value because it is primarily in the cascaded wavevector range while the damping rates are still small enough to clearly show dispersion. Dispersion curves (black dashed lines) of plasma waves at this specifick ? from the numerical solution of the full linear kinetic dispersion equation are superimposed on each of the three panels. The vicinity of the fluctuation energy in the simulation and linear dispersion curve not only confirms the whistler wave range of turbulence, but also, more importantly shows the nonlinear cascaded frequency and wavevector properties that linear dispersion theory alone could not predict. Dispersion of lower-frequency magnetosonic waves in 2D PIC simulations of turbulence has been demonstrated earlier by [91, 46]. 99 Figure 5.13: Fluctuation dispersion (real frequency! as a function of parallel wavenumberjk k j) computed over the full temporal evolution of each simulation at three different values of e as labeled. The color contours represent simulation results of magnetic fluctuation energy, log[jB(k k ;k y )j 2 =8n e k B T e;t=0 ], and the black dashed lines represent numerical solutions of the full kinetic dispersion equation for whistler fluctuations. Herejk ? e j = 0.61 for all three values of e . 100 (a) k k view (b) k ? view (c) Bottom view (d) Aerial view Figure 5.14: 3D iso-surfaces of magnetic fluctuation dispersion (real frequency! as func- tions of parallel wavenumberk k and perpendicular wavenumberk ? ) computed over the full temporal evolution of the e = 0:1 simulation. The color contours represent simulation results of magnetic fluctuation energy, log[jB(k k ;k y )j 2 =8n e k B T e;t=0 ]. Each subfigure represents a different view perspective as labeled. To gain a better understanding of the fluctuation dispersion computed the simulation, we created a series of 3D iso-surfaces plots from the e = 0:1 case as shown in Fig. 5.14. Four view perspectives, i.e., k k , k ? , bottom and aerial views, are provided in panels (a), (b), (c) and (d), respectively. At roughlyk k e < 0:4 andk ? e < 0:5, there is an anomalous energy enhancement throughout the whole frequency range. Based on its wavevector range, 101 we believe this is a numerical effect due to the initially imposed fluctuation spectrum. It may affect the computed dispersion shape since the dispersion is calculated over the whole temporal evolution of the magnetic fluctuations. Nevertheless, we are more interested in the cascaded range of the dispersions. As indicated by Fig. 5.14, at higher frequency bands, the whistler turbulence wavevectors tend to be less anisotropic, while at lower frequency bands, the wavevectors are much more anisotropic. These detailed feature as functions of frequencies could not be obtained from any previous reduced magnetic energy spectra without performing a four dimensional (three spatial dimensions + one time dimension) Fourier transform. Figure 5.15 illustratesk k k ? reduced magnetic energy spectra from all three runs at several representative frequencies. The three left-hand panels of this figure show spectra near the lower hybrid frequency with successively decreasing wavevector anisotropies with increasing e , same as reduced wavevector spectra in Fig. 5.8. The other three columns of panels show that the wavevector anisotropies continue to decrease as the frequency goes up. This trend as function of frequency is exactly the opposite to the trend as function of time in Fig. 5.8. It could be explained by the linear wave-particle damping as illustrated in Fig. 5.1 right column. As! increase, for the same!, the oblique waves are damped more strongly than the parallel propagating modes, especially for lower e cases, and therefore, the anisotropy ofk ? k k is reduced at higher frequencies. Fig. 5.1 also shows that for the same!, quasi-perpendicular and oblique propagating modes almost always have much bigger damping ratej j than parallel propagating modes, especially for e = 0.01 and 0.1. So if linear wave-particle interaction is the primary factor in determining the spectra in Fig. 5.15, the anisotropy would bejk k jjk ? j, rather than thejk ? jjk k j shown by our simulations. A more possible explanation is due to nonlinear wave-wave coupling processes. In sec- tion 1.1.2, we described that the nonlinear three-wave interaction could proceed to enhance 102 Figure 5.15: k k k ? reduced magnetic fluctuation energy spectra, log[jB(k k ;k y )j 2 =8n e k B T e;t=0 ], at four frequency bands as labeled for the three simulations with initial values of e as labeled. These spectra, which are as functions of frequencies, are computed over the same time intervals as those used in Fig. 5.3(a) 103 Parent mode 1 Parent mode 2 Daughter mode 3 Eqs. 1.7-1.8 sign: + !=j e j 0.1 0.1 0.2 k k e 0.2 -0.3 -0.1 k ? e 0.3 0.4 0.9 Eqs. 1.7-1.8 sign: !=j e j 0.1 0.075 0.025 k k e 0.3 0.4 -0.1 k ? e 0.5 -0.4 0.9 Table 5.2: A detailed listing of wave mode data from the e = 0:1 and e = 2:0 simula- tion that satisfies the wave-wave coupling Eqs. 1.7 and 1.8. All data are picked from the simulation as illustrated in the middle panel of Fig. 5.15. possible modes of higher or lower frequencies, shorter or longer wavelength if the condi- tions 1.7 and 1.8 are satisfied. A detailed analysis of the simulation results indicates that our whistler turbulences which covers over a fairly broad frequency and wavevector spec- trum can interact with each other leading to enhanced power in both the higher and lower frequency band wave modes of more oblique propagating directions than the two coupling modes. An example of this can be found, easily if not trivially, from Fig. 5.15. Table 5.2 shows the set of wave modes in the e = 0:1 and e = 2:0 simulation that interact and by using the positive and negative sign in Eqs. 1.7 and 1.8 produce the enhanced power in higher and lower frequency band more oblique propagating wave modes. A similar analy- sis for whistler simulations in the inner magnetosphere is illustrated in [82]. It is interesting to note that the parent modes seen in Table 5.2 propagate at different angles and produce a third mode which is more oblique than the two parent modes. This general trend such that wave modes become more oblique is supported by frequency mismatch test introduced in section 5.6. To conclude the frequency analysis, Fig. 5.16 quantitatively shows the wavevector anisotropy factor tan 2 B as a function of real frequency !. Consistent with 2D energy spectra in Fig. 5.15, the anisotropy decreases as e increases, and decreases as frequency 104 10 −2 10 −1 10 0 10 0 10 1 10 2 /| e | tan 2 B e = 0.01 e = 0.1 e = 1.0 Figure 5.16: Wavevector anisotropy factor of the cascaded magnetic fluctuation energy as a function of frequency during the three simulations with initial e = 0.01, 0.10, and 1.0. increases. However, after !=j e j 0:3, the larger e cases become more anisotropic than the lower e case. Besides, the high e = 1:0 case anisotropy factor stops decrease and even increase a little bit as frequency continue to increase, whereas the lower e cases anisotropy factor still decrease as frequency continue to increase. Here again, linear particle interaction is responsible for such anisotropy relations between three e cases. In Fig. 5.1 right column, at real frequency !=j e j 0:22 of the e = 1:0 case, the parallel propagating damping ratej k j have exceed the oblique propagating damping ratej oblique j. Therefore, linear damping here will reinforce the anisotropyjk ? j>jk k j instead of reducing it as the case in !=j e j 0:22 regime. In contrast, for e = 0.01 and 0.1, the ratio j k = oblique j 1 for the same!, and continues to decrease with increasing!. Thus the reduction of anisotropyjk ? j >jk k j continues as! increases. In other words,j k = oblique j is largest for e = 1:0 case, and lowest for e = 0:01 case at!=j e j 0:22. Therefore, in the higher frequency range we get anisotropy as a function of e in completely opposite order, as compared with those in lower frequency range. 105 Table 5.3: Computational parameter difference in the new run simulation Simulations = De ; = e ; = e 1:0= p 10, 0:1= p 10, 0.1 L x;y;z = De ;L x;y;z = e ;L x;y;z = e 162, 16.2, 51.2 ^ v te p 10 t! e 0.025 5.5 Electron-scale spectral break scaling As listed as one of objectives of this dissertation, we want to find some properties or quantities existing in the simulations that could effectively distinguish whistler and kinetic Alfv´ en modes in the short-wavelength turbulence. One possible approach is to investigate the detailed scaling of the electron-scale spectral break, to determine whether it scales as the inverse electron inertial length (jk ? e j' 1:0) or as the inverse electron gyroradius (jk ? e j' 1:0). The kinetic Alfv´ en wave interpretation has argued for the 1 e scaling [76], and 3D gyrokinetic simulations have demonstrated the break atjk ? e j = 1:0 [45]. How- ever, all these simulations and observations are at e ' 1:0, where e ' e . So there is no way to tell the detailed scaling from a e = 1:0 simulation or observation alone. In contrast, e = 0:1 simulations, as we have carried out, havek e = 3:16k e . Therefore, the condi- tionjk e j = 1:0 is located at 3.16 in a spectrum normalized tojk e j, while the condition jk e j = 1:0 is still at 1.0. This difference could be enough to distinguish the spectral break scaling and help identify the constituent modes of short wavelength plasma turbulence. Since the cascade process may be a function of the initial spectrum shape in the system, it is more appropriate to perform two separate simulations for comparison with initial same pattern spectrum normalized to same value of k e and k e . Thus, we renormalized the computational model to e in contrast to the previous normalization to e , and imposed the initial spectrum of same patterns, with wavevector normalized atk k=? e =0:1227, 0:2454, and0:3682, andk ? e = 0. We denote this simulation at e = 0:1 and e = 2:0 106 as the “new run”, and the corresponding previous simulation as the “old run”. The different computational parameters used in the new run are listed in Table 5.3. Figure 5.17 compares the scaling of electron-scale spectral break based on such idea. The left panel shows the magnetic fluctuation energy spectra normalized tojk ? e j from the early-time stage of the new run, with the blue dashed line indicating the power law fitting, and the vertical black dashed line marking the inverse electron inertial length (jk ? e j = 1:0). Clearly, there is no observable spectral break nearjk ? e j = 1, and the whole cascaded spectra could be described by a single power law function roughly/k 5:0 ? . The right panel of Fig. 5.17 illustrates the spectra from the corresponding old run sim- ulations at e = 0:1 and e = 2:0 where both computational grid spacing and initial spectra wavevector are normalized to e . The green solid line represent the old run from previous section 5.4. To better extend the spectrum down to shorter scale after the spectral break, we applied a spectral filter, as described in section 2.3.2, to reduce the noise by filtering the current densityJ source from PIC’s superparticles, and the parameters of the filter are labeled in the panel’s legend. The corresponding filtered result of the old run is shown in the black solid line. Besides, to better extend the spectrum up to the larger scale before the spectral break, we remove the old run’s initial fluctuation modes atk k=? e =0:2454, and 0:3682, only keep those atk k=? e =0:1227 andk ? e = 0. With fewer initial imposed modes, the decaying turbulence could be able to cascade from the largest scale possible, and form a much smoother spectrum at larger scale as opposed to our previous simula- tions’ spectra which usually oscillate at larger wavelength. The drawback of such single mode initial fluctuations is that it takes considerably longer time for the turbulence to form a quasi-steady forward-cascaded spectrum since the initial available modes for wave-wave coupling to proceed is very limited. The red solid line represent such a simulation. The blue dashed line indicating the power law fitting, and the vertical black dashed line marking the inverse electron gyroradius (jk ? e j = 1:0). Consistent with our previous results, there is a distinct spectral break atjk ? e j = 1:0, and there is no observable break around the 107 10 −1 10 0 10 1 10 −8 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 e =0.1, new runs E(k )/K e,t=0 |k e | 10 −1 10 0 10 1 =1.0 k e 10 −1 10 0 10 1 e =0.1, old runs |k e | old run, | e |t=134 old run w/ filtering J: a1=0.5, a2=20 old run w/ only largest initial modes new run, | e |t=22 new run, | e |t=67 k e =1 k −5.0 k e =1 k −5.0 k −2.8 Figure 5.17: k ? reduced magnetic fluctuation energy spectra at corresponding early time from the various simulations at e = 0:1 and e = 2:0. The blue dashed lines represent power law functions as labeled for comparison against the simulation results, and the vertical black dashed lines indicates the scaling as labeled. The left panel, normalized tojk ? e j, shows the spectra of the new run as defined before. The right panel, normalized tok ? e , shows the spectrum of the old run as defined before in the green sold line, the spectral filtered spectrum of the old run in the black solid line with the filter parameters as labeled, and the spectrum of the modified old run in the red sold line with only largest initial fluctuation modesk k=? e =0:1227 retained. 108 inverse electron gyroradius at the vertical dashed line. More importantly, the right panel’s spectral index =5:0 atjk e j > 1:0 is same as the left panel’s value. This indicates that the left panel is indeed a higher-resolution larger-range reproduction of the spectrum atjk ? e j& 3:16 from the right panel. The spectrum from the new run covers over a mag- nitude in wavenumber centered aroundjk ? e j = 1, and there are no spectral break around. Therefore, we are confident to conclude that, for whistler turbulence, the electron-scale spectral break happens atjk ? e j' 1:0, notjk ? e j' 1:0. Note the similarity to [12]’s observations that the inertial range proton-scale spectral break scales as the inverse pro- ton inertial lengthjk ? p j 1:0 but not the inverse proton gyroradius. We think it is an important analogy to our finding. An important implication from this finding is that if observations could focus on only those low e 1:0 time periods (although it is not typical), we could construct a energy spectrum with distinguishable 1 e and 1 e scaling. Assuming the kinetic Alfv´ en turbu- lence indeed has spectral break atjk ? e j' 1:0 as [76, 45] argue, we could tell whether the major identity of the short wavelength solar wind turbulence is whistlers or kinetic Alfv´ en waves by looking at the exact electron-scale spectral break location. Such a suggestion, however, would require considerable solar wind data analysis, and is beyond the purview of this dissertation. 5.6 Wave-wave coupling and the forward cascade In this section we develop a simple test [20] to examine whether wave-wave coupling implies a wavevector anisotropy for the forward cascade of whistler turbulence. This test has also been applied recently to lower-frequency turbulence consisting of Alfv´ en, magne- tosonic, and magnetosonic-whistler fluctuations [34]. Our test begins with the Eqs. 1.7 and 1.8 which represent conservation of fluctuation energy and momentum in nonlinear wave-wave coupling in plasmas [86]. Here we assume 109 !(k) is the real frequency part of the solution to the linear kinetic dispersion equation. Although the rates of energy transfer must be obtained from the full nonlinear dynamical equations of the plasma, Eqs. 1.7 and 1.8 represent necessary (but not sufficient) conditions for wave-wave coupling to proceed. In the simplest case where the dispersion relation is ! =v k (5.5) Equations 1.7 and 1.8 are trivially satisfied and turbulent cascade processes proceed directly for all frequencies and wavenumbers. The collisionless magnetized plasmas of space, how- ever, are the basis for many different normal modes, each of which exhibit anisotropic prop- agation (due to the presence of the background magnetic fieldB o ), strong dispersion (i.e., violation of Eq. 5.5), and wave-particle dissipation (e.g., Landau and/or cyclotron damp- ing) [63, 33]. Then the matching conditions are satisfied for a limited or perhaps empty range of! and k. We may gain insight into wave-wave cascade processes in plasma turbulence by com- bining linear dispersion analysis with Eqs. 1.7 and 1.8. It is difficult to satisfy both equa- tions exactly for most plasma normal modes, so we assume that wave dissipation and/or wave growth allow a small frequency mismatch ! in the satisfaction of Eq. 1.7 , i.e., ! 1 +! 2 =! 3 + ! (5.6) Furthermore, we assume the cascade process is one-dimensional in wavevector space so thatk = ^ ek where ^ e is a unit vector with fixed direction. Then Eq. 1.8 reduces to k 1 +k 2 =k 3 (5.7) 110 Mithaiwala et al. [58] examined the coupling between two relatively high-frequency whistlers and a low-frequency kinetic Alfv´ en wave. Our PIC simulations of whistler turbu- lence run for only a few proton cyclotron periods and do not provide a long-time represen- tation of kinetic Alfv´ en waves. So we here consider wave-wave coupling among three sep- arate whistler modes. Although whistler dispersion at relatively long wavelengths exhibits ! kk k dependance, e.g. Eq. 1.2, there is a wavenumber regime (0.20 < kc=! e < 1.0) on which whistler dispersion approximately satisfies the dimensionless equation ! =ak +b (5.8) where a and b are fitting parameters to the solution of the dispersion equation and are functions of e and, e.g. Fig. 5.1. Using Eq. 5.8 in the real part of Eq. 5.6 yields the dimensionless result a(k 1 +k 2 ) + 2b =ak 3 +b + ! !=j e j =b 0 (5.9) In other words, Eq. 1.7 is more likely to be satisfied if the quantityb 0 is as small as possible. Figure 5.18(a) illustrates whistler dispersion at k ? = 0 for e = 0.01, 0.10, and 1.0 obtained from the numerical solutions of the same dispersion solver described in section 5.3. The solid lines indicate the full dispersion equation solution for!(k k ), and the dashed lines indicate the damping rates. Although ! k 2 at k e 1, the dispersion approx- imately satisfies Eq. 5.8 at shorter wavelengths. The figure further suggests that the fre- quency mismatch roughly decreases with increasing e at parallel propagation modes, since the dispersion relation becomes more dispersive as e decreases. For the more general case of non-parallel propagation, the real frequencies of the whistler mode also approximately satisfy Eq. 5.8 at sufficiently large wavenumbers for a broad range of propagation direc- tions. Figure 5.18(b) shows whistler dispersion at = 75 o for the same three values of 111 Figure 5.18: Linear theory results for whistler dispersion properties as functions of wavenumber. (a) = 0 o . (b) = 75 o . Solid lines represent the real frequency ! and dotted lines show the damping rate . The blue, green and red lines correspond the three simulation sets in section 5.4 at e = 0.01, 0.1, and 1.0, respectively. e . Here, however, the dispersion relation becomes more dispersive and the frequency mismatch increases with increasing e . Finally, Fig. 5.19(a) and (b) show two fitting examples for Eq. 5.8 of quasi- perpendicular propagating whistler mode at e = 0.1 and 1.0. Clearly, e = 0:1 has much smaller frequency mismatch than that of e = 1:0 at quasi-perpendicular propa- gation. Panel (c) quantifies the frequency mismatch, showing !=j e j, calculated from 112 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.2 0.4 0.6 0.8 1 wh_be=0.10_80^o om_r/Om_ce y = -0.023072 + 0.12626x R= 0.99969 om_r/Om_ce kc/wpe ω/|Ω e | kλ e (a) Fitting at e = 0:1; = 80 o ω/|Ω e | kλ e 0 0.05 0.1 0.15 0.2 0.25 0 0.2 0.4 0.6 0.8 1 wh_be=1.0_75^o om_r/Om_ce y = -0.070516 + 0.28542x R= 0.99751 om_r/Om_ce kc/wpe (b) Fitting at e = 1:0; = 75 o 0 0.04 0.08 0 40 80 Figure 9 θ β e = 0.01, 0.10, 1.0 Δ ω |Ω e | (c) Frequency mismatch over the whole propagation range at the three e cases Figure 5.19: Linear theory results for the frequency mismatch ! for whistler fluctuations over 0.20kc=! e 1.0 as a function of propagation angle. The blue, green and red lines correspond the three simulation sets in section 5.4 at e = 0.01, 0.1, and 1.0, respectively. linear theory using the approach illustrated in panel (a) and (b), as a function of for each of the three simulation sets in section 5.4. At e = 0.01 and 0.10, the frequency mismatch is relatively large at quasi-parallel propagation, but decreases toward zero as approaches the perpendicular. This suggests that the forward cascade of whistler fluctuations is favored 113 at quasi-perpendicular propagation, consistent with our simulation results in Fig. 5.5 of a strongk ? >>k k wavevector anisotropy at e << 1. For parameters corresponding to the e = 1.0 simulation, linear dispersion of whistlers at> 75 o does not well satisfy Eq. 5.8, and at smaller angles of propagation the frequency mismatch is relatively large compared to those of lower e cases. This suggests there is less difference between the forward cas- cade rates at quasi-parallel and quasi-perpendicular propagation, and is therefore consistent with our simulation results of much weaker wavevector anisotropies at e = 1:0. We fur- ther note a possible source of the peculiar “star-shape” spectra shown at late times of e = 0.10 and 1.00 in Fig. 5.8. Panel (c) shows that, at the two larger values of e , the frequency mismatch has its largest values at intermediate values of, suggesting that the wave-wave interactions which lead to the forward cascade are weakest for such oblique wavevector directions. 5.7 Conclusions We have carried out 3D particle-in-cell simulations of an initial value problem for whistler turbulence in a homogeneous, magnetized, collisionless plasma with three initial values of e , namely, 0.01, 0.10 and 1.0. The e is varied by changes in the background magnetic field magnitudeB o , whereas the initial electron temperature is the same for all three runs. The total energy in the initial magnetic fluctuations is also held constant, so the conse- quences of electron wave-particle heating can be compared directly. The simulations begin with an ensemble of relatively long wavelength whistler modes and follow the temporal evolution of the fluctuations. At early times wave-wave interac- tions dominate and yield a forward cascade into a broadband, turbulent spectrum at shorter wavelengths with a wavevector anisotropy in the sense of k ? > k k . In addition, wave- particle interactions yield fluctuating field damping and electron heating with a temperature anisotropy in the sense ofT k > T ? . At sufficiently late times these dissipative processes 114 dominate. Larger values of e correspond to a faster forward cascade in wavenumber and to a faster rate of electron heating, as well as to a less anisotropic wavevector distribution and to a less anisotropic electron velocity distribution. The overall electron heating yields T k > T ? , which implies that the primary linear wave-particle interaction is through the Landau resonance (although the cyclotron reso- nance which should yield a T ? > T k plays a larger role as e increases). The Landau resonance for whistlers is ineffective atkB o = 0, so Landau damping alone should lead to whistler spectra which are anisotropic in the sense of k k > k ? . But Figs. 5.5 and 5.16 show the simulated spectra have the opposite wavevector anisotropy. Therefore, lin- ear wave-particle interaction must play a minor role in shaping the wavevector spectra, whereas nonlinear wave-wave interactions are overall stronger and faster processes, and ultimately determine the wavevector anisotropy. Frequency analysis and simulation dispersions indicate a transfer of wave energy to higher frequencies at early times, and to lower frequencies throughout the simulations. At late times, enhancement near the lower hybrid frequency lh ' 0:024j e j is also observed for e = 0:01; 0:1 cases. The wavevector is very anisotropic at lower frequencies whereas relatively isotropic at higher frequencies. Bearing in mind the two-phase characteristics illustrated in Fig. 5.3(b) and 5.5, we could summarize as follows. At early times, the nonlinear wave-wave coupling pumps fluctuation energy toward shorter wavelengths, anisotropick ? >k k wavevectors and both higher and lower frequency regime. The dominant force in determining the spectra at early stage is the wave-wave cou- pling. For same frequency band, linear wave-particle damping has much stronger damping of the oblique and quasi-perpendicular modes than the parallel propagating modes, espe- cially for lower e cases. The two processes collectively determine the energy!k spectra such that the anisotropy decreases at higher frequency range shown in 3D fluctuation dis- persion of Fig. 5.14. 115 At late times, wave-wave coupling gradually weakens due to the reduced amplitude of the decaying turbulence, and wave-particle interactions dominate. Linear Landau damp- ing has the strongest damping of whistlers at oblique propagation (40 o < < 70 o ). The damping of such waves supports the continuous growth in wavevector anisotropy, and con- tributes to the peculiar “star-shape” spectra shown at late times of e = 0.10 and 1.00 in Fig. 5.8. Because of the presence of the enhanced lower hybrid frequency fluctuations at the late-time stage in Fig. 5.12, nonlinear Landau damping may also contribute to the total wave-particle dissipation in our simulations. The k ? reduced energy spectra demonstrate the characteristic spectral break near the inverse electron inertial length, with asymptotically constant values of spectral index at nonlinear-process dominated large fluctuation amplitude solar wind plasmas. This indi- cates again the importance of nonlinear processes in determining the final state of the power spectrum. Although the e = 1:0 simulation wavevector spectral indices could not match closely the frequency spectral indices from observations, we interpret this difference due to the questionable applicability of Taylor’s hypothesis on whistler waves at e = 1:0 sim- ulations. In contrast, the energy frequency spectra from our simulations exhibit a close similarity with solar wind observations frequency spectra. Besides, the power spectra of electron density fluctuations share a close similarity to the solar wind observation and MHD predications, both qualitatively and quantitatively. These properties reinforce the interpre- tation that whistler fluctuations could be the substantial constituent of solar wind turbulence at higher frequencies and short wavelengths. Further investigation on the scaling of electron-scale spectral break position suggests that, for whistler turbulence, the spectral break happens atk ? e ' 1:0, instead ofk ? e ' 1:0. If observations could obtain spectra from low e 1:0 time periods, and if kinetic Alfv´ en turbulence indeed has spectral break atk ? e ' 1:0 as argued by [76, 45], we could use this distinct property to identify the constituent fluctuations of solar wind turbulence at electron scale as predominantly whistler or kinetic Alfv´ en waves. 116 We have used the frequency mismatch based on linear theory results for whistler fluc- tuations to determine the conditions favoring a forward cascade of whistler turbulence on the range 0.20kc=! e 1.0. At e = 0.01 and 0.10, the matching conditions imply that the forward cascade rate at quasi-perpendicular propagation should be much faster than at quasi-parallel propagation. At e = 1.0, the matching conditions imply that the forward cascade rate for both quasi-parallel and quasi-perpendicular propagating may be similar. Both these findings are consistent with our simulation results. Our interpretations of the simulation results, so far, have been primarily framed in terms of the quasilinear scenario. There are other possible strong nonlinear processes, e.g., dis- sipation at reconnection sites, which may also play a role here. The final question left is how much contribution to the wave energy dissipation and electron heating from all these processes, and whether the intermittent current sheets play a role in the collisionless, homo- geneous whistler turbulence dissipation process. These questions are addressed in the next chapter 6. 117 CHAPTER 6: ENERGY DISSIPATION BY WHISTLER TURBULENCE FORWARD CASCADE In this chapter, we focus on the question of fluctuation energy dissipation and electron heat- ing in the homogeneous whistler turbulence. Using data from particle-in-cell simulations described in previous chapters, we compare the total field fluctuation dissipation rate versus the dissipation due to linear collisionless damping as functions of e and e , and make a statistical diagnosis of turbulent intermittency, e.g., the current sheet. With an introduction in section 6.1, diagnostic quantities are described in section 6.2. In section 6.3, detailed simulation results and diagnoses are presented and discussed, followed by conclusions in 6.4. 6.1 Introduction An unsolved problem in plasma turbulence is how fluctuation energy is dissipated at short- wavelength scales. One way to express this question is to ask what are the relative con- tributions to fluctuation energy dissipation and plasma heating from linear and nonlinear processes. The dissipation of turbulence in neutral fluids is usually described in terms of a forward cascade to short wavelengths where collisional processes such as resistivity or viscosity convert fluctuating fluid energy into thermal energy. Turbulence in collision- less magnetized plasmas is also often described in terms of a forward cascade. However, particle collisions are too infrequent in collisionless plasmas to provide the necessary dis- sipation at short wavelengths. The conversion of fluctuating electric and magnetic field energy into plasma kinetic energy may be described in terms of the damping of linear 118 plasma modes (i.e. via linear Landau and/or cyclotron damping), by various higher order wave-particle dissipation processes (e.g., nonlinear Landau damping), and by reconnec- tion at current sheet structures in the plasma. Note that, as described in section 1.1.2, the first two kinds of processes are within the framework of the quasilinear scenario where the plasma is considered as a relatively homogeneous medium, whereas the third kind is the core of the intermittency scenario, which describes a turbulent plasma as an intrinsically inhomogeneous and intermittent medium consisting of fundamentally nonlinear, coherent structures such as current sheets. To identify the relationships between the quasilinear scenario and the intermittency scenario in describing the dissipation of whistler turbulence, we here apply particle-in-cell (PIC) simulations of whistler turbulence at variable e and e to address the question of electron heating and dissipation of the fluctuating fields. We note that, although several authors [97, 46, 93] have done similar studies using strong wave-driven turbulence kinetic simulations, our simulations are decaying turbulence cascades with imposed initial small perturbations. By doing this, we ensure that our simulations run on a relatively homogenous plasma medium, consistent with our goal to examine the very fundamental properties of the whistler turbulence. 6.2 Initial conditions We follow the simulations described in chapter 3, 4, 5 and begin with a collisionless, homo- geneous, magnetized plasma upon which is imposed an initial ensemble of relatively long whistler fluctuations. At early times, the whistler modes undergo a forward cascade to a broadband, anisotropic, turbulent spectrum at shorter wavelengths. At later times the field fluctuations are dissipated and the plasma electrons are heated, with a temperature anisotropyT k >T ? . 119 To study dissipation of whistler turbulence, we organized our previous simulation runs into two sets, one varying with electron initial e = 0:01; 0:1; 1:0 while fixing the initial fluctuation amplitude e = 2:0 as in chapter 5, the other varying the initial fluctuation amplitude e = 0:5; 1:0; 2:0; 5:0 while fixing electron initial e = 0:1 as in chapter 4. In certain diagnoses, we also include data from variable e at three different e . All plasma and computational parameters are same as those in corresponding previous chapters. We define the total dissipation rate t from the simulations by computing the change rate of fluctuating magnetic as well as electric field energy densities as follows, t = d[B 2 (t) +E 2 (t)]=8 dt (6.1) The instantaneous t can be calculated from the slope of the curve of total fluctuating energy time history. A small smoothing factor is used to calculate t , i.e., average over every 6 electron plasma periods to reduce the noise. We also define the linear damping rate of the fluctuating fields as l . As described in appendix A.1, the linear theory assumes the fluctuation amplitude is sufficiently small that the perturbation could be described as, for example,Be t , and the corresponding fluctuation energyB 2 (t)=8 is thenB 2 e 2 t =8. Therefore, using Eq. 6.1, the instantaneous l can be approximated by an integral over all wavevectorsk of the quantity as below, l = Z k 2 (k)[B 2 (k) +E 2 (k)]=8dk (6.2) where is the imaginary part of the wave frequency computed from linear kinetic dis- persion theory [33] based on plasma conditions at t = 0. The integral is computed as k 2 (k)(B 2 (k) +E 2 (k))=8 on the computational grid in the simulations. Both of these damping rates are normalized to the initial electron kinetic energy densityK e;t=0 in units of inverse time. 120 The fluctuating magnetic field increment as a function of spatial separation lengthr is denoted as r j B i , where r j B i = B i (j +r)B i (j), and i;j could be either x;y, or z direction. Kurtosis of magnetic field increments is then defined as<= 4 ij >: <= 4 ij >= < r j B 4 i > 4 (6.3) where =<j r j B i j 2 > 1=2 is the non-centralized standard deviation and< > denotes an ensemble average over three dimensional space. The correlation coefficient of two quantitiesX andY is denoted as (X;Y ) = cov(X;Y ) X Y (6.4) wherecov denotes covariance, and is the standard deviation. The correlation coefficient is a normalized value between -1 and 1. Value close to 1 indicates that there is a positive linear relationship between the data, while values close to -1 indicate a negative linear rela- tionship. Value close to or equal to 0 suggests there is no linear relationship between the data. Note that the correlation reflects the non-linearity and direction of a linear relation- ship, but not the slope of that relationship, nor many aspects of nonlinear relationships. 6.3 Particle-in-cell simulations and discussion First, we quantify the amount of fluctuating field energy dissipation is contributed from lin- ear damping and nonlinear dissipation. Figure 6.1 shows the comparison between the total damping rate t in solid lines versus the linear damping rate l in dashed lines. The differ- ence between the the total and linear damping rates indicates the dissipation due to non- linear dissipation mechanisms which may come from various higher order wave-particle dissipation and/or dissipation associated with reconnection at intermittent current sheet structures. We note that in early time of some cases, the approximated linear damping 121 10 −5 10 −4 10 −3 10 −2 Damping Rate (a) time history of total (solid) and linear (dashed) damping rate as function of e at e =2.0 10 −1 10 0 10 1 10 2 10 3 10 −4 10 −3 10 −2 Damping Rate | e |t (b) time history of total (solid) and linear (dashed) damping rate as function of e at e =0.1 10 −1 10 0 10 1 10 −4 10 −3 10 −2 Damping Rate e (c) total (solid) and linear (dashed) damping rate as function of e and e at early times e =0.01, t e =0.01, l e =0.1, t e =0.1, l e =1.0, t e =1.0, t e =0.5, t e =0.5, l e =1.0, t e =1.0, l e =2.0, t e =2.0, l e =5.0, t e =5.0, l e =0.01, t e =0.01, l e =0.1, t e =0.1, l e =1.0, t e =1.0, l Figure 6.1: Comparison of total damping rate t (solid lines) and linear damping rate l (dashed lines). (a) as functions of e andj e jt at e = 2:0; (b) as functions of e and j e jt at e = 0:1; (c) as functions of e and e at their respective early times, that is for e = 0:01 cases,j e jt = 707:1, for e = 0:1 cases,j e jt = 111:8, and for e = 1:0 cases, j e jt = 35:4. 122 rate l is larger than the measured total dissipation rate t . The reason of this unphysical phenomenon is due to the approximation of linear dissipation rate in Eq. 6.2. After all, the fluctuation fields go as/ e t only if everything is linear, and in early time, the domi- nant process is the nonlinear wave-wave interaction, and the fields do not necessarily go as /e t . Thus, the estimation of the linear damping rate via l is more valid when dissipation processes dominate the simulation. Figure 6.1(a) illustrates the time history of t and l as a function of e at e = 2:0. As e decreases, the difference between t and l increases, especially in early time when the fluctuation amplitude is large. The difference between t and l is approximately the nonlinear damping rate. It indicates that in e = 0:01 whistler turbulence where linear dissipation is weak, most of fluctuation damping comes from the nonlinear dissipation, since the fluctuation could sustain for much longer times and much larger amplitudes to facilitate the nonlinear processes; whereas in higher e = 1:0 where the linear dissipation mechanism is much stronger, linear damping rate is approximately same as the the total dissipation rate. Figure 6.1(b) shows the time history of t and l as a function of e at e = 0:1. Not like Fig. 6.1(a), the ratio of l = t remains roughly the same as e increases until a larger decrease when e reaches 5.0. This indicates that for warm and hot plasmas, the nonlinear dissipation mechanism involved in these whistler turbulence simulations is a relatively weak function of e , however, a stronger function of e . Figure 6.1(c) compares t and l as functions of both e and e at their relative early times. One could notice that, for the same e , as e increases, the linear dissipation rate increases slower than a linear function of e , indicating a breakdown of linear damping assumption. This becomes more prominent when e reaches 5.0. The difference between t and l rate becomes maximized at largest e for the same e , indicating the strongest nonlinear dissipation at e = 0.1 and 1.0 only appear when amplitudes are sufficiently large. For example, at e = 5:0, for e = 0:01, l = t ' 2%, t l = 0:00267; for e = 0:1, 123 l = t ' 50%, t l = 0:00274; for e = 1:0, l = t ' 77%, t l = 0:00369. This is consistent with the findings in Figs. 4.1(c) and 5.4 that at larger e 5:0 fully nonlinear processes gradually play a major role in energy dissipation. By choosing a specific time, it is also clearer to show that the higher e is, the less the difference between total and linear damping is. We next investigate whether intermittent current sheet structures are formed in our sim- ulations. A strong current sheet formation usually indicates a strong dissipation associ- ated with reconnection. In a collisionless plasma, magnetic reconnection in the presence of a guide magnetic field is expected to develop an electron diffusion region at the scale k ? e & 1 [8], where a strong current sheet is expected to form. In Figures 6.2 - 6.4, we plot the parallel current densityJ k in ax-y plane for both e = 2.0 and 5.0 cases at e = 0.01, 0.1, and 1.0, respectively. These contour plots are chosen at each simulation’s early time stage as labeled in their respective captions. To explore the contribution to the current from different scales, we applied spatially band-pass filterers, that is, in Figs. 6.2 - 6.4, the left columns contain J k withinjk ? e j < 0:49, the middle columns containJ k within 0:49jk ? e j< 1:51, and the right columns containJ k within 1:96jk ? e j. The middle columns of Figs. 6.2 - 6.4, whose filter are approximately in the range ofk ? e 1, shows that our 3D simulation produces current sheets at the electron inertial length scale, consistent with such development in wave-driven plasma turbulence simulations [97, 46, 93]. The relatively lack of significant current density in right column of Fig. 6.4 shows that current sheets do not form strongly at scalesjk ? e j> 2. Figures. 6.2 through 6.4, show that, in general, for the same e , the current sheet inten- sity (color of the contours) becomes stronger with increasing e . It is most obvious for lower e = 0:01 cases, and is least obvious for higher e = 1:0 cases. This trend is consis- tent with that of Fig. 6.1(b) and (c). On the other hand, for the same e , with increasing e , current sheet intensity shifts from smaller scale structure to larger scale structures. Becasue 124 Figure 6.2: e = 0:01 cases: Parallel current density, J k , for a perpendicular plane (xy) at timej e jt = 707:1, with initial value of e and different band-pass filters applied as labeled. The color contours represent simulation values ofJ k = p K e;t=0 . 125 Figure 6.3: e = 0:1 cases: Parallel current density,J k , for a perpendicular plane (xy) at timej e jt = 134:2, with initial value of e and different band-pass filters applied as labeled. The color contours represent simulation values ofJ k = p K e;t=0 . 126 Figure 6.4: e = 1:0 cases: Parallel current density,J k , for a perpendicular plane (xy) at timej e jt = 35:4, with initial value of e and different band-pass filters applied as labeled. The color contours represent simulation values ofJ k = p K e;t=0 . 127 strong current sheet is expected to form in the range ofk ? e 1 [8] and e = e = p e , for lower e cases, stronger current sheet structure should appears in the rangejk ? e j> 1. It is important to note that in our simulations the current sheets form self-consistently from the cascade of decaying turbulence without any driving force, and are not seeded or otherwise initialized. Therefore, we conclude that electron-scale current sheets form as a natural consequence of whistler turbulence, even in small perturbed, relatively homoge- neous plasmas such as these 3D PIC simulations. (n e T e ;J 2 ) e = 0:01 e = 0:1 e = 1:0 e = 1:0 0.138 0.245 0.326 e = 2:0 0.219 0.304 0.336 e = 5:0 0.247 0.368 0.412 Table 6.1: Correlation coefficients betweenn e T e andJ 2 . All grid points in the 3D simula- tion space are taken into account. The largest coefficients are listed here. To investigate the relationship between energy dissipation associated with current sheets and the total electron heating measured in the simulation, we compute the corre- lation coefficient, (n e T e ;J 2 ), between the current density squared J 2 and the measured electron temperaturen e T e in the simulations as functions of both e and e . The values are listed in Table 6.1. Note that, in general, (n e T e ;J 2 ) does not reflect Landau / cyclotron damping’s correlation with electron temperature. The current density J on its own does not necessarily correspond to Landau damping, because in the micro-picture, J means a velocity shift of all the electrons relative to the ions. The correlation Table 6.1 shows that an increase in the initial wave amplitude trans- lates into a monotonic increase in the coefficient(n e T e ;J 2 ). The increasing correlation between current density and electron temperature with increasing e is consistent with the trend that nonlinear processes involved in our whistler turbulence play increasing roles in energy dissipation with increasing e as shown in Fig. 6.1(c). Besides, the correlation also increases with increasing e . In the limit of zero e , the magnetic field is rigid and there 128 are no current sheets. As e increases, the magnetic field becomes weaker and the plasma currents act to push the magnetic field around, forming current sheets. Therefore, current sheet formation and current sheet dissipation should probably increase with increasing e , at least for e . 1. The relative contribution to electron heating from linear damping could be incorporated into correlations by calculating the coefficient(n e T e ;E 2 ). For Landau / cyclotron damp- ing to be involved, there must be a microscopic quantity involved, such asE z orE. We note that(n e T e ;E 2 ) does not uniquely represent the correlation between linear damping and electron heating, but just involves the linear dissipation effects. Figure 6.5 shows the spatial contours of J 2 , n e T e and E 2 from two typical cases, that is, the upper panels represent e = 0:1 and e = 2:0 case at timej e jt = 134:2, and the lower panels represent e = 1:0 and e = 5:0 case at timej e jt = 14:1. The correlation coefficient ratio (neTe;J 2 ) (neTe;E 2 ) from e = 0:1 case is 4.83, whereas (neTe;J 2 ) (neTe;E 2 ) from e = 1:0 case is 1.92. Since the overall total damping rate increases with e as shown in Fig. 6.1(a), this decreasing ratio with increasing e indicates that the relatively less importance of current sheet dissipation mechanism and more importance of linear damping mechanism in higher e plasmas, and this trend is consistent with that of Fig. 6.1(c) In general, as shown in Table 6.1, the current sheet dissipation processes involved in our whistler turbulence play increasing roles in energy dissipation with increasing e . However, the higher ratio e = 0:1 case in Fig. 6.5 has a smaller e = 2:0, whereas the lower ratio e = 1:0 case has a bigger e = 5:0. Therefore, the decreasing ratio (neTe;J 2 ) (neTe;E 2 ) with increasing e and e indicates that the relative contribution between current sheet dissipation and linear damping is ultimately determined by e rather than e . Nevertheless, the relatively strong correlation betweenJ 2 andn e T e indicates that current sheet formation is playing an important role in whistler turbulence simulations here. The spatial structures of J 2 resemble the electron temperature n e T e structures shown in Fig. 6.5, providing a qualitative evidence of dissipation in regions of coherent intermittent structures. 129 Figure 6.5: Current density square J 2 , electron temperature n e T e , and electric field energy E 2 from a perpendicular plane (xy). The color contours represent normalized simulation values of corresponding quantities. The left column isJ 2 , the middle column isn e T e , and the left column isE 2 . The upper panels are from e = 0:1 and e = 2:0 case at timej e jt = 134:2, and the correlation coefficient ratio (neTe;J 2 ) (neTe;E 2 ) = 4:83. The lower panels are from e = 1:0 and e = 5:0 case at timej e jt = 14:1, and the correlation coefficient ratio (neTe;J 2 ) (neTe;E 2 ) = 1:92. Note that the correlation coefficients are computed from all grid points in the 3D space, not just the 2D planes shown here. 130 −6 −4 −2 0 2 4 6 10 −8 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 r y B x / PDF Guassian r y =0.1 e r y =1 e r y =2.5 e r y =5 e r y =10 e r y =25 e (a) PDF( ry B x ) −6 −4 −2 0 2 4 6 10 −8 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 r y B z / PDF Guassian r y =0.1 e r y =1 e r y =2.5 e r y =5 e r y =10 e r y =25 e (b) PDF( ry B z ) Figure 6.6: Normalized PDF of fluctuating magnetic field increments from e = 0:1 and e = 2:0 case at timej e jt = 134:2 for different spatial separation lengthr as labeled. (a) the perpendicular component ry B x ; (b) the parallel component ry B z . Here the PDFs are normalized to the standard deviation of their respective magnetic field increments. The black circled lines indicate the Gaussian distribution. Quantitative measurements of intermittency can be framed in terms of probability den- sity function (PDF) of the fluctuating magnetic field increments r j B i as a function of spatial separation length r. Figure 6.6 shows the normalized PDFs for r y = 0:1 e , 1 e , 2:5 e , 5 e , 10 e , 25 e as labeled from e = 0:1 and e = 2:0 case at timej e jt = 134:2. Panel (a) illustrates one of the perpendicular magnetic field increments, i.e., ry B x ; Panel (b) illustrates the parallel magnetic field increments, i.e., ry B z . As the spatial lag r is decreased from large scales to smaller scales, both PDFs are found to become increas- ingly non-Gaussian, especially near the tail ends of PDFs. However, there is a noticeable anisotropic difference between the PDFs of the perpendicular and parallel field increments. PDFs of ry B z has much larger non-Gaussian tails than that of ry B x , indicating stronger nonlinearity. This anisotropy could be explained by the fact that whistler turbulence cas- cade is strong oblique propagation, which results in more contribution fromB k thanB ? in the probability density function calculation. In section 1.2.2, we note that recent solar wind observations [47] find that the PDFs of r B ? and r B k are similar for kinetic Alfv´ en 131 waves. Therefore, this anisotropic difference between whistlers from kinetic kinetic Alfv´ en waves might help in identifying the characteristics of short wavelength solar wind turbu- lence. The PDFs in Fig. 6.6 show enhanced tails at large amplitudes, but do not show more strongly peaked than Gaussians at small amplitude, as appeared in Fig. 1.8 (b) [47, 46]. We believe that the enhanced PDF tails are due to the intermittency in turbulence, but that the steepening of the peak is due to the large-scale inhomogeneities. When sum- ming over Gaussians with somewhat different temperatures, a so-called Castaing distri- bution is formed, which usually has a more peaked distribution at small amplitudes [89]. This implies that observations [47] and strong wave-driven simulations [46] may sum over many different small-scale measurements in inhomogeneous plasmas, so that their PDFs represent several different local Gaussians which result in the more peaked distributions at small amplitudes, whereas Fig. 6.6’s relatively Gaussian distributions at small amplitudes confirms that our whistler simulations run on a homogeneous medium. Another quantity as an indicator of intermittency and probably also current sheets is kurtosis as defined early in Eq. 6.3. A Gaussian distribution yields a kurtosis of 3; depar- tures from a Gaussian distribution correspond to increasing values of the kurtosis and are a measure of what has been called turbulent intermittency. Figure 6.7 shows the scale dependent kurtosis of the magnetic field increments as func- tions ofr and e for e = 0.01, 0.1 and 1.0, in left, middle, and right panels, respectively. In all three panels, the kurtosis grows from the Gaussian value of 3 or mildly non-Gaussian values of 3 4 at scales larger than 10 e to more non-Gaussian values near the electron scales e . The values of kurtosis increases with increasing e and e , and also gener- ally increase with decreasing spatial separation lengthr (increasing wavenumberk) until saturation at much smaller scales. This behavior of the increments clearly signifies the phenomenon of statistical intermittency of the magnetic field, and indicates that the inter- mittency in the whistler turbulence generally increases with e and e , the same trend as 132 10 −1 10 0 10 1 3 3.5 4 4.5 5 5.5 Kurtosis of r x B z r x / e e =0.01, | e |t=707.1 10 −1 10 0 10 1 Kurtosis of r z B y r z / e e =0.1, | e |t=134.12 10 −1 10 0 10 1 Kurtosis of r y B z r y / e e =1.0, | e |t=35.4 e =1.0 e =2.0 e =5.0 e =0.5 e =1.0 e =2.0 e =5.0 e =1.0 e =2.0 e =5.0 Figure 6.7: Kurtosis of the magnetic field increments r j B i as a function of spatial separation length r for different e and e as labeled. The left panel is from e = 0:01 cases at timej e jt = 707:1; The middle panel is from e = 0:1 cases at time j e jt = 134:12; The right panel is from e = 1:0 cases at timej e jt = 35:4. Data are low pass filtered atk x=y=z e = 5 to remove the statistical noise due to finite superparticle number of the PIC simulations. 133 the current and electron temperature correlation coefficient (n e T e ;J 2 ), shown in Table 6.1. The noticeable jump of the kurtosis values near k e between e = 5:0 cases and other lower e cases at e = 0.01 and 0.1 indicates the point when intermittent current sheet structures become relatively dominant in these whistler turbulence simulations. Note that the values of kurtosis of magnetic field increments in our simulations are much smaller than that in the recent large-scale 2D PIC simulation of sheared flow driven turbulence [97] which yields values up to 10:5. This indicates again that our whistler simulations remain relatively homogeneous. 6.4 Conclusions Using results from our 3D PIC simulations from chapters 3, 4 and 5, we have examined energy dissipation and electron heating properties of the forward cascade of homogeneous whistler turbulence. By comparing the total dissipation rate and linear damping rate as functions of e and e , we obtain physical insight into the relative importance of linear and nonlinear pro- cesses for dissipation of whistler turbulence in a collisionless, homogeneous, magnetized magnetized plasma. We find that linear damping constitutes a very small part of the total dissipation at e = 0.01, but that increasing e corresponds to an increasing fraction of linear damping. At high e = 1:0, linear damping contributes more than 75% of the total dissipation, depending on different e . On the other hand, increasing e reduces the relative contribution of linear damping to the overall dissipation, and this trend is more obvious at lower e cases. Our demonstration that linear damping represents an intermediate fraction of the total dissipation of whistler turbulence over a wide range of e and e stands in con- trast to the results of both [46], who found current sheet heating to be several orders of magnitude more efficient than damping by kinetic Alfven waves, and [93], who concluded 134 that collisionless damping of kinetic Alfven waves via the electron Landau resonance is sufficient to account for all of their measured heating. We further quantify the intermittency in our simulations, by showing the PDFs of the magnetic field increments depart more strongly from Gaussian distributions as the sepa- ration length decreases. An anisotropy is found between PDFs of parallel and perpendic- ular magnetic field increments, which is probably due to the strong oblique propagation of whistler modes. The PDFs not showing strong peaks at smaller amplitudes confirms that our simulation maintains a homogeneous plasma environment. The kurtosis of the magnetic field increments, with decreasing separation length, similarly increases from the value of three, corresponding to a Gaussian PDF, at much larger than electron scales, to non-Gaussian values ranging from 3:5 5:5 near electron scales. It suggests that the intermittency in the whistler turbulence generally increases with increasing e and e . Although our whistler simulations are decaying turbulence cascades, as opposed to stronger driven turbulence, electron-scale current sheets still form clearly in our simula- tions. This suggests that electron scale current sheet generation may be a generic feature of collisionless turbulence. The correlation coefficients(n e T e ;J 2 ) between electron tem- perature and current density show a monotonic increase with e and e , implying that the amount of current sheet dissipation should increase with increasing e and e . This trend is the same as that of the values of the computed kurtosis, and we believe that there is a strong link between the current sheet dissipation and intermittency in the whistler turbulence. We estimate the relative contribution to the overall electron heating from current sheet dissipation and linear damping by computing the correlation coefficient ratio (neTe;J 2 ) (neTe;E 2 ) . The ratio decreases with increasing e , implying a less importance of current sheet dissi- pation and a more importance of linear damping in higher e plasmas. Because the relative contribution from linear damping generally increases with e and decreases with e , the decreasing ratio (neTe;J 2 ) (neTe;E 2 ) with increasing e and e also suggests that the relative contri- bution between current sheet dissipation and linear damping is ultimately determined by 135 e rather than e . These properties maintain a consistent correlation with the trends of the relative contributions from linear and nonlinear dissipation described earlier in this section illustrated in Fig. 6.1. Therefore, we conclude that the nonlinear dissipation processes involved in our whistler turbulence simulation is primarily associated with dissipation in regions of intermittent current sheet structures. 136 CHAPTER 7: THREE-DIMENSIONAL PARTICLE-IN-CELL SIMULATIONS OF WHISTLER TURBULENCE: FORWARD CASCADE VS. INVERSE CASCADE In this chapter, we present the results of the first fully three-dimensional particle-in-cell simulations of whistler turbulence in a magnetized, homogeneous, collisionless plasma in which both forward cascades to shorter wavelengths, and inverse cascades to longer wavelengths are allowed to proceed. The simulation setup and corresponding plasma and computational parameters are given in section 7.2. The detailed simulation results are discussed and a comparison with previous forward cascade simulations are presented in section 7.3. Conclusions are given in section 7.4. 7.1 Introduction In previous chapters 3, 4, 5 and 6, we focused on the forward cascade of whistler tur- bulence where the wavelength of the initial fluctuations imposed upon the system are of the same order as the simulation dimensions. Hence, in those simulations, only the forward cascade to shorter wavelengths could proceed, and the inverse cascade pro- cesses are excluded. This approach is also adopted in most recent kinetic simulations of short-wavelength solar wind turbulence (both driven and decaying turbulence simulations) [45, 94, 97, 46, 37, 79, 80, 15, 78]. One reason for using this configuration is because we assume that the Kolmogorov picture is valid and that longer wavelength turbulence can carry fluctuations all the way down to shorter electron scales via forward cascades. 137 The “uphill” spectrum longer than the largest wavelength of our system should also form a broadband ensemble of incoherent field fluctuations like those in our system, and the presence of such a spectrum in homogeneous plasmas should not affect the results of the whistler turbulence evolution within our smaller simulation domain. Under this assump- tion, we neglect the dynamics on larger scales, and only focus on shorter scales compa- rable to whistlers. However, such a assumption also implies that the dynamics of inverse cascade to longer scale, if any, will not affect the conclusions drawn from the whistler turbulence forward cascade evolution, and thus the exclusion of inverse cascade process from our system could be justified. Another reason is due to computational constrains. The computational resources required for kinetic simulations of short-wavelength turbulence is approaching the limits of state-of-the-art supercomputers. Hence, in most cases, one sim- ply cannot afford to use a simulation domain much larger than the characteristic length of the short-wavelength solar wind turbulence. In this chapter, we extend our previous simulations by including the inverse cascade process in the simulation. We choose a simulation box larger than that used in our earlier simulations, and let the inverse cascade process freely compete against the forward cas- cade process. This fully 3D large scale PIC computation is more general than previous inverse cascade computational models [102, 32, 24] and permits a direct comparison of the competition between forward and inverse cascades for whistler turbulence. Thus, we inves- tigate the preferential evolution direction, energy cascade and dissipation rate, wavevector anisotropy of whistler turbulence and corresponding electron heating associated with the inverse cascade, and comparisons against previous corresponding forward cascade case are presented. 138 Table 7.1: Computational parameter difference in Run 1 and Run 2 simulations L x;y;z = De ;L x;y;z = e total superparticles initial spectrumk Run 1 512, 51.2 1:72 10 10 jk k=? e j = 0:1227; 0:2454; 0:3682,k ? e = 0 Run 2 2048, 204.8 5:50 10 11 0:2454jk k=? e j 0:4909 7.2 Initial conditions and simulation setup The plasma is collisionless and homogeneous with periodic boundary conditions in all three directions. We define two runs for this chapter. Run 1 is the original forward cascade case at e = 0:1 and e = 2:0. All plasma and computational parameters are same as before. The simulation domain is stillL x;y;z = e = 51:2 with a grid spacing of = 0:1 e , and the initial fluctuation spectrum att = 0 is still chosen to consist of whistler waves at k k=? e =0:1227,0:2454, and0:3682, andk ? e = 0 as before. Thus, Run 1 only allows forward cascade to proceed. Run 2 is the new large scale run which allows both forward and inverse cascades. Most of the plasma and computational parameters of Run 2 are same as Run 1. The differences are summarized in Table 7.1. Run 2 uses a larger simulation box, that is 204820482048 cells. With a grid spacing still = 0:1 e , the simulation domain is enlarged toL x;y;z = e = 204:8. For Run 2, we choose a different set of initial fluctuation modes. Our computational power is limited, and a practical maximum simulation box size 204:8 3 e is able to extend the smallest wavenumber tojk x;y;z e j = 0:0307, a factor of 4 decrease comparing to previous simulation cases that have the smallest wavenumber atjk x;y;z e j = 0:1227. Furthermore, we shift the initial fluctuation modes to 0:2454jk k=? e j 0:4909, and leave a factor of 8 increase in wavelength for inverse cascade to proceed. Additionally, the system used in Run 2 has a much higher resolution, and includes many more wavenumbers that could be resolved on the grid between the range 0:2454 jk k=? e j 0:4909. Therefore, we load the whistler waves on all these wavenumbers. The 139 system then has much more initial fluctuation modes than the 150 modes from the previous forward cascade simulation runs. We assume all modes have the same amplitude, and the condition e = 2:0 is still enforced. Because of the high resolution, the initial spectrum of Run 2 shown in upper panel of Fig. 7.1 shows a clear ring distribution instead of a low resolution square in Fig. 3.1. Note that in the upper right panel of Fig. 7.1, the inverse cascade wavevector space in the center of the spectrum (dark blue area), where there is initially only a low level of thermal fluctuations present, is an elliptical shape, and has two red lines indicating energy extending inwards. This is because that this reduced spectrum is a projection of k x on the k k k y plane, which includes fluctuations from the 45 o , 90 o and 135 o rotation about the k k axis. Using a 3D iso-surface visualization, the inverse cascade wavevector space in the center is confirmed to be a spherical shape. 7.3 Particle-in-cell simulations and discussions Run 2 is the first large scale 3D PIC simulation of decaying whistler turbulence which allows both forward and inverse cascades in a magnetized, homogeneous plasma. The initial ensemble of whistler fluctuations undergoes both forward cascade to shorter wave- lengths and inverse cascade to longer wavelengths. Due to the computational limits, we have run this simulation only toj e jt 80; however, this is long enough to let the inverse cascade reach a quasi-steady state, and to compare the difference between Run 2 and pre- vious forward cascade simulations. Figure 7.1 compares the 2D reduced magnetic fluctuation energy wavevector spectra at t = 0 andj e jt = 78:3 for Run 2. The left column is k x k y spectra representing the plane perpendicular toB o , and the right column isk k k y spectra representing fluc- tuations in a plane containingB o . As in previous forward cascade cases, Run 2 exhibits a gyrotropic spectrum in the perpendicular plane for both forward cascade and inverse 140 Figure 7.1: Reduced magnetic fluctuation energy wavevector spectra, log(jBj 2 =8n e k B T e;t=0 ), from Run 2. The left column is k x k y spectra, the right column isk k k y spectra. Upper panel: att = 0; Lower panel: at atj e jt = 78:3. 141 cascade, as illustrated in the lower left panel. The lower right panel bears the usual quasi- perpendicular wavevector anisotropy for the forward cascade to shorter wavelengths, but is much smoother than the corresponding spectra in the forward cascade such as Fig. 5.8. The much larger number of wavevector modes available in Run 2 better facilitate nonlinear wave-wave interactions, and forms a more realistic smooth spectrum. At smallerjkj in the direction of the inverse cascade, thisk k k y spectra is rather deceptive, showing that the initial empty inverse cascade wavevector space evolves into a quasi-parallel anisotropy in the sensek k > k y (orange color area in the center). This false phenomenon is still due to the projection ofk x atk k k y plane. To get the correct anisotropy for the inverse cascade, we need to turn to more quantitative approaches. Figure 7.2(a-b) show k k k y cuts at k x = 0 from 3D magnetic energy spectrum at the beginning and the end of Run 2, respectively. Note that these are not 2D reduced spectra, therefore have no projection of k x , and much clear inverse cascade wavevector space in the center of the spectra is observable. The initial empty inverse cascade space (dark blue area in the center of panel (a)) is a sphere, rather than an elliptical shape in Fig. 7.1. Through nonlinear wave-wave interactions, the initial fluctuations inverse cascade to larger wavelengths, and the initial inverse cascade space has reduced into an anisotropic and smaller space (dark blue area in the center of panel (b)), in a similar profile as that of the forward cascaded spectrum Figure 7.2(c) shows the anisotropy factor for both the forward and inverse cascade of Run 2, and compare them with Run 1. Run 1 of course only has a forward cascade anisotropy factor. Note that our original definition of fluctuation wavevector anisotropy factor tan 2 B in Eq. 3.4 is inadequate to describe the preferential direction for inverse cascade, because the wavevector origin (k k = 0;k ? = 0) is in the same direction of inverse cascade. So the more initial spectrum inverse cascades, the smaller wavenumberk is used in the summation of Eq. 3.4. In fact, for the inverse cascade, the correct cascaded k =jk o kj, where k o is the initial spectrum wavenumber at t = 0. Alternatively, it 142 (a) B(k x = 0) 2 att = 0 (b) B(k x = 0) 2 atj e jt = 78:3 0 0.2 0.4 0.6 0.8 1 0 50 100 150 200 1 2 3 4 5 6 7 anisotropy factor | e |t Run 1, tan 2 B Run 2, for. cas., tan 2 B Run 2, inv. cas., tan 2 B Run 2, inv. cas., tan 2 inv (c) Anisotropy factor Figure 7.2: (a-b) A k k k y cut at k x = 0 from 3D magnetic energy spectrum, log[jB(k k ;k y ;k x = 0)j 2 =8n e k B T e;t=0 ], at t = 0 andj e jt = 78:3 from Run 2. (c) The anisotropy factor tan 2 B and tan 2 inv of the cascaded energy as a function of time. The green line is tan 2 B from Run 1. The black and blue lines are tan 2 B from the for- ward and inverse cascade of Run 2. The red line is tan 2 inv from Run 2 inverse cascade. Evaluations are taken on spectra over 0:65jk e j 3:0 for Run 1, 0:52jk e j 3:0 for forward cascade of Run 2, andjk e j 0:22 for inverse cascade of Run 2. might be better to calculate the anisotropy factor of the wavevector space in the center of the spectrum (dark blue area in center of panel (b)) which only has a low level of thermal 143 fluctuations. The reciprocal of that factor then indicates the initial fluctuations (red ring in panel (a)) preferential inverse cascade direction. Thus, I propose to define the anisotropy factor tan 2 inv for the inverse cascade as follows, tan 2 inv k k 2 ? jB(k)j 2 k k 2 k jB(k)j 2 (7.1) An isotropic spectrum still corresponds to tan 2 inv = 1:0. A value larger than unity indi- cates that there is less energy ink ? than ink k , suggesting the inverse cascade preferential direction isk k . Panel (c) shows that for the forward cascade, Run 1 and Run 2 have no qualitative difference, i.e., a strong anisotropic wavevector withk ? >k k . Both anisotropies approach to the same level, although Run 2 has a more continuous and much smoother time history of tan 2 B due to its much larger number of modes and broader spectrum band for wave-wave interactions. For the inverse cascade of Run 2, the values of original tan 2 B are roundly around 1.3, indicating there is no preferential wavevector for the cascaded fluctuation. The newly proposed tan 2 inv gives values around 2.5, indicating the inverse cascade pumps slightly more fluctuation energy into smaller k k than smaller k ? . Values of both tan 2 B and tan 2 inv are not large enough to be considered as a strong anisotropy. This indicates that the inverse cascade process does not develop a strong anisotropic wavevector spectrum in Run 2. Figure 7.3 plotsjkj reduced magnetic energy spectrum E(jkj) atj e jt = 67:1 from both Run 1 (red line) and Run 2 (black line). For comparison purpose, the k ? reduced spectrumE(k ? ) of Run 1 is also superimposed on the figure in green line. The definition ofE(jkj) is slightly different toE(k ? ) in Eq. 3.2, such that jBj 2 =8 Z d 3 kjB(k)j 2 =8 = Z djkjjkjE(jkj) (7.2) 144 10 −1 10 0 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 |k| e E(|k|)/K e,t=0 Run 1, E(k ) Run 1, E(|k|) Run 2, E(|k|) |k| 4.6 |k| −5.5 |k| −3.7 Figure 7.3: jkj reduced magnetic fluctuation energy spectraE(jkj) atj e jt = 67:1. The red and black lines are from Run 1 and Run 2, respectively. The green line isk ? reduced spectrumE(k ? ) from Run 1. The grey dashed line represents the initial spectrum of Run 2 att = 0. The blue dashed lines represent power law functions as labeled for comparison against the simulation results. For Run 2, we have to choose thejkj direction instead ofk ? to perform the data reduction, otherwise, the inverse cascade wavevector area would be blocked by the fluctuations from k k ’s non-inverse wavevector space. The figure shows that, for the forward cascade, again, Run 1 and Run 2 have no qualitative difference. They both show a characteristic spectral break atjkj e 1, and have the same power law spectral index. 145 Because of the decaying turbulence nature of our whistler simulation, the inverse cas- cade of Fig. 7.3 shows a positive spectral indexed power law dependance rather than a negative spectral index dependance shown in driven turbulence simulations such as [84]. The inverse cascade spectrum is steeper than the forward cascade spectrum atjkj e 1, but shallower than the spectrum atjkj e 1. Similar spectrum property are shown in 3D EMHD whistler turbulence [102, 24]. Integrating the area below the the cascaded spec- trum curves for both inverse and forward cascade, that is, under the power law dependance /jkj 4:6 for inverse cascade, and under the power law dependance/jkj 3:7 and/jkj 5:5 for forward cascade, we getE frd (jkj)=K e;t=0 = 0:00318 for forward cascade energy, and E inv (jkj)=K e;t=0 = 3:578 10 5 for inverse cascade energy. The ratio between the two is about a factor of 90. Therefore, the total forward cascaded energy is much larger than the total inverse cascaded energy in our whistler turbulence simulations here. Figure 7.4 (a) shows the total magnetic fluctuation energy density as a function of time for Run 1 and Run 2. It demonstrates, once again, that the total fluctuating field energy is dissipated over time, and has no qualitative difference between Run 2 and the corresponding Run 1. Because that the overall damping rate does not increase in the presence of the inverse cascade, it indicates that the dominating process to facilitate the inverse cascade is same as the forward cascade, that is, the nonlinear wave-wave interaction, rather than nonlinear wave-particle interactions, at least for the early time stage where our simulation ends. Figure 7.4(b) illustrate the cascaded magnetic fluctuation energy for both the forward and inverse cascade of Run 2, and the corresponding Run 1, as a function of time in log- arithmic and linear scales, respectively. The wavenumbers satisfying the cascaded range is described in the caption of Fig. 7.2. Again, there has no qualitative difference between the the forward cascade of Run 2 and Run 1: they grow at the same pace and saturate at the same pace. The forward cascaded energy of Run 2 is slightly larger than that of Run 1, which is due to the large number of wave modes in this large scale simulation which serve as the channels for the wave-wave interaction. The difference between them begins 146 0 20 40 60 80 100 0.5 0.6 0.7 0.8 0.9 1 1.1 | B| 2 /8 K e,t=0 (a) magnetic fluctuation energy density 10 20 40 60 80 100 10 −3 10 −2 10 −1 | B| 2 /8 K e,t=0 (b) cascaded magnetic fluctuation energy density in log scale Run 1 Run 2 Run 1, forward cascade Run 2, forward cascade Run 2, inverse cascade 0 20 40 60 80 100 0 0.1 0.2 0.3 d| B| 2 /dt/| B| 2 | e |t (c) cascaded magnetic fluctuation energy density rate of change run 1, forward cascade rate run 2, forward cascade rate run 2, inverse cascade rate Figure 7.4: Magnetic fluctuation energy as functions of time from both Run 1 and Run 2 as labeled. (a) Total fluctuation energy density. (b) Cascaded fluctuation energy density normalized to initial electron kinetic energy density. (c) The rate of cascaded fluctuation energy density. The cascaded range is described in the caption of Fig. 7.2(c). 147 to grow atj e jt 10, reaches maximum atj e jt 60 and continues to decrease after that. This correlates with the same difference shown in Fig. 7.4(a), and implies that the slightly smaller total damping rate for Run 2 in the early time of Fig. 7.4(a) is due to the stronger wave-wave interactions and thus weaker wave-particle damping presented in this large scale simulation. The early-time stage ends at the maximum value of the cascaded energy; at this time the overall damping rate becomes faster than the overall cascade rate so that dissipation dominates the rest of the simulations. Besides, consistent with Fig. 7.3, the forward cascaded energy of Run 2 is much larger than the inverse cascaded energy, such that atj e jt = 67:1 the ratio is roughly 90. Note that although forward cascades saturate atj e jt 60, there is no sign of saturation for inverse cascade. Thus, till the end of our simulation, the overall inverse cascade rate is still faster than the overall damping rate in the inverse cascaded wavevector space. This implies that the wave-wave interaction is still a stronger process than wave-particle damping for inverse cascade. This result is consistent with our linear dispersion prediction in Fig. 5.1: at such large wavelengths, the damping rate of whistler wave is very small for e = 0:1 plasmas. Figure 7.4(c) shows that the forward cascade rate of Run 1 and Run 2 are qualitatively similar. Besides, the inverse cascade rate and the forward cascade rate of Run 2 are also similar. At the end of the early-time stage, forward cascade rates becomes negative, while the inverse cascade rate remains still positive. This is consistent with Fig. 7.4(b). The similar energy transport rate in the early time between the forward and inverse cascade implies that the overall energy difference is likely due to the more wavevector space in shorter wavelengths than that in longer wavelengths which initially only has a low level of thermal fluctuations. If our data sampling rate for Run 2 is as high as Run 1, the forward cascade rate (the black line) will probably exceed the inverse cascade rate (the red line) in the very early time, e.g.,j e jt. 20. Figure 7.5(a) and (b) illustrates the electron kinetic energy density and the kinetic energy anisotropy as a function of time for Run 2 and the corresponding Run 1. The figure 148 1 1.1 1.2 1.3 1.4 1.5 K e /K e,t=0 (a) electron kinetic energy density 0 20 40 60 80 100 1 1.2 1.4 1.6 1.8 (K ||e /K e ) / (K ||e /K e ) t=0 | e |t (b) electron kinetic energy anisotropy Run 1 Run 2 Figure 7.5: (a) Normalized electron kinetic energy density and (b) normalized electron kinetic energy anisotropy as functions of time from both Run 1 and Run 2 as labeled. demonstrate two points. First, consistent with all previous results, there has no qualitative difference between Run 2 and Run 1. Second, the electron heating rate and anisotropy for Run 2 is slightly lower than Run 1 in the early time, consistent with Fig. 7.4(a). This implies again, due to much more wave modes available, Run 2 has stronger wave-wave interaction and thus relatively weaker wave-particle damping than the corresponding Run 1 for the early time stage. When dissipation dominates the late time of the simulations, there is no difference between the two runs. 149 0 0.1 0.2 0.3 0.4 v ! / √ 2v te o f e (v || ) 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 v ! / √ 2v te o 0 0.2 0.4 0.6 0.8 1 v || /v te,t=0 0 0.2 0.4 0.6 0.8 1 v || /v te,t=0 −6 −4 −2 0 2 4 6 0 0.1 0.2 0.3 0.4 v ⊥ / √ 2v te o f e (v ) −6 −4 −2 0 2 4 6 10 −5 10 −4 10 −3 10 −2 10 −1 10 0 v ⊥ / √ 2v te o 0 0 0.2 0.4 0.6 0.8 1 v /v te,t=0 0 0 0.2 0.4 0.6 0.8 1 v /v te,t=0 t=0 Run 1, | e |t=67.1 Run 2, | e |t=67.1 Figure 7.6: Reduced electron velocity distributions as functions ofv k (top panels) and ofv ? (bottom panels) atj e jt = 67:1 from both Run 1 and Run 2 as labeled. The distributions are shown in linear scales in the left panels, and logarithmic scales in the right panels. 150 Figure 7.6 illustrates the late-time reduced electron velocity distributions f e (v k ) and f e (v x ) from Run 2 and the corresponding Run 1. As in Fig. 7.5(b),T k >T ? for both Run 2 and Run 1, implying that the primary energy transfer channel of linear dissipation from whistler fluctuations to the electrons is via Landau damping, with or without the inverse cascade process presented in the simulation. The electron distribution of the two cases are close to each other, demonstrating that there has no difference in terms of electron heating between Run 2 and the corresponding Run 1. 7.4 Conclusions The first fully 3D PIC simulation of decaying whistler turbulence within a simulation box large enough to permit both forward and inverse cascades has been carried out. We impose an initial isotropic spectrum of whistlers upon the system, with wavelengths smaller than the simulation box size, in an initial e = 0:1 and e = 2:0 plasmas. The system exhibits both a forward cascade to shorter wavelengths and inverse cascade to longer wavelengths, and yields transfer of fluctuation energy to a broadband, turbulent spectrum with fluctuation wavevectors preferentiallyk ? >k k in the forward cascade direction, and to a spectrum with fluctuation wavevectorsk k &k ? in the inverse cascade direction. Qualitatively, all properties of the forward cascade in the new run are similar to pre- vious results obtained from a corresponding 3D simulation which only permits forward cascades. These properties include magnetic fluctuation energy wavevector anisotropy, reduced energy spectrum and its electron-scale spectral break and spectral indices, mag- netic energy total dissipation rate and cascaded rate, electron heating rate and kinetic energy anisotropy, and Landau damping as primary linear dissipation channel. These results indicate that the presence of inverse cascade process does not affect qualitative con- clusions previously established from the whistler turbulence forward cascade simulations [19, 35, 20, 21]. 151 Quantitatively, in the new run, turbulence exhibits a smaller energy dissipation rate and electron heating rate, and a larger energy cascaded rate than the corresponding forward- cascade only run at the early-time stage, and such difference keeps reducing at the late- time stage when dissipation dominates. It indicates that relatively stronger wave-wave interactions and weaker wave-particle processes are developed in new run as compared to the forward-cascade only run. However, such difference are not related to the presence of inverse cascade in the new run, but rather due to the its large number of wave modes which facilitate a more rapid transfer of fluctuation energy and stronger wave-wave interactions. Since the overall dissipation rate does not increase in the presence of the inverse cascade, it also indicates that, at least in the early time, the inverse cascade process is primarily driven by the nonlinear wave-wave interaction rather than nonlinear wave-particle interactions, which is same as the forward cascade process. Our results ofjkj reduced magnetic energy spectrum shows that the inverse cascade spectrum is steeper than the forward cascade spectrum atjkj e 1, but shallower than the forward cascade spectrum atjkj e 1. Further computations from the cascaded energy regime suggest that the total energy cascaded forwards is about 90 times as much as the energy cascaded inwards. However, the energy cascades rate of the forward and inverse cascades are similar at early times, indicating that the forward cascade process is not nec- essarily faster than the inverse cascade process in our whistler turbulence evolution simu- lations presented here. 152 CHAPTER 8: CONCLUSIONS AND FUTURE WORK 8.1 Conclusions This thesis presents the first fully 3D PIC simulations of whistler turbulence forward cas- cade in a homogeneous, collisionless plasma with a uniform background magnetic field B o , and the first 3D PIC simulation of whistler turbulence with both forward and inverse cascades. Such computationally demanding research is made possible by the use of mas- sively parallel, high performance EMPIC simulations on state-of-the-art supercomputers. To facilitate this research, our parallel 3D-EMPIC code have been upgraded and extensively optimized. We impose an initial relatively isotropic spectrum of long-wavelength whistler fluctua- tions upon the system, and compute the the temporal evolution of the free decay of these modes to a broadband shorter wavelength regime of turbulence. We have done compar- isons between previous 2D whistler turbulence simulations and our own 3D computations. Effects of initial fluctuation amplitude e and initial electron beta have been studied exten- sively. Relative contributions to energy dissipation and electron heating in whistler turbu- lence from the quasilinear scenario and the intermittency scenario have been estimated. With a simulation box large enough to permit both forward and inverse cascade processes, whistler turbulence evolution direction and wavevector anisotropy have been investigated and compared with corresponding forward-cascade only simulations. This research work has advanced the understanding of whistler turbulence physics and its role in energy cascade and dissipation in the solar wind. It contributes to the ongoing short-wavelength solar wind turbulence and dissipation challenge [67]. Specifically, it provides answers and insights to the following objectives which are also listed in section 153 1.3. What are the 3D characteristics of whistler turbulence evolution and their effects on energy forward cascade, electron velocity distributions and dissipation in solar wind plasmas compared with those in 2D systems? The 3D whistler turbulence exhibits a forward cascade to shorter wavelengths and yields transfer of fluctuation energy to a broadband, anisotropic, turbulent spectrum with wavevectors preferentially quasi-perpendicular to B o . Qualitatively, the forward cascade and wavevector anisotropy properties are similar to those of 2D simulations. However, quantitatively, 3D turbulence develops a stronger anisotropic cascade more rapidly to a larger value, and exhibits much faster rates of energy dissipation than that in 2D simulations. The more modes that are used as initial conditions, the more rapid the cascade and the more strongly anisotropic the turbulence becomes. The biggest difference is that, thek ? reduced magnetic energy spectrum of 3D whistler simulations show a characteristic spectral break at the inverse electron inertial length k ? e ' 1 which is similar to recent observations of short-wavelength solar wind turbulence, whereas the 2D simulation does not show a clear spectral break. We attribute the difference between 3D and 2D runs to the much larger volume of perpendicular wavevector space in the 3D runs, which facilitates the transfer of fluctuation energy toward perpendicular directions. How do 3D whistler turbulence characteristics, such as the energy cascade rates, the electron heating rates, the wavevector anisotropies, the fluctuation spectrum depen- dences, and the different mode couplings, vary under variable fluctuation amplitude e and electron beta e ? At e 0:2 where fluctuation amplitudes are relatively weak, increasing e leads to more anisotropic magnetic wavevector spectra, and stronger anisotropic electron heating in parallel directions. Further increase of e results in a decrease of wavevector and electron 154 temperature anisotropy. This suggests two regimes of distinct dominating dissipation in response to various e . At lower o 0:2, linear Landau damping is the primary dissipation mechanism, whereas at higher o 0:5, nonlinear processes is more likely a dominant dissipation mechanism. Our results show that larger values of e correspond to a faster forward cascade in wavenumber and to a faster rate of electron heating, as well as to a less anisotropic wavevec- tor distribution and to a less anisotropic electron velocity distribution. The overall electron heating yields T k > T ? for all e values, which implies that the primary linear wave- particle interaction is through the Landau resonance, although the cyclotron resonance which should yield aT ? >T k plays a larger role as e increases. However, our results suggest that linear wave-particle interactions play a minor role in shaping the wavevector spectrum, whereas nonlinear wave-wave interactions are overall stronger and faster processes, and ultimately determines the wavevector anisotropy. At early times, the nonlinear wave-wave coupling pumps energy toward shorter wavelengths, anisotropick ? > k k wavevectors, and toward both higher and lower frequencies regime. For the same frequency band, linear wave-particle damping has much stronger damping of the oblique and quasi-perpendicular modes than the parallel propagating modes, especially for lower e cases. The two processes collectively determine the energy spectra (!k) such that the anisotropy decreases with increasing frequency. At late times, wave-wave coupling gradually weakens due to the reduced amplitude of the decaying turbulence, and wave-particle interactions dominate. Linear Landau damping has the strongest damping of whistlers at oblique propagation, which continuously supports the growth of the wavevector anisotropy at late times. All e and e values demonstrate a characteristic spectral break in the k ? reduced energy spectra at k ? e ' 1, with asymptotically constant values of spectral index at nonlinear-process dominated large fluctuation amplitude solar wind plasmas. Because of the inapplicability of Taylor’s hypothesis on whistler waves at e = 1:0 simulations, 155 the spectral indices from our simulated energy wavevector spectrum do not match the frequency spectral indices from observations closely. In contrast, the energy frequency spectra from our simulations exhibit a close similarity to solar wind observations. Electron density fluctuations power spectra also share a close similarity to solar wind observations and MHD predications, both qualitatively and quantitatively. These properties indicate that whistler fluctuations could be the substantial constituent of solar wind turbulence at higher frequencies and short wavelengths, and support the magnetosonic-whistler interpretation of the quasilinear scenario. Is there any effective way to distinguish whistler modes and KAWs in short-wavelength solar wind turbulence and therefore guide future spacecraft measurements and observa- tions? The scaling of the electron-scale spectral break in our simulations suggests that, for whistler turbulence, the spectral break happens atk ? e ' 1:0, instead ofk ? e ' 1:0. If observations could obtain spectra from low e 1:0 time periods, and if kinetic Alfv´ en turbulence indeed has spectral break at k ? e ' 1:0 as argued by [76, 45], we could use this distinct property to distinguish the constituent fluctuations of solar wind turbulence at electron scales. What are the relative contributions to fluctuating energy dissipation from the quasilinear scenario and the intermittency scenario in whistler turbulence? We find that linear damping constitutes a very small part of the total dissipation at e = 0.01, but that increasing e corresponds to an increasing fraction of linear damping. At e = 1:0, linear damping contributes more than 75% of the total dissipation depending on different e . In contrast, increasing e reduces the relative contribution of linear damping to the overall dissipation, and it is more evident at lower e cases. Our results that linear damping represents an intermediate fraction of the total dissipation of whistler turbulence 156 over a wide range of e and e stands in contrast to the results of both [46] and [93], who either found current sheet heating to be several orders of magnitude more efficient, or claimed that the electron Landau resonance is sufficient to account for all electron heating. PDFs and kurtosis of the magnetic field increments in our whistler simulations suggest that intermittency in whistler turbulence generally increases with increasing e and e . Our results indicate that electron scale current sheet generation may be a generic feature of collisionless turbulence. Correlation coefficient calculations imply that the current sheet dissipation should increase with increasing e and e , and suggest that the relative contri- bution between current sheet dissipation and linear damping decrease with increasing e and is not a strong function of e . These properties suggests that the nonlinear dissipation processes involved in our whistler turbulence simulation are primarily associated with dissipation in regions of intermittent current sheet structures. What is the preferential evolution direction, energy cascade rate and wavevector anisotropy of whistler turbulence in a system large enough to permit both forward and inverse cascades? How does the forward cascade process compete against the inverse cascade process, and compare against previous forward cascade process? An even larger scale 3D whistler turbulence simulation exhibits both a forward cascade to shorter wavelengths and an inverse cascade to longer wavelengths, and yields transfer of fluctuation energy to a broadband, turbulent spectrum with fluctuation wavevectors pref- erentially k ? > k k in the forward cascade direction, and to a spectrum with fluctuation wavevectors k k & k ? in the inverse cascade direction. Thejkj reduced magnetic energy spectrum shows that the inverse cascade spectrum is steeper than the forward cascade spec- trum atjkj e 1, but shallower than the forward cascade spectrum atjkj e 1. For our specific computations, the energy cascaded forwards is about 90 times as much as the energy cascaded to longer wavelengths. However, the energy inverse cascade rate is simi- lar to the energy forward cascade rate at early times. The difference of the total cascaded 157 energy is likely due to the more wavevector space in shorter wavelengths than that in longer wavelengths which initially only has a low level of thermal fluctuations. Our results also suggest that, at least in the early time, the inverse cascade process is primarily driven by the nonlinear wave-wave interaction rather than nonlinear wave-particle interactions, which is same as in the forward cascade process. All properties of the forward cascade in this large scale simulation are qualitatively sim- ilar to previous forward cascade results obtained from the corresponding forward-cascade only simulation. Quantitatively, in the new run, turbulence exhibits a smaller energy dissi- pation rate and electron heating rate, and a larger energy cascaded rate than at the early-time stage, indicating a relatively stronger wave-wave interactions and weaker wave-particle processes. Such difference are due to the large number of modes in this large scale simu- lation which facilitate a more rapid transfer of fluctuation energy and stronger wave-wave interactions. Nevertheless, we conclude that the presence of inverse cascade process does not affect qualitative conclusions established from the whistler turbulence forward cascade simulations [19, 35, 20, 21]. 8.2 Future work What is the role of whistler turbulence in short wavelength solar wind observations? Our simulation results support the magnetosonic-whistler interpretation of the quasilinear sce- nario, which hypothesizes that short-wavelength turbulence immediately above the proton- scale spectral break is populated by KAWs, but that whistler fluctuations dominate turbu- lence at still higher frequencies and shorter wavelengths, as illustrated in Fig. 8.1. To directly compare the magnetosonic-whistler and KAWs, it is preferable to expand our current PIC simulations of whistler turbulence to even larger computational grids and thereby encompass the physics of lower frequency, longer wavelength turbulence, including 158 Figure 8.1: A schematic picture of a magnetic energy spectrum of short wavelength solar wind turbulence based on the magnetosonic-whistler interpretation under the framework of quasilinear scenario. both magnetosonic-whistler and KAWs. Such an ensemble simulation of the forward cas- cades of both kinetic Alfven and magnetosonic-whistler turbulence could let magnetosonic- whistler waves and KAWs to directly compete against each other, and conclusions drawn from these should be mostly similar to the actual solar wind condition. We recognize that KAWs propagate at highly oblique angles to the background magnetic field requiring very large simulation box sizes in the direction parallel to B o . So in order to do comparative studies of magnetosonic-whistler waves and KAWs, we must either use large scale sim- ulation box like we did in chapter 7, or compromise by making approximations such as limiting ourselves to 2D simulations. It is also necessary to use artificial proton/electron mass ratio smaller than 1836. Wave-driven turbulence feature should also be implemented besides the current decay- ing turbulence. As has have been done in [45], using an antenna to drive the turbulence cascade should be utilized to create a more quasi-steady state for data analyses than our 159 current two phases in turbulence evolution. This technique will in general lead to simula- tion results which more closely represent the actual solar wind turbulence. 160 APPENDIX A: APPENDIX A.1 Background in plasma waves and instabilities Electromagnetic waves propagate at speed of light in vacuum. The frequency! as a func- tion of wavevector k, which is usually called a “dispersion relation”, is simple! 2 =k 2 c 2 . This indicates electromagnetic waves of different frequency travel at same phase velocity v =!=k, a property called “dispersionless”. In contrast, when electromagnetic waves propagate in a plasma, there can be strong interactions between the plasma and the fluctuating fields, creating dispersion. Dispersion relations can been categorized into many types and become rather complex, depending on the plasma properties such as hot or cold, unmagnetized or magnetized, and on various assumptions made on the system. In general, each type of dispersion relation represents a certain type of “plasma wave”, and the phase velocity can depend on wave frequency, that is “dispersive”. The wave amplitude now could be attenuated or excited depending on the range of!. As an simple example, consider an electromagnetic wave in an unmagnetized plasma propagating in the x direction. The dispersion relation is now ! 2 = k 2 c 2 +! 2 e . The phase velocity is larger than the speed of light in vacuum. The wave is generally dispersive, and must have frequency larger than! e to propagate in a plasma with field pro- portional to exp(ikx), or it will exhibits a “cut-off” at! e with field attenuated proportional to exp(kx). Damping in plasmas refers to transferring energy from waves to plasmas, so the fluctua- tions amplitudes are attenuated and the plasma particles are heated. In contrast, microinsta- bilities refer to transferring free kinetic energy from the plasma to the fluctuations, and the term “instability” is usually applied only to the case that a rapid growth of wave amplitude and transfer of energy and momentum from sources to waves are present. 161 Solution of the linear Vlasov equation in a homogeneous plasma usually involves a Laplace transform in time, which yields a complex frequency! expressed as! =! r +i , where ! r and are real and imaginary part of the frequency. Assuming the fluctuation amplitude is sufficiently small, the magnetic field could be expressed as single frequency small amplitude waves, such asB =B 0 +Be i(kr!t) , and the amplitude of small fluctu- ation term isBe t . If the wave is either to be damped or to grow, it must satisfy certain resonance condition that the wave and particles are able to interact. If is negative, the wave is damped; and if is positive, the fluctuation will grow exponentially with time. is then called the damping rate or the growth rate depending on whether it is a damp- ing or an instability. We describe several important wave-particle damping mechanisms that are related to short-wavelength solar wind turbulence in section 1.1.2, and here only discuss plasma instability. The above linear plasma wave approach is valid only at small amplitudes, the exponentially-growing wave due to instabilities could only satisfies this approximation for a short time compared with 1= ; at longer time, the wave amplitude becomes too large that makes the linear theory above no longer valid [17]. When the wave amplitude is large enough, nonlinear effects are introduced to the sys- tem, and waves with different frequencies and wavenumbers start to interact. That leads to excitation of waves with frequencies and wavelengths in wide intervals, not necessarily in resonance with the external source. This is the onset of plasma turbulence. As turbulence is always accompanied by particle dissipation, the energy cascade process will eventually limit the continuous growth of waves. Fluctuations are damped, the medium is heated, and the system is brought back to steady state. Therefore, linking the plasma instability and turbulence evolution may gain insight on energy cascade process in the solar wind. 162 A.2 3D-EMPIC code sequential algorithm and normaliza- tion The basic sequential algorithm for our EMPIC code is as follows 1. Set the initial conditions of the particles and fields. The initial conditions must satisfy the two divergence Maxwell Eqs. 2.4 and 2.5. 2. Particle move: a) Interpolate the electromagnetic field on the particle position to obtain the force on each particle (gather); b) Update the particle velocity and position from Eqs. 2.11 and 2.12. 3. Current deposit (scatter): Calculate the charge carried by particles across cell surfaces within the time step to obtain the current fluxJdt through each cell surface. Steps (2) and (3) together form the particle push stage of the code. 4. Field update: Solve the two curl Maxwell equations by finite-difference leapfrogging Eqs. 2.9 and 2.10 to update the electromagnetic field Finally, the EMPIC model needs to be normalized for minimizing the round-off error. The inverse electron plasma frequency ! e is the unit of time, and the unit of length usually equals to Debye length De . The full normalization units are x = ^ x t = ^ t !e v = ! e ^ v De = ^ v te q =e^ q m =m e ^ m n =n o ^ n =en o ^ E = me! 2 e e ^ E B = me! 2 e e ^ B J =en o ! e ^ J e =! e ^ e ^ B = ^ c ^ e where ^ indicates the normalized value. 163 A.3 Some definitions used in turbulence diagnostics The fluctuating magnetic field increment as a function of spatial separation length r is denoted as r j B i =B i (j +r)B i (j) (A.1) wherei;j could be eitherx;y, orz direction. Kurtosis of magnetic field increments is then defined as<= 4 ij >: <= 4 ij >= < r j B 4 i > 4 (A.2) where =<j r j B i j 2 > 1=2 is the non-centralized standard deviation and< > denotes an ensemble average over three dimensional space. A Gaussian distribution yields a kur- tosis of 3; departures from a Gaussian distribution correspond to increasing values of the kurtosis and are a measure of what has been called turbulent intermittency. The probability density function (PDF) of a continuous distribution is a function that describes the relative likelihood for this random variable to take on a given value, and is defined as PDF[aXb] = Z b a f X (x)dx (A.3) wheref X is density function and is a non-negative Lebesgue integrable function. The correlation coefficient of two quantitiesX andY is denoted as (X;Y ) = cov(X;Y ) X Y (A.4) wherecov denotes covariance, and is the standard deviation. The correlation coefficient is a normalized value between -1 and 1. Value close to 1 indicates that there is a positive linear relationship between the data, while values close to -1 indicate a negative linear rela- tionship. Value close to or equal to 0 suggests there is no linear relationship between the 164 data. Note that the correlation reflects the non-linearity and direction of a linear relation- ship, but not the slope of that relationship, nor many aspects of nonlinear relationships. The electron temperature is calculated as T je =m e v 2 te j =m e ( 1 N N X i=1 v 2 i j v j ) (A.5) where N is number of electron superparticles, j could be either x, y or z direction, and v j = 1 N P N i=1 v i j . The common method of applying spectral techniques to equations that depend on space or time is to work in Fourier domain as a function of wavenumber or frequency, using Fourier transforms on the space and time variable. This is a well known technique which provides for a reliable solution to problems with periodic functions. Any physical process can be described either in the physical space domain, by the values of some quantityu as a function of physical space x, e.g., u(x), or else in Fourier domain, where the process is specified by giving its amplitudeU as a function of wavenumberk, that isU(k), with 1 < k <1. To go back and forth between two domains, the 1D Fourier and inverse Fourier transform equations are introduced [73]: U(k) = Z 1 1 u(x)e ikx dx (A.6) u(x) = 1 2 Z 1 1 U(k)e ikx dk (A.7) These equations are continuous and in computational applications, we have to use the dis- cretized version of the Fourier transform. A common situation is to have a function sampled 165 at evenly spaced discrete intervals, x. For any sampling interval x, the Nyquist critical wavenumber is given by: k c = 1 2x (A.8) This lead to an important theorem for discrete sampling known as the sampling theorem: If a continuous functionu(x), sampled at an interval x, is bandwidth limited to wavenum- bers smaller in magnitude thank c , then the functionu(x) is completely determined by its samplesu n . The sampling theorem is significant because it allows for the entire informa- tion content of a signal to be reproduced provided it is sampled at sufficient rate. Suppose we have an even function ofN + 1 consecutive samples: u j u(x j );x j =jx;j = 0; 1; 2;:::;N (A.9) such that the sampling interval is x, the set of pointsx j are referred to as nodes or grid points and the period of the quantityu j is just the domain lengthNx. The wavenumber is defined as: k n = 2n Nx ; where mode numbern = N 2 ;:::; 0;:::; N 2 (A.10) The discrete Fourier transform (DFT) is of theN + 1 pointsu j , denoted byU n U n = N X j=0 u j e iknx j = N X j=0 u j e i 2n N j (A.11) Likewise, the formula for the inverse discrete Fourier transform can be derived as: u j = 1 N N=2 X n=N=2 U n e iknx j = 1 N N=2 X n=N=2 U n e i 2n N j (A.12) 166 Similar to the one-dimensional case, we define the two-dimensional DFT and inverse DFT as following [16]: U nx;ny = Nx X jx=0 Ny X jy =0 u jx;jy e iknx x jx e ikny y jy (A.13) u jx;jy = 1 N x 1 N y Nx=2 X nx=Nx=2 Ny=2 X ny =Ny=2 U nx;ny e iknx x jx e ikny y jy (A.14) where k nx = 2nx Nx ; k ny = 2ny Ny are corresponding wavenumbers in ^ k x and ^ k y direc- tion, respectively, and x jx = j x x, y jy = j y y are discrete grid points. Here n x = Nx 2 ;:::; Nx 2 ;n y = Ny 2 ;:::; Ny 2 ;j x = 0; 1;:::;N x ;j y = 0; 1;:::;N y For 1D case, the DFT maps N + 1 complex number U n into N + 1 complex num- bers u j . As the number of input data points increases, the output frequency spectrum becomes closer to continuous spectrum. The amount of computation involved in comput- ing an N-point DFT is O(N 2 ) process. With N getting bigger, the computational time will be dramatically increased. In modern scientific computing (also in this thesis), a far more efficient algorithm known as fast Fourier transform (FFT) which is derived from DFT is usually applied instead of DFT. The basis idea of FFT is that a DFT of length N can be expressed as the sum of two DFTs, each of lengthN=2, cutting the processing time in half. This process can be repeated recursively on the DFT of the length N=2 and so on until only two components left. The computational time for an N-point FFT isO(log N 2 ), which is much more efficient comparing with DFT when dealing with a largerN. How- ever, it is usually necessary thatN be an integer power of 2 for a FFT in order to do binary manipulation. The detailed derivation of the FFT algorithm can be found in [73]. 167 Reference List [1] National Aeronautics and Space Administration. 2010 science plan. 40, 2010. [2] O. Alexandrova, V . Carbone, P. Veltri, and L. Sorriso-Valvo. Small-scale energy cascade of the solar wind turbulence. Astrophys. J., 674(1153), 2008. [3] O. Alexandrova, C. Lacombe, A. Mangeney, R. Grappin, and M. Maksimovic. Solar wind turbulent spectrum at plasma kinetic scales. Astrophys. J., 760(121), 2012. [4] O. Alexandrova, J. Saur, C. Lacombe, A. Mangeney, J. Mitchell, S. J. Schwartz, and P. Robert. Universality of solar wind turbulent spectrum from mhd to electron scales. Phys. Rev. Lett., 103(165003), 2009. [5] S. D. Bale, Kellog P. J., F.S. Mozer, T. S. Horbury, and H. Reme. Magnetic fluc- tuation power near proton temperature anisotropy instability thresholds in the solar wind. Phys. Rev. Lett., 103(211101), 2009. [6] H. J. Beinroth and F. M. Neubauer. Properties of whistler mode waves between 0.3 and 1.0 au from helios observations. J. Geophys. Res., 86(7755), 1981. [7] C. K. Birdsall and A. B. Langdon. Plasma physics via computer simulation. McGraw-Hill, 1991. [8] J. Birn and E. R. Priest. Reconnection of Magnetic Fields: Magnetohydrodynamics and Collisionless Theory and Observations. Cambridge University Press, 2006. [9] D. Biskamp and W.-C. M¨ uller. Scaling properties of three-dimensional isotropic magnetohydrodynamic turbulence. Phys. Plasmas, 7(4889), 2000. [10] D. Biskamp, E. Schwarz, A. Zeiler, A. Celani, and J. F. Drake. Electron magnetohy- drodynamic turbulence. Phys. Plasmas, 6(751), 1999. [11] J. E. Borovsky and M. H. Denton. No evidence for heating of the solar wind at strong current sheets. Astrophys. J., 739(L61), 2011. [12] S. Bourouaine, O. Alexandrova, E. Marsch, and M. Maksimovic. On spectral breaks in the power spectra of magnetic fluctuations in fast solar wind between 0.3 and 0.9 au. Astrophys. J., 749(102), 2012. 168 [13] K. J. Bowers, B. J. Albright, L. Yin, B. Bergen, and T. J. T. Kwan. Ultrahigh per- formance three-dimensional electromagnetic relativistic kinetic plasma simulation. Phys. Plasmas, 15(055703), 2008. [14] O. Buneman. TRISTAN: The 3-D Electromagnetic Particle Code. Computer Space Plasma Physics: Simulation Techniques and Software, Edited by H. Matsumoto and Y. Omura, pages 67–84. Terra Scientific Publishing Company, 1993. [15] E. Camporeale and D. Burgess. The dissipation of solar wind turbulent fluctuations at electron scales. Astrophys. J., 730(114), 2011. [16] G. Chae. Numerical Simulation of Ion Waves in Dusty Plasmas. PhD thesis, Virginia Polytechnic Institute and State University, 2000. [17] O. Chang. Numerical simulation of ion-cyclotron turbulence generated by artifi- cial plasma cloud release. Master’s thesis, Virginia Polytechnic Institute and State University, 2009. [18] O. Chang. CSCI 590 directed research report: Streaming algorithm of three dimen- sional electromagnetic particle-in-cell method on graphical processing units. Tech- nical report, University of Southern California, 2013. [19] O. Chang, S. P. Gary, and J. Wang. Whistler turbulence forward cascade: Three- dimensional particle-in-cell simulations. Geophys. Res. Lett., 38(L22102), 2011. [20] O. Chang, S. P. Gary, and J. Wang. Whistler turbulence at variable electron beta: Three-dimensional particle-in-cell simulations. J. Geophys. Res., 118(1), 2013. [21] O. Chang, S. P. Gary, and J. Wang. Whistler turbulence dissipation: Three- dimensional particle-in-cell simulations. Geophys. Res. Lett., 2013, in preparation. [22] C. H. K. Chen, T. S. Horbury, A. A. Schekochihin, R. T. Wicks, O. Alexandrova, and J. Mitchell. Anisotropy of solar wind turbulence in the dissipation range. Phys. Rev. Lett., 104(255003), 2010. [23] C. H. K. Chen, C. S. Salem, J. W. Bonnell, F. S. Mozer, and S. D. Bale. Density fluctuation spectrum of solarwind turbulence between ion and electron scales. Phys. Rev. Lett., 109(035001), 2012. [24] J Cho. Magnetic helicity conservation and inverse energy cascade in electron mag- netohydrodynamic wave packets. Phys. Rev. Lett., 106(191104), 2011. [25] J. Cho and A. Lazarian. The anisotropy of electron magnetohydrodynamic turbu- lence. Astrophys. J., 615(L41), 2004. [26] J. Cho and A. Lazarian. Simulations of electron magnetohydrodynamic turbulence. Astrophys. J., 701(236), 2009. 169 [27] C. Crabtree, L. Rudakov, G. Ganguli, M. Mithaiwala, V . Galinsky, and V . Shevchenko. Weak turbulence in the magnetosphere: Formation of whistler wave cavity by nonlinear scattering. Phys. Plasmas, 19(032903), 2012. [28] B. Eliasson and K. Papadopoulos. Numerical study of mode conversion between lower hybrid and whistler waves on short-scale density striations. J. Geophys. Res., 113(A09315), 2008. [29] S. Galtier. Wave turbulence in incompressible hall magnetohydrodynamics. J. Plasma Phys., 72(721), 2006. [30] S. Galtier and A. Bhattacharjee. Anisotropic weak whistler wave turbulence in elec- tron magnetohydrodynamics. Phys. Plasmas, 10(3065), 2003. [31] S. Galtier and E. Buchlin. Multiscale hall-magnetohydrodynamic turbulence in the solar wind. Astrophys. J., 656(560), 2007. [32] G. Ganguli, L. Rudakov, W. Scales, J. Wang, and M. Mithaiwala. Three dimensional character of whistler turbulence. Phys. Plasmas, 17(52310), 2010. [33] S. P. Gary. Theory of space plasma microinstabilities. Cambridge University Press, 1993. [34] S. P. Gary. Test for wavevector anisotropies in plasma turbulence cascades. Astro- phys. J., 769(36), 2013. [35] S. P. Gary, O. Chang, and J. Wang. Forward cascade of whistler turbulence: Three- dimensional particle-in-cell simulations. Astrophys. J., 755(142), 2012. [36] S. P. Gary, K. Liu, and D. Winske. Whistler anisotropy instability at low electron: Particle-in-cell simulations. Phys. Plasmas, 18(082902), 2011. [37] S. P. Gary, S. Saito, and H. Li. Cascade of whistler turbulence: particle-in-cell simulations. Geophys. Res. Lett., 35(L02104), 2008. [38] S. P. Gary and C. W. Smith. Short-wavelength turbulence in the solar wind: Linear theory of whistler and kinetic alfv´ en fluctuations. J. Geophys. Res., 114(A12105), 2009. [39] S. P. Gary and J. Wang. Whistler instability: electron anisotropy upper bound. J. Geophys. Res., 101(10749), 1996. [40] S. Ghosh, E. Siregar, D. A. Roberts, and M. L. Goldstein. Simulation of high- frequency solar wind power spectra using hall magnetohydrodynamics. J. Geophys. Res., 101(2493), 1996. 170 [41] J. He, C. Tu, E. Marsch, and S. Yao. Do oblique alfv´ en/ion-cyclotron or fast- mode/whistler waves dominate the dissipation of solar wind turbulence near the proton inertial length? Astrophys. J., 745(L8), 2012. [42] T. S. Horbury, M. A. Forman, and S. Oughton. Spacecraft observations of solar wind turbulence: an overview. Phys. Plasmas, 47(B703), 2005. [43] G. G. Howes, S. C. Cowley, W. Dorland, G. W. Hammett, E. Quataert, and A. A. Schekochihin. A model of turbulence in magnetized plasmas: Implications for the dissipation range in the solar wind. J. Geophys. Res., 113(A05103), 2008. [44] G. G. Howes, W. Dorland, S. C. Cowley, G. W. Hammett, E. Quataert, A. A. Schekochihin, and T. Tatsuno. Kinetic simulations of magnetized turbulence in astrophysical plasmas. Phys. Rev. Lett., 100(065004), 2008. [45] G. G. Howes, J. M. TenBarge, W. Dorland, E. Quataert, A. A. Schekochihin, R. Numata, and T. Tatsuno. Gyrokinetic simulations of solar wind turbulence from ion to electron scales. Phys. Rev. Lett., 107(035004), 2011. [46] H. Karimabadi, V . Roytershteyn, M. Wan, W. H. Matthaeus, W. Daughton, P. Wu, M. Shay, B. Loring, J. Borovsky, E. Leonardis, S. C. Chapman, and T. K. M. Nakamura. Coherent structures, intermittent turbulence, and dissipation in high- temperature plasmas. Phys. Plasmas, 20(012303), 2013. [47] K. H. Kiyani, S. C. Chapman, F. Sahraoui, B. Hnat, O. Fauvarque, and Yu. V . Khotyaintsev. Enhanced magnetic compressibility and isotropic scale invariance at sub-ion larmor scales in solar wind turbulence. Astrophys. J., 763(10), 2013. [48] K. H. Kiyani, Yu. V . Khotyaintsev S. C. Chapman, M. W. Dunlop, and F. Sahraoui. Global scale-invariant dissipation in collisionless plasma turbulence. Phys. Rev. Lett., 103(075006), 2009. [49] V . Krishan and S. M. Mahajan. Magnetic fluctuations and hall magnetohydrody- namic turbulence in the solar wind. J. Geophys. Res., 109(A11105), 2004. [50] R. J. Leamon, C. W. Smith, N. F. Ness, W. H. Matthaeus, and H. K. Wong. Obser- vational constraints on the dynamics of the interplanetary magnetic field dissipation range,. J. Geophys. Res., 103(4775), 1998. [51] D. Lengyel-Frey, R. A. Hess, R. J. MacDowall, R. G. Stone, N. Lin, A. Balogh, and R. Forsyth. Ulysses observations of whistler waves at interplanetary shocks and in the solar wind. J. Geophys. Res., 101(27555), 1996. [52] A. T. Y . Lui, P. H. Yoon, Chinook Mok, and Chang-Mo Ryu. Inverse cascade feature in current disruption. J. Geophys. Res., 113(A00C06), 2008. 171 [53] B. Marder. A method for incorporating Gauss’ law into electromagnetic PIC codes. J. Comput. Phys., 68(48), 1987. [54] W. H. Matthaeus, M. L. Goldstein, and D. A. Roberts. Evidence for the presence of quasi-two-dimensional nearly incompressible fluctuations in the solar wind. Phys. Rev. Lett., 95(20673), 1990. [55] W. H. Matthaeus, S. Oughton, S. Ghosh, and M. Hossain. Scaling of anisotropy in hydromagnetic turbulence. Phys. Rev. Lett., 81(2056), 1998. [56] W. H. Matthaeus and M. Velli. Who needs turbulence? a review of turbulence effects in the heliosphere and on the fundamental process of reconnection. Space Sci. Rev., 160(145), 2011. [57] R. Meyrand and S. Galtier. A universal law for solar wind turbulence at electron scales. Astrophys. J., 721(1421), 2010. [58] M. Mithaiwala, L. Rudakov, C. Crabtree, and G. Ganguli. Co-existence of whistler waves with kinetic alfven wave turbulence for the high-beta solar wind plasma. Phys. Plasmas, 19(102902), 2012. [59] M. Mithaiwala, L. Rudakov, G. Ganguli, and C. Crabtree. Weak turbulence theory of the nonlinear evolution of the ion ring distribution. Phys. Plasmas, 18(055710), 2011. [60] B. Moon, H.V . Jagadish, C. Faloutsos, and J. H. Saltz. Analysis of the clustering properties of the hilbert space-filling curve. IEEE Trans. Knowledge and Data Eng., 13(124), 2001. [61] A. Nakano, R. K. Kalia, K. Nomura, A. Sharma, P. Vashishta, F. Shimojo, A. C. T. van Duin, W. A. Goddard, R. Biswas, and D. Srivastava. A divide-and- conquer/cellular-decomposition framework for million-to-billion atom simulations of chemical reactions. Comp. Mat. Sci., 38(642), 2007. [62] Y . Narita, S. P. Gary, S. Saito, K. H. Glassmeier, and U. Motschmann. Dispersion relation analysis of solar wind turbulence. Geophys. Res. Lett., 38(L05101), 2011. [63] D. R. Nicholson. Introduction to plasma theory. McGraw-Hill, 1983. [64] K. T. Osman, W. H. Matthaeus, A. Greco, and S. Servidio. Evidence for inhomoge- neous heating in the solar wind. Astrophys. J., 727(L11), 2011. [65] K. T. Osman, W. H. Matthaeus, B. Hnat, and S. C. Chapman. Kinetic signatures and intermittent turbulence in the solarwind plasma. Phys. Rev. Lett., 108(261103), 2012. 172 [66] K. T. Osman, W. H. Matthaeus, M. Wan, and A. F. Rappazzo. Intermittency and local heating in the solar wind. Phys. Rev. Lett., 108(261102), 2012. [67] T. N Parashar and C. Salem. Turbulent dissipation challenge: A community driven effort. Phys. Plasmas, 2013, in press. [68] S. Perri, M. L. Goldstein, J. C. Dorelli, and F. Sahraoui. Detection of small- scale structures in the dissipation regime of solar-wind turbulence. Phys. Rev. Lett., 109(191101), 2012. [69] J. J. Podesta. Solar wind turbulence: Advances in observations and theory. Advances in Plasma Astrophysics, 274(295), 2010. [70] J. J. Podesta. The need to consider ion bernstein waves as a dissipation channel of solar wind turbulence. J. Geophys. Res., 117(A07101), 2012. [71] J. J. Podesta, J. E. Borovsky, and S. P. Gary. A kinetic alfv´ en wave cascade subject to collisionless damping cannot reach electron scales in the the solar wind at 1 au. Astrophys. J., 712(685), 2010. [72] J. J. Podesta and J. M. TenBarge. Scale dependence of the variance anisotropy near the proton gyroradius scale: Additional evidence for kinetic alfv´ en waves in the solar wind at 1 au. J. Geophys. Res., 117(A10106), 2012. [73] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical recipes in FORTRAN: the art of scientific computing. Cambridge University Press, 2nd edition, 1992. [74] V . Przebinda. Vertical optimization of particle in cell simulation. Master’s thesis, University of Colorado, 2003. [75] L. Rudakov, M. Mithaiwala, G. Ganguli, and C. E. Crabtree. Linear and non-linear landau resonance of kinetic alfv´ en waves: Consequences for electron distribution and wave spectrum in the solar wind. Phys. Plasmas, 18(012307), 2011. [76] F. Sahraoui, M. L. Goldstein, G. Belmont, P. Canu, and L. Rezeau. Three dimen- sional anisotropick spectra of turbulence at subproton scales in the solar wind. Phys. Rev. Lett., 105(131101), 2010. [77] F. Sahraoui, M. L. Goldstein, P. Robert, and Yu. V . Khotyaintsev. Evidence of a cascade and dissipation of solar-wind turbulence at the electron gyroscale. Phys. Rev. Lett., 102(231102), 2009. [78] S. Saito and S. P. Gary. Beta dependence of electron heating in decaying whistler turbulence: Particle-in-cell simulations. Phys. Plasmas, 19(012312), 2012. 173 [79] S. Saito, S. P. Gary, H. Li, and Y . Narita. Whistler turbulence: particle-in-cell simu- lations. Phys. Plasmas, 15(102305), 2008. [80] S. Saito, S. P. Gary, and Y . Narita. Wavenumber spectrum of whistler turbulence: particle-in-cell simulation. Phys. Plasmas, 17(122316), 2010. [81] C. S. Salem, G. G. Howes, D. Sundkvist, S. D. Bale, C. C. Chaston, C. H. K. Chen, and F. S. Mozer. Identification of kinetic alfv´ en wave turbulence in the solar wind. Astrophys. J., 745(L9), 2012. [82] D. Schriver, M. Ashour-Abdalla, F. V . Coroniti, J. N. LeBoeuf, V . Decyk, P. Travnicek, O. Santol´ ık, and D. Winningham. Generation of whistler mode emis- sions in the inner magnetosphere: An event study. J. Geophys. Res., 115(A00F17), 2010. [83] D. Shaikh. Whistler wave cascades in solar wind plasma. Mon. Not. R. Astron. Soc., 395(2292), 2009. [84] D. Shaikh and G. P. Zank. Driven dissipative whistler wave turbulence. Phys. Plas- mas, 12(122310), 2005. [85] D. Shaikh and G. P. Zank. Spectral features of solar wind turbulent plasma. Mon. Not. R. Astron. Soc., 400(1881), 2009. [86] J. V . Shebalin, W. H. Matthaeus, and Montgomery D. Anisotropy in mhd turbulence due to a mean magnetic field. Phys. Plasmas, 29(515), 1983. [87] C. W. Smith, K. Hamilton, B. J. Vasquez, and R. J. Leamon. Dependence of the dissipation range spectrum of interplanetary magnetic fluctuations on the rate of energy cascade. Astrophys. J., 645(L85), 2006. [88] C. W. Smith, B. J. Vasquez, and J. V . Hollweg. Observational constraints on the role of cyclotron damping and kinetic alfv´ en waves in the solar wind. Astrophys. J., 745(8), 2012. [89] L. Sorriso-Valvo. Understanding kappa distributions in space science. In Presenta- tions of 2013 American Geophysical Union Fall Meeting, 2013. [90] O. Stawicki, S. P. Gary, and H. Li. Solar wind magnetic fluctuation spectra: Disper- sion versus damping. J. Geophys. Res., 106(A5), 2001. [91] V . A. Svidzinski, H. Li, H. A. Rose, B. J. Albright, and K. J. Bowers. Particle in cell simulations of fast magnetosonic wave turbulence in the ion cyclotron frequency range. Phys. Plasmas, 16(122310), 2009. [92] G. I. Taylor. The spectrum of turbulence. Proc. R. Soc. A, 164(476), 1938. 174 [93] J. M. TenBarge and G. G. Howes. Current sheets and collisionless damping in kinetic plasma turbulence. Astrophys. J., 771(L27), 2013. [94] J. M. TenBarge, J. J. Podesta, K. G. Klein, and G. G. Howes. Interpreting magnetic variance anisotropy measurements in the solar wind. Astrophys. J., 753(107), 2012. [95] O. P. Verkhoglyadova, B. T. Tsurutani, and G. S. Lakhina. Properties of obliquely propagating chorus. J. Geophys. Res., 115(A00F19), 2010. [96] J. Villasenor and O. Buneman. Rigorous charge conservation for local electromag- netic field solvers. Comput. Phys. Commun., 69(306), 1992. [97] M. Wan, W. H. Matthaeus, H. Karimabadi, V . Roytershteyn, M. Shay, P. Wu, W. Daughton, B. Loring, and S. C. Chapman. Intermittent dissipation at kinetic scales in collisionless plasma turbulence. Phys. Rev. Lett., 109(195001), 2012. [98] J. Wang, Y . Cao, R. Kafafy, J. Pierru, and V . K. Decyk. Simulations of ion thruster plume-spacecraft interactions on parallel supercomputer. IEEE Trans. Plasma Sci., 34(2148), 2006. [99] J. Wang, D. Kondrashov, P. Liewer, and S. R. Karmesin. Three-dimensional deformable-grid electromagnetic particle-in-cell for parallel computers. J. Plasma Phys., 61(367), 1999. [100] J. Wang, P. Liewer, and V . Decyk. 3d electromagnetic plasma particle simulations on a mimd parallel computer. Comput. Phys. Commun., 87(35), 1995. [101] J. Wang, P. Liewer, and E. Huang. Three-dimensional electromagnetic particle-in- cell with monte carlo collision simulations on three mimd parallel computers. J. Supercomputing, 10(331), 1997. [102] C. J. Wareing and R. Hollerbach. Cascades in decaying three-dimensional electron magnetohydrodynamic turbulence. J. Plasma Phys., 76(117), 2010. [103] J. Weiland and H. Wilhelmsson. Coherent Nonlinear Interaction of Waves in Plas- mas. Pergamon Press, 1977. [104] www2.cisl.ucar.edu/resources/yellowstone, 2013. [105] www.mathworks.com/help/wavelet/ref/scal2frq.html, 2013. [106] www.top500.org. Top500 list - June 2013. [107] J. S. Yee. Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media. IEEE Trans. Antennas Propagat., 14(302), 1966. 175

Abstract (if available)

Abstract The objective of this dissertation is to study the physics of whistler turbulence evolution and its role in energy transport and dissipation in the solar wind plasmas through computational and theoretical investigations. This dissertation presents the first fully three-dimensional (3D) particle-in-cell (PIC) simulations of whistler turbulence forward cascade in a homogeneous, collisionless plasma with a uniform background magnetic field Βₒ, and the first 3D PIC simulation of whistler turbulence with both forward and inverse cascades. Such computationally demanding research is made possible through the use of massively parallel, high performance electromagnetic PIC simulations on state-of-the-art supercomputers. ❧ Simulations are carried out to study characteristic properties of whistler turbulence under variable solar wind fluctuation amplitude (εₑ) and electron beta (βₑ), relative contributions to energy dissipation and electron heating in whistler turbulence from the quasilinear scenario and the intermittency scenario, and whistler turbulence preferential cascading direction and wavevector anisotropy. ❧ The 3D simulations of whistler turbulence exhibit a forward cascade of fluctuations into broadband, anisotropic, turbulent spectrum at shorter wavelengths with wavevectors preferentially quasi-perpendicular to Βₒ. The overall electron heating yields T∥ > Τ⊥ for all εₑ and βₑ values, indicating the primary linear wave-particle interaction is Landau damping. But linear wave-particle interactions play a minor role in shaping the wavevector spectrum, whereas nonlinear wave-wave interactions are overall stronger and faster processes, and ultimately determine the wavevector anisotropy. ❧ Simulated magnetic energy spectra as function of wavenumber show a spectral break to steeper slopes, which scales as κ⊥λₑ ≃ 1 independent of βₑ values, where λₑ is electron inertial length, qualitatively similar to solar wind observations. Specific spectral indices from simulated wavevector energy spectra do not match the frequency spectral indices from observations due to the inapplicability of Taylor's hypothesis. In contrast, the direct comparison of simulated frequency energy spectra and solar wind observations shows a closer similarity. Electron density fluctuations power spectra also exhibit a close similarity to solar wind observations and MHD predications, both qualitatively and quantitatively. ❧ Linear damping represents an intermediate fraction of the total dissipation of whistler turbulence over a wide range of βₑ and εₑ. The relative importance of linear damping by comparison to nonlinear dissipation increases with increasing βₑ but decreases with increasing εₑ. Correlation coefficient calculations imply that the nonlinear dissipation processes in our simulation are primarily associated with dissipation in regions of intermittent current sheet structures. The simulation results suggest that whistler fluctuations could be the substantial constituent of solar wind turbulence at higher frequencies and short wavelengths, and support the magnetosonic-whistler interpretation of the quasilinear scenario. ❧ An even larger scale 3D whistler turbulence simulation exhibits both a forward cascade to shorter wavelengths with wavevectors preferentially κ⊥ > κ∥, and an inverse cascade to longer wavelengths with wavevectors κ∥ ≳ κ⊥. The inverse cascade process is primarily driven by the nonlinear wave-wave interaction. It is shown that the energy inverse cascade rate is similar to the energy forward cascade rate at early times although the overall energy in the two cascades is very different. The presence of inverse cascade process does not affect qualitative conclusions established from the whistler turbulence forward cascade simulations.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

Grid-based Vlasov method for kinetic plasma simulations

PDF

Numerical and experimental investigations of dust-plasma-asteroid interactions

PDF

Large eddy simulations of turbulent flows without use of the eddy viscosity concept

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

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

PDF

Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction

PDF

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

PDF

A subgrid-scale model for large-eddy simulation based on the physics of interscale energy transfer in turbulence

PDF

Modeling and analysis of parallel and spatially-evolving wall-bounded shear flows

Asset Metadata

Creator Chang, Ouliang (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Astronautical Engineering

Publication Date 08/31/2013

Defense Date 08/02/2013

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

Tag kinetic simulation,OAI-PMH Harvest,particle-in-cell,plasmas,solar wind,whistler turbulence

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Wang, Joseph (committee chair), Gary, S. Peter (committee member), Gruntman, Michael (committee member), Nakano, Aiichiro (committee member)

Creator Email ouliang@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-327291

Unique identifier UC11293925

Identifier etd-ChangOulia-2028.pdf (filename),usctheses-c3-327291 (legacy record id)

Legacy Identifier etd-ChangOulia-2028.pdf

Dmrecord 327291

Document Type Dissertation

Format application/pdf (imt)

Rights Chang, Ouliang

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