Particle Simulation of Lower Hybrid Waves and Electron-ion Hybrid
Instability
by
Lei Qi
A dissertation submitted to the Graduate Faculty of
Auburn University
in partial ful?llment of the
requirements for the Degree of
Doctor of Philosophy
Auburn, Alabama
May 03, 2014
Keywords: Particle simulation, Landau damping, Lower Hybrid Wave, parametric
instability, Electron-ion hybrid (EIH) instability
Copyright 2014 by Lei Qi
Approved by
Yu Lin, Chair, Professor of Physics
J.D. Perez, Professor of Physics
Thomas H. Pate, Professor of Mathematics
Edward Thomas, Jr., Professor of Physics
Kaijun Liu, Assistant Professor of Physics
Abstract
Lower hybrid wave (LHW) has been of great interest to laboratory plasma physics
for decades due to its important applications in particles heating and current drive
in plasmas devices. There are two fundamental characteristics of LHWs, Landau
damping and parametric instability (PI), both of which play key roles in particles
heating and current drive problems. Linear physics of LHWs has been studied in great
details in analytical theories, while nonlinear physics is usually too complicated to be
resolved analytically. Computer kinetic simulation technique has developed to be one
of the best tools for the investigation of kinetic physics, especially in the nonlinear
stage. Although a great amount of theoretical work has been done in the investigation
of LHWs, little particle simulation work can be found in published literatures. In
this thesis, an electrostatic gyro-kinetic electron and fully kinetic ion (GeFi) particle
simulation scheme is utilized to study the linear and nonlinear physics of LHWs. GeFi
model is particularly suitable for plasma dynamics with wave frequencies lower than
the electron gyrofrequency, and for problems in which the wave modes ranging from
Alfv?en waves to lower-hybrid/whistler waves that need to be handled on an equal
footing with realistic electron-to-ion mass ratio.
Firstly, Interactions of LHWs with both electrons and ions through nonlinear
Landau damping, as well as linear Landau damping, is investigated by utilization
of GeFi particle simulation in electrostatic limit. Landau damping, a wave-particle
interaction process, provides a way for LHWs to exchange energy and momentum
with electrons and ions, thus to heat plasmas and generate electric currents in the
plasmas. Unlike most other wave modes, LHWs can resonantly interact with both
electrons and ions, with the former being highly magnetized and the latter being
ii
nearly unmagnetized around lower hybrid frequency. Direct interactions of LHWs
with electrons and/or ions are investigated for cases with various k?/k, Ti/Te, and
wave amplitudes. Here, k is wave vector, k? is parallel (to static magnetic ?eld)
wave vector, Ti and Te are ion and electron temperatures, respectively. In the linear
electron Landau damping (ELD), real frequencies and damping rates obtained from
our kinetic simulations have excellent agreement with the analytical linear dispersion
relation. As wave amplitude increases, the nonlinear Landau e?ects are present, and a
transition from strong decay at smaller amplitudes to weak decay at larger amplitudes
is observed. In the nonlinear stage, LHWs in a long time evolution ?nally exhibit a
steady BGK mode, in which the wave amplitude is saturated above noise level. While
resonant electrons are trapped in the wave electric ?eld in the nonlinear ELD, resonant
ions are untrapped in LHWs time scales. Ion Landau damping is thus predominantly
in a linear fashion, leading to a wave saturation level signi?cantly lower than that
in the ELD. On long time scales, however, ions are still weakly trapped. Simulation
results show a coupling between LHW frequency and ion cyclotron frequency during
the long-time LHW evolution.
Secondly, our investigation extends to a transverse sheared ?ow driven instability,
electron-ion hybrid (EIH) instability, whose frequency is in the lower hybrid frequency
range. EIH instability is studied by the GeFi model in a magnetized plasma with a
localized electron cross-?eld ?ow. Macroscopic ?ows are commonly encountered in
various plasmas, such as plasmas in tokamak devices, laser-produced plasmas, Earth?s
magnetopause, plasma sheet boundary layer, and the Earth?s magnetotail. As a
benchmark, linear simulations of EIH are ?rstly performed in both slab geometry and
cylindrical geometry with kz = 0 in either uniform plasmas or nonuniform plasmas,
and the results are compared with linear theories in a slab geometry. Here for the
slab geometry, static magnetic ?eld is at z axis, and electron shear ?ow as a function
of x is put at y axis. And in the cylindrical geometry, static magnetic ?eld is also
iii
at z direction, while electron shear ?ow as a function of radial position r is in ?
(poloidal) direction. Linear eigen mode structures and growth rates of EIH instability
are calculated for various kyL, ?1 = V0E/L?e, and L/Ln. Here ky (kz) is the wave
vector in y (z) direction, L the shear length of electric ?eld, V0E the peak value of
the E ? B drift velocity, ?e the electron gyro frequency, Ln the scale of density
gradients. The results have very good agreement with the theoretical predications.
Nonlinear simulations are performed to investigate the nonlinear evolution of EIH
instabitly. It is found that the EIH instability nonlinearly evolves from a short wave
length (kx?i ? 12) mode to a long wave length (kx?i ? 3) mode with frequency
? ?LH. Simulation results under realistic plasma conditions of the Auburn Linear
Experiment for Instability Studies (ALEXIS) device are discussed and compared with
ALEXIS experimental results as well.
Finally, parametric instability of LHWs is investigated using the GeFi model.
Parametric instability is a nonlinear process that involves wave-wave interactions, and
is of great interest. In the propagation of LHWs, a pump LHW (?0,k0) decays into
a low frequency wave mode (?,k) and two high frequency sidebands (???0,k?k0),
where ? represent the wave frequency, k represent the wave vector, and the subscript
?0? indicates the pump wave. Two di?erent cases with parametric instability are
discussed in great details. The parametric decay process is found to occur very fast
within several lower hybrid wave periods. Growth rates of the excited modes are es-
timated as well, and results are compared with the analytical theory. The simulation
shows that the parametric instability process is complicated due to multiple decay
channels. These multiple parametric instability processes usually occur simultane-
ously. The corresponding electron and ion particle distributions are investigated in
the decay process. Finally, a discussion is presented for the future study of electron
and ion nonlinear physics of the PI instability.
iv
Acknowledgments
I would like to express my great gratitudes to my supervisor Dr. Yu Lin for
her patient guidance in my research work. Without her, I would never be able to
complete my research program. I am also greatly appreciative of helpful discussions
with my collaborator Dr. Xueyi Wang.
v
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii
1 An Introduction of Lower Hybrid Waves and Their Landau Damping and
Parametric Instability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Introduction of Lower Hybrid Waves . . . . . . . . . . . . . . . . . . 1
1.1.1 Electromagnetic dispersion relation in a cold plasma . . . . . . 2
1.1.2 Electrostatic dispersion relation in a cold plasma . . . . . . . 4
1.1.3 LHW dispersion relations in warm plasmas . . . . . . . . . . . 7
1.2 Landau damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.1 Liner Landau damping . . . . . . . . . . . . . . . . . . . . . . 9
1.2.2 Nonlinear Landau damping . . . . . . . . . . . . . . . . . . . 10
1.3 Landau damping of electrostatic LHWs . . . . . . . . . . . . . . . . . 14
1.3.1 Coe?cient of the linear electron Landau damping . . . . . . . 14
1.3.2 Physical mechanisms of the Landau damping of LHWs . . . . 16
1.4 Parametric instability . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.4.1 The PI dispersion relation near lower hybrid frequency . . . . 17
1.4.2 Decay channels . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2 GeFi Simulation Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.1 GeFi Scheme Algebras . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2 GeFi model in the electrostatic limit . . . . . . . . . . . . . . . . . . 37
vi
2.3 Benchmark of the GeFi Scheme . . . . . . . . . . . . . . . . . . . . . 38
2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3 GeFi Simulation of Landau Damping of Lower Hybrid Waves . . . . . . . 41
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3 Simulation of linear ELD of LHWs . . . . . . . . . . . . . . . . . . . 42
3.4 Nonlinear Landau Damping of Lower Hybrid Waves . . . . . . . . . . 44
3.4.1 Case 1, Electron Landau Damping of LHWs . . . . . . . . . . 46
3.4.2 Case 2, Ion Landau Damping of LHWs . . . . . . . . . . . . . 48
3.4.3 case 3, Electron And Ion Landau Damping of LHWs . . . . . 51
3.5 Nonlinear Damping Rates and Saturation of LHWs . . . . . . . . . . 54
3.5.1 Landau damping rate of LHWs . . . . . . . . . . . . . . . . . 54
3.5.2 Saturated Electric Field, Es . . . . . . . . . . . . . . . . . . . 57
3.6 Calculation of Driven Current by LHWs . . . . . . . . . . . . . . . . 59
3.6.1 Fast Electrons Drive . . . . . . . . . . . . . . . . . . . . . . . 60
3.6.2 Electric Currents obtained from GeFi Simulation of LHWs . . 61
3.7 Ions Phase Bunching . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4 GeFi Simulation of Electron-ion Hybrid instability . . . . . . . . . . . . . 70
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.2 Linear theory of EIH instability . . . . . . . . . . . . . . . . . . . . . 71
4.2.1 General physical model . . . . . . . . . . . . . . . . . . . . . . 72
4.2.2 EIH instability dispersion equations . . . . . . . . . . . . . . . 76
4.2.3 Numerical methods for solving the dispersion equation . . . . 77
4.3 Theoretical numerical solutions of the EIH instability . . . . . . . . . 79
4.3.1 Numerical Investigation of EIH Instability in Uniform Plasmas 79
4.3.2 Numerical investigation of EIH instability in nonuniform plasma 81
vii
4.4 Linear GeFi simulation of the EIH instability . . . . . . . . . . . . . 84
4.4.1 Simulation model . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.4.2 EIH instability in a slab geometry and uniform plasma with
kz = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.4.3 EIH instability in a slab geometry and nonuniform plasma with
kz = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.4.4 EIH instability in cylindrical geometry and uniform plasma
with kz = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.4.5 EIH instability in cylindrical geometry and nonuniform plasma
with kz = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.4.6 EIH instability with kz ?= 0 . . . . . . . . . . . . . . . . . . . . 95
4.5 Nonlinear evolution of EIH instability . . . . . . . . . . . . . . . . . . 97
4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5 Nonlinear Parametric Decay Instability of Lower Hybrid Wave . . . . . . 104
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.2 GeFi Simulation of Parametric Instability of LHWs . . . . . . . . . . 105
5.2.1 Single Mode simulation of Pump Wave . . . . . . . . . . . . . 107
5.2.2 GeFi simulation of PI of LHWs . . . . . . . . . . . . . . . . . 108
5.2.3 Particles distribution in PI simulation . . . . . . . . . . . . . . 120
5.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
5.4 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
viii
List of Figures
1.1 n2? versus N for ?xed n?, NS (NF) denotes the cuto? density to the slow
(fast) wave.[Figure is from: Paul Bonoli, IEEE Transactions on plasma
science, Vol. PS-12, NO. 2, June 1984.] . . . . . . . . . . . . . . . . . . . 5
1.2 Phase trajectories of the resonant electrons . . . . . . . . . . . . . . . . 12
1.3 Manfredi?s 3 runs of the time evolution of the amplitude of the electric ?eld 13
1.4 Manfredi?s shaded plot of the distribution function in the resonant region 14
2.1 Comparison between the dispersion relations obtained from the kinetic
GeFi simulation and the corresponding analytical linear dispersion rela-
tions for various k?/k?. Top plot is ?Bz in the fast magnetosonic/whistler
branch. Bottom plot is ?By in the shear Alfv?en/kinetic Alfv?en mode
branch. The dashed line shows the analytical dispersion relation of the
MHD shear Alfv?en mode. . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1 Perturbed electrostatic potential ?1 (in logarithm scale) vs. ?it for a case
of linear electron Landa damping of LHWs . . . . . . . . . . . . . . . . 43
3.2 Diamonds show the real frequency and the linear ELD rate as a function
of the wave number for LHWs with B0x = 0.0659 and B0y = 0.999, which
corresponds to a ?xed k?/k = 0.066, and Ti/Te = 1.0. The solid lines are
based on the analytical theory. . . . . . . . . . . . . . . . . . . . . . . . 45
ix
3.3 Time evolution of the spatial Fourier mode of the electric ?eld (in the
natural logarithm scale) for case 1 with k?/k = 0.066 and Ti/Te = 1.0,
and k?e = 0.2255, E0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.4 (a) Contour plots of the electron distributions in the particle phase space
(x,ve?) and (b) the corresponding parallel velocity distribution functions
averagedoverx, shownwiththeblacksolidcurves, attimes?it = 0,0.22,2.41
and 8.78 obtained from Case 1. All the distribution functions are plotted
in the logarithm scales The red solid lines show the ion distribution func-
tions, and the black dashed lines show the electron resonant phase velocity
ver = 3.785Vte based on the theoretical prediction. The ion velocities are
normalized to the ion thermal speed. . . . . . . . . . . . . . . . . . . . . 49
3.5 Time evolution of the spatial Fourier mode of the electric ?eld (in the
natural logarithm scale) for case 2 with k?/k = 0.001, Ti/Te = 1.0, k?e =
0.2255, and E0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.6 Time (normalized to 1/?i) evolution of the electron phase-space distribu-
tion in case 2: contour plots of the (a) electron and (b) ion phase-space
distributions and (c) the corresponding electron (black solid lines) and ion
(red solid lines) distribution functions. The ion velocities are normalized
to the ion thermal speed. The black (red) dash line shows the theoretical
electron (ion) resonant phase-velocity vir = 3.264Vti. . . . . . . . . . . . 52
3.7 Time evolution of the spatial Fourier mode of the electric ?eld (in the
natural logarithm scale) for case 3 with k?/k = 0.0404, Ti/Te = 4.0,
k?e = 0.2255, and E0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.8 Resulting particle distributions obtained from case 3 in the same format
as Fig. 3.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
x
3.9 Damping rate, ?/?i, as a function of the initial wave amplitude E0 in the
cases with k?e = 0.2255 and Ti/Te = 1, for k||/k = 0.0904, 0.066, 0.0404,
and 0.001. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.10 Damping rate vs. E0 in the cases with k?e = 0.2255 and k?/k = 0.0404,
for Ti/Te = 4.0, 1.0 and 0.33. . . . . . . . . . . . . . . . . . . . . . . . . 58
3.11 Saturated electric ?eld Es/E0 vs. Ti/Te in the cases with k?e = 0.2255
and k?/k = 0.066, 0.0404,0.001 . . . . . . . . . . . . . . . . . . . . . . . 59
3.12 The plateau in the parallel velocity distribution function in the electron
Landau damping of LHWs (V par is the parallel velocity, v?, normalized
in Vte) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.13 Time evolution of J? (normalized to en0Vte) obtained from case 1. . . . 64
3.14 Resulting parallel currents as a function of the initial LHW amplitude E0. 65
3.15 Scatter plots of ions distribution for selected particles in velocity phase
space (Vix, Viz) at ?it = 0.0, 0.8075, 1.6231, 6.0275, 6.3537, 6.8431. . . . 66
3.16 Ions velocity distribution function in terms of Vix at ?it = 5.7012, 6.0275,
6.1906, 6.8431. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.1 A typical eigenfunction electrostatic potential (?(x,?)) of EIH instability.
Here ? = 1.0, ?1 = 0.3 and kyL = 0.5. For this case, the frequency ? =
(4.35+2.84i)?LH. The black line shows the real part of the eigenfunction,
and the red line represents the imaginary part of the eigenfunction. . . . 80
4.2 The real parts (dashed lines) and imaginary parts (solid lines) of the eigen-
values (?) of Eq. 4.17 as a function ofkyLfor? = 0.5 (red lines), 1.0 (black
lines) and 5.0 (blue lines). Here mi/me = 1836 and ?1 = 0.3. . . . . . . . 81
xi
4.3 The real parts (dashed lines) and imaginary parts (solid lines) of the eigen-
values (?) of Eq. 4.17 as a function of ?1 for ? = 0.5 (red lines), 1.0 (black
lines). Here mi/me = 1836 and kyL = 0.5. . . . . . . . . . . . . . . . . . 82
4.4 The real parts (black solid line) and imaginary parts (red dashed line) of
the eigenvalue (?, ?) of Eq. 4.18 as a function of kyL for ? = 1.0, ?1 = 0.1
and L/Ln = 0.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.5 The real parts (black solid lines) and imaginary parts (red dashed line)
of the eigenvalue (?, ?) of Eq. 4.18 as a function of L/Ln for ? = 1.0,
kyL = 0.4 and ?1 = 0.1, L is ?xed as 5.0?e. . . . . . . . . . . . . . . . . 83
4.6 Plot of the eigenfunction of the perturbed electrostatic potential. Solid
black (red) line is the real (imaginary) eigenfunction from numerically
solving the dispersion equation. Dashed black (red) line is the real (imag-
inary) eigenfunction from the GeFi simulation. . . . . . . . . . . . . . . 86
4.7 Contour plot of perturbed physical quantities in (x,y) real space. . . . . 87
4.8 Plot of real frequency (?) and growth rate (?) normalized by lower hybrid
frequency (?LH), as a function of kyL. Dashed black line is the real fre-
quency and solid red line is the growth rate, both of which are from the
theory. The black dots (red dots) are real frequency (growth rate) from
the GeFi simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.9 Plot of real frequency (?) and growth rate (?) normalized to lower hybrid
frequency (?LH), as a function of ?1. Dashed black line is the real fre-
quency and solid red line is the growth rate, both of which are from the
theory. The black dots (red dots) are real frequency (growth rate) from
the GeFi simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
xii
4.10 Plot of real frequency and growth rate as a function of kyL. . . . . . . . 91
4.11 Plot of real frequency and growth rate as a function of the ratio L/Ln
with L = 0.35cm and kyL = 0.22 being ?xed. . . . . . . . . . . . . . . . 91
4.12 Contour plot of perturbed physical quantities in (r,?) real space . . . . . 92
4.13 Plot of real frequency (?) and growth rate (?) normalized to lower hybrid
frequency (?LH), as a function of k L. . . . . . . . . . . . . . . . . . . . 93
4.14 Real frequency and growth rate as a function of k L . . . . . . . . . . . 94
4.15 Real frequency and growth rate as a function of ?1 . . . . . . . . . . . . 95
4.16 Growth rates vs kzL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.17 A measurement from the ALEXIS experiment, in which no instability is
found. The spikes are the noise in the circuit. . . . . . . . . . . . . . . . 97
4.18 Frequencies and growth rates vs k L for kzL = 0.005 and kzL = 0.000. . 98
4.19 Contour plot of electrostatic potential in (?,x) space in a nonlinear sim-
ulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.20 Contour plot of electrostatic potential in (?,x) space in a linear simulation. 99
4.21 Plot of electrostatic potential as a function of ?. . . . . . . . . . . . . . . 100
4.22 Contour plots of physical quantities in real space from a nonlinear simu-
lation with ? = 1.0,?1 = 0.3,kyL = 0.35 at ?it = 0.0599. . . . . . . . . 100
4.23 Contour plots of physical quantities in real space from a nonlinear simu-
lation with ? = 1.0,?1 = 0.3,kyL = 0.35 at ?it = 0.8497. . . . . . . . . . 101
xiii
4.24 Nonlinear EIH frequency as a function of kyL for cases with ? = 1.0,?1 =
0.1. Black solid line is nonlinear GeFi particle simulation results, and red
dashed line is for linear theory. . . . . . . . . . . . . . . . . . . . . . . . 102
5.1 Time evolution of the Fourier electric ?eld of the pump wave in a single
mode simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.2 Contour plot of the electric ?eld in the (kx,kz) space with the mode num-
ber in y direction my being 0. . . . . . . . . . . . . . . . . . . . . . . . 109
5.3 Contour plot of the electric ?eld in the (kx,kz) space with the mode num-
ber in y direction my being 1. Wave modes with (0,1,0), (4,1,4) and
(?4,1,?4) are found to have the largest amplitudes of electric ?eld. . . . 110
5.4 Time evolution of Fourier electric ?eld of the 3 wave modes: pump wave
(4,0,4) (black line), wave mode (0,1,0) (red line) and wave mode (4,1,4)
(green line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.5 Frequencies of pump wave (4,0,4), wave mode (0,1,0) and wave mode
(4,1,4) from top to bottom. . . . . . . . . . . . . . . . . . . . . . . . . 112
5.6 Time evolution of Fourier electric ?eld of the 3 wave modes: pump wave
(4,0,4)(blackline), wavemode(0,1,0)(redline)andwavemode(?4,1,?4)
(green line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.7 Frequencies of pump wave (4,0,4), wave mode (0,1,0) and wave mode
(?4,1,?4) from top to bottom. . . . . . . . . . . . . . . . . . . . . . . . 115
5.8 Contour plot of the electric ?eld in (kx,kz) space with mode number in y
direction my being 2. Wave modes with (0,2,0), (4,2,4) and (?4,2,?4)
are found to have the largest amplitudes of electric ?eld. . . . . . . . . . 116
xiv
5.9 Time evolution of Fourier electric ?eld of the 4 wave modes: pump wave
(4,0,4) (black line), wave mode (0,2,0) (red line), wave mode (4,2,4)
(yellow line) and wave mode (?4,2,?4) (green line). . . . . . . . . . . . 116
5.10 Frequencies of pump wave (4,0,4), wave mode (0,2,0), wave mode (4,2,4)
and wave mode (?4,2,?4) from top to bottom. . . . . . . . . . . . . . . 117
5.11 Contour plot of the electric ?eld in (kx,kz) space with mode number in y
direction my = 3. Wave modes with (0,3,1), (4,3,5) and (?4,3,?3) are
found to have the largest amplitudes of electric ?eld. . . . . . . . . . . . 118
5.12 Contour plot of the electric ?eld in (kx,kz) space with mode number in
y direction my = 3 at ?it = 0.58824. Wave modes (1,3,0) and (?1,3,0)
are found to have large amplitudes of electric ?eld. . . . . . . . . . . . . 118
5.13 Contour plot of the electric ?eld in (kx,kz) space with mode number in
y direction my = 3 at ?it = 0.54466. Wave modes with (5,3,3) and
(?5,3,?3) are found to have large amplitudes of electric ?eld. . . . . . . 119
5.14 Time evolution of Fourier electric ?eld of the 4 wave modes: pump wave
(4,0,4) (black line), wave mode (0,3,1) (red line), wave mode (4,3,5)
(yellow line) and wave mode (?4,3,?3) (green line). . . . . . . . . . . . 119
5.15 Frequencies of pump wave (4,0,4), wave mode (0,3,1), wave mode (4,3,5)
and wave mode (?4,3,?3) from top to bottom. . . . . . . . . . . . . . . 120
5.16 Time evolution of Fourier electric ?eld of 3 wave modes: pump wave
(4,0,4) (black line), wave mode (1,3,0) (red line), and wave mode (5,3,4)
(green line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
xv
5.17 Time evolution of Fourier electric ?eld of 3 wave modes: pump wave
(4,0,4) (black line), wave mode (?1,3,0) (red line), and wave mode
(?5,3,?4) (green line). . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
5.18 Frequencies of pump wave (4,0,4), wave mode (1,3,0) and wave mode
(5,3,4) from top to bottom. . . . . . . . . . . . . . . . . . . . . . . . . . 122
5.19 Frequencies of pump wave (4,0,4), wave mode (?1,3,0) and wave mode
(?5,3,?4) from top to bottom. . . . . . . . . . . . . . . . . . . . . . . . 122
5.20 Time evolution of electron and ion distributions in Case 1 for ?it = 0.0218
(top), ?it = 0.5229 (middle), ?it = 0.6972 (bottom). Column (a) is
electron distribution in phase space (Ve?,x). (b) is ion distribution in
phase space (Vi,x). (c) are plots of electron (solid lines) and ion (red
lines) distribution functions. . . . . . . . . . . . . . . . . . . . . . . . . . 123
xvi
List of Tables
xvii
Chapter 1
An Introduction of Lower Hybrid Waves and Their Landau Damping and
Parametric Instability
The existence of a damping mechanism by which plasma particles absorb wave
energy even in a collisionless plasma was found by L.D. Landau[1], under the condition
that the plasma is not cold and the velocity distribution is of ?nite extent. Energy-
exchange processes between particles and waves play important roles in plasma heat-
ing by waves and in the mechanism of instabilities. Previous investigation of Landau
damping mechanism mainly used Langmuir waves, typical electrostatic longitudinal
plasma waves. Electrostatic lower hybrid wave (LHW) is a preferable source to heat
both electrons and ions to the thermonuclear temperature and generate electric cur-
rents in plasmas under Landau damping. Thus the Landau damping of LHWs is
fundamentally important. In this chapter, I will introduce and give an overview of
linear dispersion relation of LHWs, Landau damping theory, theory of Landau damp-
ing of LHWs, and theory of parametric instability (PI).
1.1 Introduction of Lower Hybrid Waves
Lower hybrid wave is of great interest for decades due to its important ap-
plications in the magnetic controlled fusion devices. Electrostatic LHW is a po-
tential source to heat both electrons and ions to thermonuclear temperature and
generate electric currents[50] in plasma fusion devices. LHWs have been inves-
tigated extensively in linear theories[16, 17, 18], nonlinear parametric instability
theories[19, 20, 21], and experiments[22]. Some experiments were performed to heat
electrons in plasmas by the direct electron Landau interactions in the lower hybrid
1
waves[23, 24, 25], and a su?ciently high increase of the electron temperature was
obtained. In this section, an introduction of the lower hybrid waves are presented, in-
cluding electromagnetic dispersion relation in a cold plasma, electrostatic dispersion
relation in a cold plasma, and dispersion relations in warm plasmas.
1.1.1 Electromagnetic dispersion relation in a cold plasma
Consider a plasma with density gradient in the x direction N(x), an external
static magnetic ?eld in the z direction B0 = B0ez, unperturbed electric ?eld E0 = 0,
and unperturbed electron and ion ?uid velocities ve0 = 0 and vi0 = 0. The perturbed
electric ?eld E, magnetic ?eld B, and the current density J are then described by
Maxwell?s equations
??B = ?0?0?E?t +?0J (1.1)
??E = ??B?t . (1.2)
Assuming that the wave ?elds vary as a wave, i.e., exp[i(k?x??t)] and the WKB
approximation (|k|?|?/?x|). Eqs. 1.1 and 1.2 can be analyzed by Fourier methods
in space and time
?k?k?E = (?/c)2E+i??0J. (1.3)
Where J = Ne(vi ? ve), with the electron (ion) velocity ve (vi) governed by the
equations of motion
me?ve?t = ?e(E+ve ?B0) (1.4)
mi?vi?t = e(E+vi ?B0) (1.5)
Combining Eqs. 1.3, 1.4 and 1.5, it yields the following equations
k?k?E+ (?/c)2???E = 0, (1.6)
2
?? =
?
??
??
??
?
?? ?i?xy 0
i?xy ?? 0
0 0 ??
?
??
??
??
?
. (1.7)
?? is so called the cold plasma dielectric tensor. Taking the approximations that the
LHW?s frequency ?i ? ? ? ?e, where ?e = eB0/me (?i = eB0/mi) is the electron
(ion) cyclotron frequency, one can obtain the elements of the tensor
?? = 1 + (?pe/?e)2 ?(?pi/?)2, (1.8)
?? = 1?(?pe/?)2 ?(?pi/?)2, (1.9)
?xy = ?2pe/(??e). (1.10)
Note that ?pe =
?
Ne2/?0me (?pi =
?
Ne2/?0mi) is the electron (ion) plasma fre-
quency. Eq. 1.6 can be written in another form as
?D?E = 0, (1.11)
?D = k?k??I+ (?/c)2??. (1.12)
Let the determinant of ?D vanish to obtain the nontrivial solutions to Eq. 1.11. This
results in the LHW electromagnetic dispersion relation in a cold plasma
D0(x,k,?) = det(?D) = P4n4? +P2n2? +P0 = 0 (1.13)
P0 = ??[(n2? ???)2 ??2xy]
P2 = (?? +??)(n2? ???) +?2xy
P4 = ??,
3
where n? = k?c/?, n? = k?c/?, k? (k?) is the component of the wave vector parallel
(perpendicular) to the external static magnetic ?eld. Here and through this thesis, ?
(?) denotes the component parallel (perpendicular) to the static magnetic ?eld. Eq.
1.13 indicates two wave modes in terms of n?. By solving n2? in Eq. 1.13, one can
obtain
n2? = 12P
4
(?P2 ??1=2), (1.14)
with ? = (P2)2 ? 4P0P4. The positive sign in Eq. 1.14 corresponds to the ?slow?
wave branch (small?/k?), and the minus sign corresponds to the ?fast? wave branch.
The slow wave branch has a cold plasma resonance condition P4 ? 0, from which the
lower hybrid wave frequency can be de?ned as ?LH = ?pi/
?
1 +?2pe/?2e. The slow
wave branch of the dispersion relation, is typically associated with the lower hybrid
wave.
Fig 1.1 shows the plots of n2? versus density N(x) for three di?erent values of
n?. Note that there is a critical value n? = nc. Fig 1.1(a) shows there exists mode
conversions between the slow wave mode and the fast wave mode with the condition
n? <nc. A wave launched at the slow wave cuto? NS will mode convert into the fast
wave mode at a density denoted asNT1. Withn? >nc, as plotted in Fig 1.1(c), a wave
launched at the slow wave cuto? will propagate to its resonance without converting
to the fast wave. n? = nc is the lowest value for the slow wave mode to propagate
from its cuto? to the maximum plasma density in the presence of density gradient
without mode conversion to the fast wave, as plotted in Fig 1.1(b).
1.1.2 Electrostatic dispersion relation in a cold plasma
The wave polarization is essentially electrostatic (kc/? ? 1) when the LHW
propagates at densities much greater than cuto? densities NS and NF. Thus, to
consider propagations of the electrostatic LHW nearly perpendicular to the magnetic
4
Figure 1.1: n2? versus N for ?xed n?, NS (NF) denotes the cuto? density to the slow
(fast) wave.[Figure is from: Paul Bonoli, IEEE Transactions on plasma science, Vol.
PS-12, NO. 2, June 1984.]
5
?eld, the dispersion relation can be obtained from Eq. 1.13
D0(x,k,?) = k2??? +k2??? = 0. (1.15)
Given the condition k? ?k?, the above dispersion relation can be reduced to
?2 = ?2LH(1 + mim
e
k2?
k2). (1.16)
Let us now understand the main features of Eqs. 1.15 and 1.16. Evaluation of
??/?k? shows that the wave is backward in the perpendicular direction (to static
magnetic ?eld) [i.e., (??/?k?)(?/k?) < 0]. The dispersion relation also shows that,
k? increases as the wave propagates into regions of increasing plasma density, when
there exists a density gradient in the ? direction. In particular, for a ?xed k?, k?
becomes very large as the wave reaches regions, where the local lower hybrid frequency
?LH is close to the wave frequency, ?. Mode conversion can occur near the lower
hybrid layer.[26]
Resonance cones
Another interesting characteristic of the propagation of the electrostatic LHW
is that point sources excite LHWs to propagate in resonance cones. To consider that
a point source excites an electrostatic disturbance in a plasma with frequency and
wave number (?,k), satisfying the dispersion relation Eq. 1.15. For simplicity, take
k = k?ex +k?ez and E = ???, where ?? exp[i(k?x+k?z??t)]. Then the Fourier
transformation in Eq. 1.15 my be inverted to give
?2?
?x2 = ?
??
??
?2?
?z2. (1.17)
6
This equation is similar to the typical wave di?erential equation
?2?
?t2 = C
2?
2?
?z2. (1.18)
The solutions of Eq. 1.18 propagate along the characteristics z?Ct. In analogous
fashion, the solutions of Eq. 1.17 propagate along characteristics z?(???/??)1=2x.
These are the resonance cones, in which there exists the propagation of a singular
electrostatic disturbance excited by a point source and obeying the LHW dispersion
relation.
1.1.3 LHW dispersion relations in warm plasmas
Both the electromangetic and electrostatic LHW dispersion relations in a cold
plasma were discussed in the previous introductions. Here, to include the warm
plasma e?ects, both the electromagnetic and electrostatic LHW dispersion relations
will be brie?y reviewed. It is assumed that ?pe ? ?e ? ? ? ?i, k2V2te/?2e ? 1
with Vte =
?
?Te/me the electron thermal speed, and the ions are singled charged,
me ?mi. ? is de?ned to be the angle between k and B, thus cos(?) = k?/k ?me/mi
In the dispersion relation, terms are kept to the ?rst order of k2V2te/?2e and to all
orders of the electromagnetic term ?2pe/k2c2. The LHW dispersion relation with both
electromagnetic and warm e?ects included is found to be[27]
?2 = ?
2
LH
1 +?2pe/k2c2[1 +
(mi/me)cos2(?)
1 +?2pe/k2c2 +W
k2V2te
?2e ], (1.19)
W = 3TiT
e
(1 + ?
2
pe
k2c2) +
?2pe
2k2c2 +
9
2 ?
15 + 21?2pe/k2c2
4(1 +?2pe/k2c2)2
? [3 ?
2
pe
k2c2 +
1?6?2pe/k2c2
4(1 +?2pe/k2c2)2]
mi
me cos
2(?)
7
+ 3(mi/me)cos
2(?) +?2
pe/k
2c2 ?Ti/Te
1 +?2pe/k2c2 + (mi/me)cos2(?) (1 +
?2pe
k2c2)
mi
me cos
2(?). (1.20)
In the limit that ?2pe/k2c2 tends to zero and (mi/me)cos2(?) is small, the LHWs
become longitudinal and nearly electrostatic, and W = 3Ti/Te +3/4 to lowest order.
The dispersion relation for electrostatic LHW in a warm plasma becomes
?2 = ?2LH(1 + mem
e
cos2(?) + 3(TiT
e
+ 14)k
2V2
te
?2e ). (1.21)
To include the electromagnetic e?ects, in the limit that ?2pe/k2c2 and (mi/me)cos2(?)
are both small, the dispersion relation for the electromagnetic LHW in a warm plasma
is found to be
?2 = ?2LH(1 + mem
e
cos2(?)? ?
2
pe
2k2c2 + 3(
Ti
Te +
1
4)
k2V2te
?2e ). (1.22)
There are several more expressions of the approximate dispersion relation of the
LHW. If the plasma is treated as being cold and EM e?ects are included, the resulting
dispersion relation is[28]
?2
?2LH =
1
1 +?2pe/k2c2(1 +
mi
me)
cos2(?)
1 +?2pe/k2c2. (1.23)
This expression makes no assumption on the size of ?2pe/k2c2.
Another LHW dispersion relation, adopted by Bingham et al.[29], includes both
warm plasma and EM e?ects, but assumes that ?2pe/k2c2 is small and includes terms
only to ?rst order in (mi/me)cos2(?), ?2pe/k2c2 and k2V2te/?2e.
?
?LH = 1 +
mi
2me cos
2(?)? ?
2
pe
k2c2 + (
3Ti
2Te + 1)
k2V2te
?2e . (1.24)
8
Note that the last term in the above dispersion relation is not equal to half of the last
term of Eq. 1.21 (or Eq. 1.22), although in the limit of ?2pe/k2c2 ? 0 they should be
equal.
1.2 Landau damping
Landau damping is a fundamental plasma process in which small perturbations
in a uniform, Maxwellian, electrostatic plasma are exponentially damped, even when
no dissipative terms are present. Landau damping process can be categorized as
linear Landau damping and nonlinear Landau damping according to di?erent time
scales, which will be discussed in details later. The nonlinear Landau damping has
distinctive characteristics that are quite di?erent from the linear Landau damping.
The linear Landau damping of Langmuir waves was extensively investigated in 50?s
and 60?s and well understood. The nonlinear Landau damping, however, was not well
understood. Recent Landau damping experiments[7] and plasma simulations[5, 6]
revealed new important physics of Landau damping of Langmuir waves for long time
evolution.
1.2.1 Liner Landau damping
The Landau damping process was found by Landau in 1946. Since then, linear
Landau damping has been extensively con?rmed in both experiments[8] and com-
puter simulations[9]. The Landau damping, especially the linear Landau damping is
a standard topic in most plasma physics textbooks (e.g., [10, 11]), mainly for Lang-
muir wave. The linear Landau theory leads to the damping coe?cient [1, 2] for the
collisionless damping of longitudinal electron plasma oscillations
?L
? =
?
2
?2pe
k2
?f
?v
  
 v=!=k, (1.25)
9
where ? is real frequency of the electron plasma wave, ?pe is the electron plasma fre-
quency, k is the wave number, and f is the velocity distribution function of electrons.
The physical interpretation of this collisionless damping is that the electron plasma
wave resonates with the electrons, which possess phase velocity v = ?/k, and loses
energy to the electrons. However, Landau?s treatment of the problem is rigorous, but
strictly linear.
The linear theory will break down after a time ? =
?
me/eEk, where me is
the electron mass, e is electron charge, and E is the amplitude of the electric ?eld.
In other words, if the linear damping rate ?L is comparable to 2?/?, the initial
decay of the electron plasma wave will soon turn into nonlinear Landau damping,
characterized as a nonlinear oscillation due to particle trapping with the oscillation
amplitude somewhat lower than the initial value[12, 13].
1.2.2 Nonlinear Landau damping
The linear Landau damping was researched thoroughly in theory, experiments
and kinetic simulations. Nevertheless, the nonlinear Landau damping is still poorly
understood. There are no exact analytical solutions of the nonlinear Landau damp-
ing, but only some approximate ones[2]. The early theory, presented by O?Neil, pre-
dicted phase mixing and amplitude oscillations for the electric ?eld, which have indeed
been demonstrated in simulations and experiments[14] of Langmuir waves. However,
O?neil?s treatment becomes invalid for a long wave time evolution. It is now found[10]
that nonlinear plasma waves undergo a few amplitude oscillations and eventually ap-
proach a Bernstein-Greene-Kruskal (BGK) steady state[4] instead of decaying away.
This conclusion could also be drawn from the recent Landau damping experiments[7]
and plasma simulations for a long wave time evolution[5, 6] of Langmuir waves.
O?Neil?s theory[2] provided an exact solution of the Vlasov equation for the reso-
nantelectrons using Jacobi elliptic functions, and the damping coe?cientas a function
10
of time is given by
?(t) = ?L
??
n=0
64
?
? 1
0
d?{ 2n?
2 sin(  nt
 F )
?5F2(1 +q2n)(1 +q?2n) +
(2n+ 1)?2?sin((2n+1) t2F )
F2(1 +q2n+1)(1 +q?2n?1)},
(1.26)
? t
0
?(t?)d(t?) = ?L?
??
n=0
64
?
? 1
0
d?{ 2?(1?cos(
 nt
 F ))
?4F(1 +q2n)(1 +q?2n)+
2??(1?cos((2n+1) t2F ))
F(1 +q2n+1)(1 +q?2n?1)}.
(1.27)
Where ?2 = 2eE/(kW+eE) with W being the conserved energy, F = F(?, 12?), with
? taking the sign of 12?, is the elliptic integral of the ?rst kind , and q = exp(?F?/F).
By taking into account the ? dependence of F and ?F, it can be shown that the
integrals over sin[(2n+1)?t/2F?], sin[?nt/?F?], cos[(2n+1)?t/2F?], cos[?nt/?F?],
phase mix to zero as t/? approaches in?nity. Thus, Eqs. (1.26) and (1.27) become
?(t = ?) = 0, (1.28)
? ?
0
?(t)dt = ?L?
??
n=1
64
?
? 1
0
d?{ 2??4F(1 +q2n)(1 +q?2n) + 2??F(1 +q2n?1)(1 +q?2n+1)}.
(1.29)
After the series in Eq. 1.29 are summed
? ?
0
?(t)dt = ?L?64?
? 1
0
d?{ 1?4[E? ? ?4F] + ??[E + (?2 ?1)F]} = O(?L?). (1.30)
O?Neil?s treatment, as presented above, is valid only if the nonlinear e?ects prevent the
amplitude from decaying by a signi?cant amount. From Eq. 1.30, the above condition
is satis?ed when?L? ? 1. Meanwhile, Eq. 1.29 shows that ast/? approaches in?nity,
?(t) phase mixes to zero. The wave evolves into a BGK[4] equilibrium, in which the
wave amplitude oscillates with a ?nite constant value.
Note that ?(t) is the damping rate, thus the value of ?L? is also the critical
value of the linear Landau damping and the nonlinear Landau damping. In general,
11
?L? ? 1, we consider the damping process mainly as nonlinear Landau damping,
otherwise, the process is treated as linear Landau damping.
The physical interpretation of the nonlinear Landau damping process lies basi-
cally on the particles trapping mechanism. As the Langmuir wave propagates in the
plasma, some electrons can be trapped in the potential well of the wave?s electric ?eld.
As the trapped electrons bounce back and forth in the potential well with a period of
order ?, the trapped electrons keep exchanging energy with the the wave. The phase
mixing of the trapped electrons causes the wave amplitude to reach an asymptotic
constant, forming a stable equilibrium with an undamped nonlinear plasma wave, i.e.,
the BGK equilibrium. Fig 1.2 shows the physical picture of the electron distribution
function f in the phase space (x,v). There is a vortex structure in the phase space of
the distribution, where is dominated by the trapped electrons with cyclic trajectory.
Figure 1.2: Phase trajectories of the resonant electrons
Recently, Giovanni Manfredi[5] presented the complete behavior of nonlinear
Landau damping of Langmuir waves in a particle simulation. The simulation work
showed explicitly the characteristics of the nonlinear Landau damping in the electron
plasma wave. Fig 1.3 shows the evolution of the amplitude of the electric ?eld in
Manfredi?s 3 runs. It is seen from the ?gure, as discussed in O?Neil?s theory, the wave
12
amplitude decreases within the ?rst few oscillations, but eventually settles down to
a steady oscillation with an undamped constant amplitude. According to O?Neil?s
theory, some electrons can be trapped and there will be a vortex structure in the phase
space of the distribution function at the resonance region, as shown by Manfredi (Fig.
1.4).
Figure 1.3: Manfredi?s 3 runs of the time evolution of the amplitude of the electric
?eld
13
Figure 1.4: Manfredi?s shaded plot of the distribution function in the resonant region
1.3 Landau damping of electrostatic LHWs
In this section the theory of Landau damping in an electrostatic LHW in the
cold plasma is introduced. The coe?cient of the electron Landau damping (ELD) of
the electrostatic LHWs will be derived. Di?erent mechanisms between the electron
Landau damping and the ion Landau damping of electrostatic LHWs will be discussed
as well as the ion Landau damping rate.
1.3.1 Coe?cient of the linear electron Landau damping
In Eq. 1.15, consider that the wave vector k is ?xed, if D0(x,k,?) has an
imaginary part due to some damping mechanism, then ? must become complex
? = ?r +i?. (1.31)
If the damping is weak, then ?/?r ? 1 and the dispersion relation may be written as
D0(x,k,?) ?D0r(x,k,?r +i?) +iD0i(x,k,?r) = 0. (1.32)
14
Taking real and imaginary parts of Eq. 1.32, one can obtain the general formula to
calculate the coe?cient of the linear ELD of LHWs
D0r(x,k,?r) ? 0
? ? ?D0i(x,k,?r)?D
0r/??r
. (1.33)
In the presence of the linear electron Landau damping, the dispersion relation of
the LHWs can be written as[16]
D0 = k
2
?
k2?? +
k2?
k2?? +
i??
k2?2D
??
2k?Vte exp(?
?2
2k2?V2te) = 0, (1.34)
where ?D = Vte/?pe is the electron Debye length. It can be seen from Eqs. 1.15 and
1.34 that the linear ELD only contributes to the imaginary part of the dispersion
relation
D0i =
??
k2?2D
??
2k?Vte exp(?
?2
2k2?V2te). (1.35)
Let us rewrite Eq. 1.15 as
D0r = k
2
?
k2?? +
k2?
k2??. (1.36)
Plug the D0r and D0i into Eq. 1.33. The expression of the coe?cient of the electro-
static LHWs? electron Landau damping is found to be
?/?r ? ?D0i?
r?D0r/??r
=
?
? 
k2 2D
!r?
2k?Vte exp(?
!2r
2k2?V 2te)
2[k2?k2 !2pi!2
r
+ k
2
?
k2
!2pe
!2r ]
. (1.37)
Since lots of di?erent approximations are made in the derivation of the damping
rate, there are few other expressions that were adopted by others[21]. From Eq. 1.33,
15
the following expression can be acquired
? ? ?D0i(x,k,?r)(?D
0r/?k)(?k/?r)
= ?D0i(x,k,?r)?D
0r/?k
(??r/?k), (1.38)
where, ?r takes the form of Eq. 1.16. From this respect, one can obtain the following
expression of the damping rate
? =
??
2?2
?2LH
k2c2s?
2
r
1
k?Vte exp(?
?2
2k2?V2te), (1.39)
with cs =
?
?Te/mi. Although the two expressions of the ELD rate of electrostatic
LHW (Eqs. 1.37 and 1.39) are di?erent, the numerical results from them are pretty
close to each other.
1.3.2 Physical mechanisms of the Landau damping of LHWs
Unlike most other wave modes, lower hybrid waves (LHWs) that interact reso-
nantly with electrons can also undergo resonant interactions with ions. This allows
LHWs to mediate the transfer of energy between the two species and possibly lead
to plasma heating or particle acceleration, making them of considerable interest in
many di?erent laboratory[26, 30, 31] and space[32, 33, 34] plasma environments. As
?i ? ?LH ? ?e, when considering interactions with LHWs, the electrons must
be treated as being magnetized, their gyro-motions cannot be neglected, and they
interact resonantly with the waves when the condition
? = k?ve? (1.40)
is satis?ed. However, the ions are often treated as unmagnetized, so their gyro-
motions may be neglected and they interact resonantly with the waves when the
16
condition
? = k?vi (1.41)
is satis?ed. Here ve and vi are the velocities of an electron and ion, respectively.
Since k? ?k? and ve ?vi for typical particles if the ion and electron tempera-
tures are similar, LHWs with the same k can satisfy both Eqs. 1.40 and 1.41 provided
that vi has a signi?cant component perpendicular to B. Hence, LHWs can interact
resonantly with both ions and electrons and transfer energy between the parallel mo-
tions of electrons and the perpendicular motions of ions. Thus, if the ion e?ects are
included, the linear Landau damping rate of LHWs can be obtained as
?
?r =
??
2
?2LH
k2c2s{
?r?
2k?Vte exp(?
?2r
2k2?V2te) +
Te
Ti
?r?
2kVti exp(?
?2r
2k2V2ti )}, (1.42)
where Vti =
?
?Ti/mi is the ion thermal speed.
1.4 Parametric instability
In this section, general theories of PI of LHWs are introduced[63, 21]. The dis-
persion relation equation near lower hybrid wave frequency is reviewed in detailed
algebras. By including the e?ects of ion nonlinearity, the various channels of para-
metric decay of LHWs are discussed.
1.4.1 The PI dispersion relation near lower hybrid frequency
Consider the system with a static magnetic ?eld inz direction and a lower hybrid
wave (pump wave) with frequency and wave number (?0,k0) to propagate in x-z
plane. A low frequency wave mode (?,k) and two sidebands (? ??0,k ? k0) are
considered to be the decay waves.
17
Equations of ions response
For the ions, with the presence of the electrostatic potentials of the pump and
decay waves, the distribution can be decomposed to be a unperturbed part (f00i), a
response part at the pump (f0i), responses at sidebands (f1;2i) and a response at low
frequency (fi).
Fi = f00i +f0i +f1;2i +fi, (1.43)
where f00i is Maxwellian distributed. For the high frequency response, ions can be
assumed to be unmagnetized, and thus one can obtain
f0i = e?0k
BTi
k0 ?v
?0 (1 +
k0 ?v
?0 )f
0
0i exp[?i(?0t?k0 ?x)], (1.44)
f1i = fl1i +fnl1i, (1.45)
fl1i = e?1k
BTi
k1 ?v
?1 (1 +
k1 ?v
?1 )f
0
0i exp[?i(?1t?k1 ?x)], (1.46)
fnl1i = e2m
i?1
(1 + k1 ?v?
1
)[k0 ? ?fi?v??0 ?k? ?f
?
0i
?v ?]. (1.47)
Here k1 = k?k0, ?1 = ???0, k1vi <?1, k0vi <?0. ?0 = ?0 exp[?i(?0t?k0?x)] is
the electrostatic potential of the pump wave. In the same wave, ? and ?1;2 denote the
electrostatic potentials of the decay waves. From the distributions, one can obtain
the ion density ?uctuations
n1i = nl1i +nnl1i, (1.48)
nl1i = ek
2
1?1
mi?21 n
0
0, (1.49)
nnl1i = e2m
i?1
[?k0 ?k1?
1
??0ni + k?k1?
1
?n?0i], (1.50)
where ni is the response at the low frequency and n0i at the frequency of the pump,
n0i = ek20?0n00/mi?20. (1.51)
18
The evolution ofFi is governed by the Vlasov equation in terms of guiding-center
coordinates Xg, magnetic momentum ?, polar angle ? of the perpendicular velocity,
and the parallel momentum p?.
?Fi
?t + ??
?Fi
?? +
???Fi
?? ?e
?
?z?
?Fi
?p? +
?Xg ? ?Fi
?Xg = 0, (1.52)
where Xg = (x +?sin?)ex + (y??cos?)ey + zez = (xg,yg,zg), ? = v?/?i, ?i =
eB/mic is the ion gyro-frequency, ? = miv2?/2?i, and the over-dots denote the time
derivatives. All the over-dots components can be easily derived by using the equations
of motions
?? = ? e?
i
???? ?v? = ?H?? , (1.53)
?? = ?(?i + ?
??e?
?) = ??H
??, (1.54)
?Xg = ? e
mi?i
???
?ygex +
e
mi?i
???
?xgey +
p?
miez. (1.55)
The Hamiltonian H = ?i?i + e?? + p2?/2mi and ?? = ?0 + ?1;2 + ? is the total
electrostatic potential of the system. As we can see, (?,?), (xg,yg), and (z,pz) are
canonical variables. Eq. 1.52 follows directly from the continuity equation of ion
density in these six dimensional space.
To linearize Eq. 1.52, one can obtain the low frequency response
fli = ?exp[?i(?t?k?Xg)]?
i
e?exp[?il(???)]Jl(k??)
??k?v? ?l?i (l
?f00i
?? +k?
?f00i
?p? ). (1.56)
Here, the identity exp[?i(?t?k?x)] ? exp[?i(?t?k?Xg)]?i exp[?il(???)Jl(k??)]
is used. ? is the angle between k? and k0? (similarly ?1 is the angle between k1?
and k0?), and v? = p?/mi is the parallel velocity of ion. For the nonlinear part of the
19
distribution associated with the low frequency response, one can obtain
?i(??k?v?)fnli ??i(?f
nl
i
?? ) = ?1 +?2, (1.57)
where ?1 and ?2 are from the contributions of lower sideband and upper sideband,
respectively.
?1 = 12iek0??0?cos??f1i?? + 12iek1??1?cos(???1)?f0i??
? 12iek0? sin?m
iv?
?0?f1i?? ? 12iek1? sin(???1)m
iv?
?1?f0i??
+ 12ik0?e?0?f1i?p
?
+ 12ik1?e?1?f0i?p
?
+ eik1y2m
i?i
?1(ik0?f0i)? iek0?2m
i?i
?0(ik1yf1i). (1.58)
?2 can be obtained by replacing ?0, f0i, ?0, and k0 by ??0, f?0i, ??0 and ?k0, re-
spectively. Taking the approximation k1? ? k1 and k1 ? v < ?1, Eq. 1.58 can be
simpli?ed
?1 = ie
2?0?1f0
0
2kBTimi [?
?
?20k0 ?k1 +
k0 ?k1
?20 (k?v?
2?
?0 k0 ?v)
? mik
BTi
k0 ?vk1 ?v(???2
0
+ k?v?2
0
? 2??3
0
k0 ?v)]exp[?i(?t?k?x)]. (1.59)
By solving Eq. 1.57, the following expression for the low frequency component of the
ion density can be obtain by reusing the identity.
ni = nli +nnli , (1.60)
nli = ? k
2
4?e?i?exp[?i?t?k?x]
= ? ?
2
pi?
2?ev2i [1 +
?
l
?
k?viZ(
??l?i
k?vi )Il(bi)exp(?bi)]exp[?i(?t?k?x)],(1.61)
20
and
nnli = (?1?0?1 +?2??0?2)exp[?i(?t?k?x)], (1.62)
where bi = k2??2i/2, Il(bi) is the Bessel function of the second kind, Z is the plasma
dispersion function. ?1 can be obtained to be the following expression
?1 = n
0
0e
2k0?k1?
2kBTimi?20
?
k?vi
?
i
Z(??l?ik
?vi
)((?1 + 2l?i? ? 4k0??k
??0
cos?)cos?1I1 exp(?bi)
+ i2k0?k?vi?i?
0
sin?cos?1 ddb
i
[Il exp(?bi)] + 2k0?l?ik
??0
d
dbi[biIl exp(?bi)][cos(???1)
+ 12 cos(?+?1)? 12 cos(3???1)] + (1? l?i? )cos?1 ddb
i
[Il exp(?bi)]
? ik0?k?vi?i?
0
d2
db2i [biIl exp(?bi)][sin(???1) +
1
2 sin(?+?1)?
1
2 sin(3???1)]
+ 12{?2?i? isin(2???1) + (1? l?i? )cos(2???1) + 2k0??ik
??0
[lcos(3???1) + 2isin(3???1)]}
? {4l
2Il exp(?bi)
k2??2i ?2Il exp(?bi)?2(1 +bi)
d
dbi[Il exp(?bi)]}
+ 12{?2?i? cos(2???1) + (1? l?i? )isin(2???1) + 2k0??ik
??0
[2cos(3???1) +ilsin(3???1)]}
? {?2l ddb
i
[Il exp(?bi)] + 4lIl exp(?bi)k2
??2i
}). (1.63)
By replacing ?0, k0 by ??0, ?k0, and the subscript 1 by 2, one can obtain ?2.
Response of electrons
Since the electron Larmor radius is much shorter than the wave lengths involved
in the parametric process, drift kinetic theory can be used to describe the response of
electrons. Referring to Tripathi et al.[63], the nonlinear electron density ?uctuations
at the decay wave frequencies can be written in the following expressions.
nnle = (?1?0?1 +?2??0?2)exp[?i(?t?k?x)],
nnl1e = ?3??0?exp?i[(???0)t?(k?k0)?x],
nnl2e = ?4?0?exp?i[(?+?0)t?(k+k0)?x], (1.64)
21
where
?1 = ?1E?B + ?, (1.65)
?1E?B = ? ie
2n0
0
?em2ev3e
k0?k1? sin?1
?k0? ??0k? (?Z??0Z0 ??1Z1), (1.66)
? = e
2n0
0
m2ev2e?2e[k
2
? ?k0? ?k1??1Z1 +k? ?k1??0Z0]
+ e
2k0?k1?k?n0
0
m2ev3e(?k0? ??0k?)2(?Z??0Z0 ??1Z1). (1.67)
Here ? = ?/k?ve, ?0;1 = ?0;1/k0;1?ve, ve =
?
2kBTe/me the electron thermal speed,
?e = eB/mec the electron gyro frequency, Z = Z(?) and Z0;1 = Z(?0;1). ?2 can be
obtained by replacing ?0, k0 by ??0, ?k0 and the subscript 1 by 2, and
?3 = ?e ??(?)?
e(?)
?1(?e ???e), (1.68)
?4 = ?e ??(?)?
e(?)
?2(?e ???e). (1.69)
? denotes the contributions from the parallel motions and polarization drift, and
?1E?B denotes the E ? B nonlinearity arising through the parallel ponderomotive
force ?mevE?B ??v?, ?1E?B ?
?
mi/me?.
Now, one can obtain the nonlinear dispersion relation by combining Eqs. 1.62,
1.64 and the Poisson?s equation
1 = ?1??
1
+ ?2??
2
, (1.70)
where
?1 = ?e(?)??(?)?
e(?)
?2pi
?20
?2pi
4k2c2s(1 +
?
k?veZ)
2 sin2?1U
2
c2s , (1.71)
22
for ? ?k?ve, and
?1 = ?e(?)??(?)?
e(?)
?2pi
?20
?2pi
k2
(mi/me)2
4(k? ? !!0k0?)2(
k3?
?2 ?
k30?
?20 ?
k31?
?21 )U
2, (1.72)
for ? ? k?ve. ?2 = ?1(?1 ??2), cs =
?
kBTe/mi the ion sound speed, and U =
ek0?0/me?e is the E?B drift velocity of electrons. ? and ?1;2 are the linear dielectric
functions at ? and ?1;2, respectively.
1.4.2 Decay channels
We have reviewed the dispersion relation in the previous introduction, here the
various decay channels of lower hybrid wave through parametric instability will be
summarized. The pump LHW with (?0,k0) decays into a low frequency wave mode
(?,k) and two high frequency sidebands (???0,k?k0).The sidebands can be lower
hybrid wave modes or ion Bernstein wave mode (? ? n?i). The various decay
channels are classi?ed into three regions according to the magnitude of k?ve/?.
Case A: k?ve/? ? 0.3
In this case, the possible channels can be
Channel 1: pump wave decays into two lower hybrid waves. ?0 > 2?LH.
Channel 2: pump wave decays into ion Bernstein and lower hybrid waves. ? ? ?i
and k??i ? 1, or ? ? 3?i with arbitrary k??i.
Channel 3: pump LHW decays into ion Bernstein reactive quasi-mode and lower
hybrid wave mode (four-wave decay).
Channel 4: the decay follows modulational instability.
The growth rates for all the channels in case A will be introduced in the following
discussions.
23
Channel 1. Decay into two lower hybrid waves. In this channel, the requirement
?0 > 2?LH needs to be satis?ed since both of the decay waves are lower hybrid waves,
which have frequencies larger than ?LH. The growth rate, dominated by the E?B
electron nonlinearity, is given by
?20 = [?
2
LH
4?0
?1=2(?0 ??)1=2mi/me
(k? ??k0?/?0)k (
k3?
?2 ?
k30?
?20 ?
k31?
?21 )]
2U2 sin2?1. (1.73)
Note that ?0 is the growth rate without taking the damping e?ects into account.
Channel 2. Ion Bernstein and lower hybrid waves. There are two cases needed
to be considered in this channel.
(1) k??i < 1. For ? ? ?i, the linear dielectric function simpli?es to
?2 = ?2i k
2
?mi
k2me(1 +
k2?mi
k2me)
?1, (1.74)
and thus one can obtain the growth rate
?0 ? (?0?i)1=2(1? ?
2
LH
?20 )
1=4( Te
2Ti)
1=2k?i?LH
?0
U
4cs sin?1. (1.75)
(2) k??i ? 1 (? ?n?i, n? 3). In this case, the dielectric function is
?? 1 + ?
2
pe
?2e ?
?2pi
?2
k2?mi
k2me +
2?2pi
k2v2i (1?
?In(bi)exp(?bi)
??n?i ). (1.76)
The growth rate turns out to be in the following expression
?0 = ?LH(k??i)
1=2
2?1=4?0 (
Te
2Ti)
1=2U
cs[
(2???0)?3i
?2(?0 ??)2 (
?21
?2LH ?1)sin?1]
1=2
? [k2??2i?20(2???0)?i?2(?
0 ??)2
( ?
2
1
?2LH ?1)sin?1 +icos?1]
1=2. (1.77)
24
Note that icos?1 term is due to the ion nonlinearity. This growth rate is comparable
to that for the decay into ion cyclotron waves. The condition for Eq. 1.77 to be
valid is ?0 <??n?i. The mode will become a reactive quasi-mode if ?0 >??n?i,
however, as the growth rate and frequency mismatch increases the ion contribution
is drastically reduced.
Channel 3. Ion Bernstein Reactive quasi-mode and lower hybrid waves (four-wave
decay). It becomes a reactive quasi-mode if the low frequency mode has ? ? k?ve
and k??i ? 1 for higher values of U/cs, or if ? ? ?i can no longer satisfy ? ? ?,
and the two sidebands are equally important. In this case, one have
? = ?r +i?, (1.78)
?r = ?2 + ?2/2, (1.79)
?2 = ??2/4 + 1/2(A?)1=2. (1.80)
Here ? = ?0 ??LH[1 + (k2?/k2)(mi/me)]1=2, k>k0, and
A = ?
2
LH
4?20 (?
2
0 ??
2
LH)
k2c2s
?0
U2
c2s . (1.81)
Note that this growth rate is less than that for the oscillating two-stream instability.
Channel 4. Modulational instability. In this channel, both high frequency side-
bands become equally important even for ? ? ?, whenever ?/k = ??0/?k0 can be
satis?ed, and the pump wave?s amplitude gets modulated.
??0
?k0? =
?0
k0?(1?
?2LH
?20 ) ?
?
k
k?
k, (1.82)
??0
?k0? = ?
?0
k20?(1?
?2LH
?20 ) ?
?
k
k?
k . (1.83)
25
Since k0? ?k0?, these conditions are satis?ed for k? ?k?, that is the low frequency
wave mode propagates mainly along the magnetic ?eld direction. The growth rate
can be obtained as
? = k?ve?
0
?(mem
i
)3=2(1? ?
2
LH
?20 )(4?
3?2LH
?20 )sin?1
U
cs. (1.84)
It can been seen that the growth rate is extremely low. This channel evolves into case
A, channel 3 for a larger value of U/cs.
Case B. 0.5 <k?ve/?< 2.0
This case corresponds to the nonlinear Landau damping by electrons and has the
following decay channels.
Channel 1: lower hybrid quasi-mode and lower hybrid wave.
Channel 2: ion cyclotron or ion Bernstein quasi-modes and lower hybrid wave.
Ion cyclotron mode has parallel phase velocity smaller thanve, while the ion Bernstein
quasi-mode?s parallel phase velocity is larger than ve.
The growth rates for all the channels in this case are estimated as well.
Channel 1. Lower hybrid quasi-mode and lower hybrid wave. In this channel
the pump wave and the lower sideband can interact resonantly with the electrons
under the resonance condition ? = ?k?ve. This condition can be satis?ed only when
?i ?? ??0 for?0;1 ?k0;1?ve. Thus, the growth rate in this Chan can be calculated
as
? ???1l + 0.6(?
2
LH
8?0 )sin
2?1(U
2
c2s ). (1.85)
The growth rate is very high and is believed to dominate the region with high density.
Channel 2. Ion cyclotron quasi-mode and lower hybrid wave. This channel may
occur when the low frequency quasi-mode ? ? ?i or ? ? n?i. In this channel, the
26
growth rate for the resonance condition ? ?k?ve takes the form
?0 = i?02 ?
4
pi
?20
sin2?
4k2c2s
U2
c2s [1 +Z(?1)]
2 (?e ??)/?e
(1 +?2pe/?2e)?. (1.86)
To simplify this expression, one can obtain
?0 = ?
2
LH
?20 sin
2?10.6
8
U2
c2s . (1.87)
The quasi-modes decay channel can occur with another possibility when???i ?
k?vi, the nonlinear ion cyclotron damping condition. In this channel, the growth rate
can be estimated by taking the limit ? <???i ?k?vi ??<k?ve.
?0 = 0.68 ?
2
LH
?0 sin
2?1U
2
c2s
Te
Ti
?I1 exp(?bi)
k?vi ?, (1.88)
where
??1 = 1 + TeT
i
[1?I0 exp(?bi) + 0.3?I1 exp(?bi)k
?vi
]2 + [0.6TeT
i
?I1 exp(?bi)
k?vi ]
2. (1.89)
Case C. k?ve/?> 0.3
In this case, there exists mainly three channels.
Channel 1: ion acoustic and lower hybrid waves (only for ?0 > 4?LH and Te >
4Ti)
Channel 2: ion cyclotron and lower hybrid waves for k??i > 1
Channel 3: oscillating two-stream instability
For all the channels, the growth rates can be obtained.
Channel 1. Ion acoustic wave and lower hybrid wave. The condition for this
channel to occur is that k?ve ?? ?kcs ?kvi and, hence, ?0 > 4?LH and Te > 4Ti.
27
For this channel, one can obtain the growth rate to be
?20 = [?LH4?
0
(??0)1=2 sin?1Uc
s
]2. (1.90)
In this case, the low frequency mode is a resonant quasi-mode.
Channel 2. Ion cyclotron and lower hybrid waves. In this case, the growth rate
is
? = ?LH?
0
[I1 exp(?bi)?i?0]1=2
2?2(1 +Ti/Te) sin?1
U
cs. (1.91)
The approximations ? ? ?i, ?? ?i ? k?vi, k? ? ? (k??i > 1), ? = ? +i? and
? ????i are adopted in deriving the growth rate.
Channel 3. Oscillating two-stream instability. There are two possibilities in this
decay channel.
(a) ? ?k?ve, k?vi. In this case, the dielectric function is
?? (?2pi/k2c2s)(1 +Te/Ti). (1.92)
To consider the both high frequency sidebands, one can obtain the nonlinear disper-
sion relation
?2 = ?2 + ?
2
LH
?20
?0?
4(1 +Ti/Te)
U2
c2s , (1.93)
where k0? ?k?, ? ??0 and
? = ?0 ??LH(1 + k
2
?
k2
mi
me)
1=2. (1.94)
For the purely growing mode to occur, ? < 0 and
?2LH
4?0(1 +Ti/Te)
U2
c2s >|?|. (1.95)
28
One can obtain the maximum growth rate to be
?max = ?
2
LH/?0
8(1 +Ti/Te)
U2
c2s . (1.96)
(b) k?vi ? ? ? k?ve, compatible with ? ? k?ve, and Im? ? ?i, the dielectric
function is
?? ?
2
pi
k2c2s(1 +
Te
Ti +i
Te
Ti
??/?i
(2?bi)1=2). (1.97)
If
Te
Ti
?
?i
?
(2bi?)1=2 ? 1, (1.98)
then the dispersion relation becomes the same as Eq. 1.93 and the growth rate is
given by Eq. 1.96.
1.5 Summary
In this chapter, the physics of Landau damping and parametric instability of
lower hybrid waves are reviewed. The theory of Landau damping, including both the
linear and nonlinear Landau damping rates of Langmuir waves, is introduced. The
main characteristics of the nonlinear Landau damping, such as the particles trapping
and long time evolution into a ?nal BGK equilibrium, etc, are discussed by referring
to the recent published literatures. For the lower hybrid wave, the electromagnetic
dispersion relation as well as the electrostatic dispersion relation in both cold and
warm plasmas is introduced in great details. The coe?cient of the linear electron
Landau damping of LHWs is derived given the condition that the damping rate is
much smaller than the real frequency. The important property of the Landau damping
of LHWs, that unlike most other wave modes the LHWs can interact with both
electrons and ions directly, is ?nally discussed. In the last section, the theory of PI is
reviewed including both the electron and ion nonlinear response e?ects. The nonlinear
29
dispersion relation equation of the PI of LHWs is given as Eq. 1.70. Based on the
theory and dispersion equation, various decay channels are categorized according to
the quantity k?ve/?. For the various possible decay channels, the growth rates are
also estimated.
30
Chapter 2
GeFi Simulation Scheme
Numerical simulation is an advanced technique to understand the kinetic physics
of various fundamental plasma processes, especially in the investigation of the non-
linear plasma dynamics under realistic conditions, for which the analytical theory is
unable to be acquired. Full-particle simulations have been utilized for decades[37,
38, 39, 40]. In the full-particle codes, both electrons and ions are treated as fully
kinetic particles, and thus the kinetic physics of both electrons and ions can be in-
cluded in it. However, an arti?cial, small ion-to-electron mass ratio is often assumed
in the full-particle codes in order to accomodate the computation resources. Another
kinetic simulation approach is the hybrid simulation[41, 42], in which the ions are
treated as fully kinetic particles, while the electrons are treated as a massless resis-
tive ?uid. Thus, the electron kinetic e?ects are absent in the hybrid model. Neither
full-particle code scheme nor hybrid code scheme is appropriate to solve the physics
of LHWs, which must include the dynamics of both electrons and ions, and requires a
realistic ion-to-electron mass ratio. Note that the frequency of LHWs is in the range
? ??LH ? ?i, well above the range of the hybrid models.
A novel new simulation code model, GeFi model, has been developed by Lin, et.
al.[35], recently. In the GeFi code model, the electron dynamics is determined by the
gyrokinetic theory[43, 44, 45, 46, 47], and the ions are treated as fully kinetic particles
governed by the Vlasov equation. In the GeFi model, the rapid electron gyro-motion
is removed while ?nite Larmor radius e?ects are retained. This treatment allows us to
deal with realistic ion-to-electron mass ratio and ?nite Larmor radii. This new model
requires that the electron gyrokinetic approximation is valid, and thus is particularly
31
suitable for the dynamics with wave frequency ? < ?e and k? < k?. It can be seen
from the discussions in chapter 1 that the lower hybrid wave falls exactly in this
range. It is therefore suitable to simulate the dynamic physics of the LHWs, such as
its Landau damping, parametric instability, and the electron-ion hybrid instability,
using the GeFi code model.
2.1 GeFi Scheme Algebras
The GeFi simulation scheme treats ions as fully kinetic (FK) particles and elec-
trons as gyrokinetic (GK) particles.[35] For the FK ions, the dynamics is governed
by Vlasov equation in the six-dimensional phase space (x, v).
?fi
?t +v?
?fi
?x +
qi
mi(E+
1
cv?B)?
?fi
?v = 0. (2.1)
mi is the ion mass, qi is the ion charge, E is the electric ?eld, B is the magnetic ?eld,
and fi is ion distribution function. fi is represented by a group of particles
fi(x,v,t) = ?
j
?[x?xj(t)]?[v?vj(t)], (2.2)
where the index j represents individual particles. The evolution of fi is determined
by ion motion under self-consistent electromagnetic ?elds, i.e.
dv
dt =
qi
mi(E+v?B),
dx
dt = v. (2.3)
The number density, ni and current density, Ji can be obtained
32
ni =
?
fid3v = ?
j
?(x?xj),
Ji = qi
?
fivd3v = qi?
j
vj?(x?xj). (2.4)
The electrons are treated as gyrokinetic particles. The following GK ordering for
electrons is adopted
?
?e ?
?e
L ?k??e ?
?B
B ??,
k??e ? 1. (2.5)
Here, ?e = vte/?e denotes the electron Larmor radius, L is the macroscopic back-
ground plasma scale length, ?B represents the perturbed magnetic ?eld, and ? is a
smallness parameter. The coordinates (x,v) are transformed to the gyrocenter coordi-
nates (R,p?,?,?), where R = (x??) is the gyrocentre position with? = (b?v?e)/?e
being the gyroradius vector, p? = meve? +qeA?/c the parallel canonical momentum
of electrons, ve? and ve? the parallel and perpendicular velocities of electrons, re-
spectively, qe the electron charge, ? the magnetic moment, b = B/B, B = ?B +?B
with ?B being the background magnetic ?eld averaged over the spatial and temporal
scales of wave perturbations, and ?B = ??A the perturbed magnetic ?eld. The
parallel direction is de?ned along the background magnetic ?eld ?B. The following
GK equation can be obtained by averaging the Vlasov equation over the gyrophase
angle ?
?Fe
?t +
dR
dt ?
?Fe
?R +
dp?
dt
?Fe
?p? = 0, (2.6)
33
where Fe(R,p?,?) is the distribution function of electron in the ?ve-dimensional gy-
rocenter phase space. The gyrocenter equations of motion for p? and R are
dp?
dt = ?b
? ?[qe <??? > +???B],
R
dt = ve?b
? + c
qe ?B
?b?[qe <??? > +???B], (2.7)
where b? = ?b + (ve?/?e)?b ? (?b ? ?)?b, ?b = ?B/?B, ?? = ? ? v ? A/c, ? and A
are perturbed scalar and vector potentials, respectively, and the operator < ??? >
represents gyro-averaging. The electron gyro-averaged guiding centre charge density
and p? current are
<Ne > =
?
Fed3v,
<Je? > = qem
e
?
p?Fed3v. (2.8)
In GK simulations, the gyro-averaging is carried out numerically on a discretized
gyro-orbit in real space[48].
In order to advance Fe and fi, we need to calculate the perturbed potentials and
?elds from the Maxwell equations. The Poisson?s equation becomes
?2?? = ?4?(qini +qene), (2.9)
where ne is the electron number density, and for w = v2/2,
ne = qem
e
?
d3v(?
?fe
?w)[??<?> +
1
c <v? ?A>]+ <Ne >. (2.10)
Assuming |?2?| ? |?2?|, we have replaced ?2? by ?2?? in Eq. 2.9 to suppress the
undesirable high-frequency Langmuir oscillation along B. Eq. 2.9 along with 2.10
34
then become the following generalized GK Poisson?s equation
(1 + ??
2
pe
??2e )?2??+ 4??neqe
?B?
?B = ?4?(qini +qe <Ne >), (2.11)
where ?ne is the spatially averaged electron density, ?B? = ?b ??B, and the second
and third terms on the left-hand side correspond to the electron density due to its
perpendicular guiding-centre polarization drift of the electrostatic electric ?eld and
E?B drift associated with inductive electric ?eld ?A?/?t, respectively.
Since ? ? ?e, the electron force balance instead of the usual perpendicular
Ampere?s law is used to calculate ?B?.
??(neqeE) = ??[??Pe ? 1cJe ?B], (2.12)
where
Pe = (?neqe?2e?2??+ 2?neTe?
?B?
?B )(I?
1
2
?B?B)+ <Pg >,
<Pg > =
?
mevvFed3v, (2.13)
and, analogous to the derivation of ne, the ?rst two terms in the expression for Pe are
due to the electron perpendicular guiding-centre drifts. Assuming zero background
electric ?eld, E = ?E. Noting that ?E = ??? ? (1/c)?A/?t, ? ? A = 0, the
electron current density Je = (c/4?)??B?ji with ji being the ion current density,
|?2?| ? |?2?|, and ignoring corresponding higher-order terms, Eq. 2.12 can then be
shown as
?2? = ???(??Pg + 1cJi ?B), (2.14)
where, noting ?neqe = ??niqi,
? = (1 +
??e) ?B?B?
4? ? ?niqi(1 +?
2
e?
2
?)?. (2.15)
35
Expressing ?B? in terms of ?, given by Eq. 2.15, the GK Poisson?s equation, Eq.
2.11, can ?nally be expressed as
[(1 + ??e + ??
2
pe
??2e )?2? ?
??2pi
?V2A ]? = ?4?[(1 + ??e)(qini +qe <Ne >)?
4??niqi
?B2 ?], (2.16)
where ??pi and ?VA are the background ion plasma frequency and the Alfv?en speed,
respectively.
The perturbed potential A?, meanwhile, is given by the following parallel Am-
pere?s law
(?2 ? ?
2
pe
c2 )A? = ?
4?
c (Ji?+ <Je? >). (2.17)
Given A? and ?B?, the vector potential A can be calculated. Let us decompose
A into three locally orthogonal components, i.e. A = A? +A??b+???. A? is then
determined by the perpendicular Ampere?s law
?2A? = ?4?c J?, (2.18)
where J? = (c/4?)???B?. ???, meanwhile, can be determined by the Coulomb
gauge ??A = 0 or ?2?? = ???(A??b).
With A being completely speci?ed, the perturbed magnetic ?eld ?B is simply
?B = ?A. Since ? ? ?e, the electric ?eld E that goes into the ion equation motion
can be calculated from the following electron force balance equation,
neqe?E = ???Pe ? 1cJe ?B. (2.19)
The calculated ?E and B = ?B+?B can then be used to advance ions.
36
2.2 GeFi model in the electrostatic limit
In this thesis, all work is done with electrostatic GeFi particle simulations, al-
though the GeFi model is electromagnetic. When considering the GeFi model in the
electrostatic limit, equations can be simpli?ed in the following.
For ions, the dynamics is governed by the Vlasov equation in six-dimensional
phase space. (Here, all symbols have the same meanings as that in the last section.)
?fi
?t +v?
?fi
?x +
qi
mi(E+
1
cv?
?B)? ?fi
?v = 0. (2.20)
The evolution of fi is determined by ion motion under self-consistent electromagnetic
?elds, i.e.
dv
dt =
qi
mi(E+v?
?B),
dx
dt = v. (2.21)
In the gyrokinetic approximations for electrostatic electrons, the gyrocenter equa-
tions of motion for parallel electron momentum p? = meve? and the electron gyrocen-
ter position R are given by
dp?
dt = ?b
? ?[qe <??> +???B],
R
dt = ve?b
? + c
qe ?B
?b?[qe <??> +???B], (2.22)
where b? = ?b+(ve?/?e)?b?(?b??)?b, ?b = ?B/?B, ? is electrostatic potential, and the
operator <???> represents gyro-averaging.
Assuming the gyrokinetic condition, k?/k ? 1 and thus ?2? ? ?2?, the electric
?eld is solved by the generalized gyrokinetic Poisson equation for the perturbed scalar
potential in electrostatic limit
37
(1 + ??
2
pe
??2e )?2?? = ?4?(qini +qe <Ne >), (2.23)
2.3 Benchmark of the GeFi Scheme
The above GeFi scheme has been benchmarked for a one-dimensional uniform
system. The background magnetic ?eld is in the xz plane and the wave vector k is
assumed to be along x direction.
The top plot of Fig. 2.1 shows a comparison between the dispersion relations
of ?Bz for the fast magnetosonic/whistler branch obtained from the kinetic GeFi
simulation (open dots), and the corresponding analytical linear dispersion relations
(solid lines). The parameters in this benchmark are ?e = ?i = 0.04, mi/me = 1836,
and k?/k? = 0.2, 0.06 and 0 are plotted. In the case with k? = 0 the electromag-
netic mode approaches the quasi-electrostatic lower hybrid mode, and the frequency
?/?i =
?
mi/me = 42.8. The bottom plot of Fig. 2.1 is the comparison between the
dispersion relations of ?By for the shear Alfv?en/kinetic Alfv?en mode branch obtain
from GeFi simulation and the analytical theory for k?/k? = 0.06. The analytical
solution of the MHD shear Alfv?en mode, ?/?i = k?VA/?i = 0.42, is also shown as
the dashed line. It can be seen from the ?gure, the GeFi simulation results are in
excellent agreement with the theoretical analysis for k? ?k?.
2.4 Summary
In this chapter, the original new GeFi kinetic simulation scheme is introduced.
This new scheme, in which the electrons are treated as GK particles and ions are
treatedasFKparticles, isparticularlyapplicabletoproblemsinwhichthewavemodes
ranging from magnetosonic and Alfv?en waves to lower-hybrid/whistler waves need to
be handled on an equal footing. To utilize this code, the simulated physical processes
should be dominated by wave frequencies ?< ?e, and wave numbers k? <k?. With
38
Figure 2.1: Comparison between the dispersion relations obtained from the kinetic
GeFi simulation and the corresponding analytical linear dispersion relations for var-
ious k?/k?. Top plot is ?Bz in the fast magnetosonic/whistler branch. Bottom plot
is ?By in the shear Alfv?en/kinetic Alfv?en mode branch. The dashed line shows the
analytical dispersion relation of the MHD shear Alfv?en mode.
39
fast electron gyro-motion and Langmuir oscillations removed from the dynamics, the
GeFi model can readily employ realistic mi/me mass ratio. As discussed above, the
GeFi model has been improved and modi?ed to allow the existence of modes with
wavelengths on the same scale of the background nonuniformity. The novel GeFi
scheme has already been benchmarked for uniform plasmas to resolve the physics of
wave modes ranging from Alfv?en waves to lower hybrid/whistler waves, and for the
tearing mode instability in the magnetic reconnection in a Harris current sheet.
40
Chapter 3
GeFi Simulation of Landau Damping of Lower Hybrid Waves
3.1 Introduction
In this chapter, the GeFi model is used to simulate the Landau damping of
LHWs. Both linear and nonlinear Landau damping of LHWs in the electrostatic
limit will be studied. The results, based on an initial value problem, will provide a
basic understanding of the Landau damping of LHWs from the lower-hybrid to the
ion cyclotron time scales. The real frequency and linear damping rates from the linear
Landau damping simulations will be compared with the theoretical predications from
the dispersion relation. The nonlinear Landau damping characteristics, including the
long time evolution, nonlinear Landau damping rates and the motions of both reso-
nant electrons and ions will be studied in great details. In either linear or nonlinear
Landau damping of LHWs, electric current can be driven through the process that
waves and particles exchange energy and momentum. Although our investigations of
Landau damping of LHWs are based on an initial value problem, it is of interest for
us to calculate the current driven in the Landau damping. For various parameters,
the driven currents will be presented in the last section of this chapter.
3.2 Simulation Model
Here, it is time to introduce the normalizations and some main parameters in
the simulation of LHWs. In the calculation, the time is normalized to the inverse of
the electron cyclotron frequency ?e, lengths are normalized to the electron Larmor
radius ?e, velocities are normalized to the electron thermal speed Vte, electric ?eld
41
is normalized to En = VteB0/c (with B0 being the background magnetic ?eld and c
the light speed), temperatures are normalized to the electron temperature Te, and
densities are normalized to the equilibrium electron particle density n0. The particle
initial distribution is assumed to be Maxwellian. The simulation is carried out as an
initial value problem in a one-dimensional uniform domain in the x direction (i.e., the
wave vector k = kxcex). Periodic boundary conditions are used for both boundaries.
The main parameters are set as the following. The ion-to-electron mass ratio is
mi/me = 1836. The grid number Nx varies from 256 to 2048 for di?erent cases, and
the particle number per cell varies from 1000 to 10,000 for both electrons and ions.
The time increment is ?e?t = 0.5. The ratio c/VA = 42.8, and VA is the Alfv?en
speed. The electron plasma frequency is thus ?pe = ?e(c/VA)
?
me/mi = 1.0, and
the ion plasma frequency ?pi = ?pe/?1836 = 0.023. The ion-to-electron temperature
ratio Ti/Te varies from 0.33 to 4.0. The background magnetic ?eld, which has a ?xed
magnitude |B0| = 1.0, is assumed to be in the xy plane, with the ratio B0x/B0y
varying from 0.001 to 0.1. Accordingly, the ratio k?/k varies from 0.001 to 0.1.
3.3 Simulation of linear ELD of LHWs
The linear electrostatic LHWs in a uniform plasma have been simulated using
the GeFi scheme, under a small perturbed electric amplitude E0 = 0.001 (E0 is the
initial amplitude). As discussed in Chapter 1.1, ?L? ? 1, we consider the process
mainly as nonlinear Landau damping, otherwise, we consider the damping as linear
Landau damping. Here, ? =
?
me/eE?k? is the electron trapping frequency in LHWs.
Thus, in order to satisfy the linear Landau damping condition, we need to set up a
small electric amplitude. In our case, E0 = 0.001 is small enough to go through a
linear ELD. Moreover, since the electric ?eld is small, it is necessary for us to lower
down the noise level to observe obvious decay. About 10,000 particles are used in
each grid in order to simulate the LHWs in the small amplitude, at low noise levels.
42
In each run, only a mode with a single k is kept in the domain, while all other modes
are ?ltered in order to single out the damping source in the electron wave-particle
interaction.
Figure 3.1: Perturbed electrostatic potential ?1 (in logarithm scale) vs. ?it for a case
of linear electron Landa damping of LHWs
Fig. 3.1 is a plot of perturbed electrostatic potential ?1, in a logarithm scale,
as a function of time ?it for a typical case of linear ELD of LHWs. In this case
k?/k = 0.066, an electron Landau damping case. The red line shows approximately
the noise level. In order to obtain a clear linear damping in logarithm, the grid
number Nx = 2048 and the particle number in cell 10,000 are used. It is seen from
Fig. 3.1, the amplitude of the potential decays into the noise level directly, and never
bounce back. Since the time for the amplitude to decay into noise level is very short,
the nonlinear e?ects can not be involved, thus the purely linear ELD happens.
Here, it is necessary to clarify that our simulations are performed only 1D in x
direction, although the LHWs propagate nearly perpendicular to the external static
43
magnetic ?eld. In order to simulate the LHWs, we put the magnetic ?eld at x and y
directions, and put the wave vector k in only x direction. Thus the important ratio,
k?/k in the real frequency calculation, is just B0x/B0. The LHWs are embedded into
our simulation system by perturbing the ion density, ni0 = a0 cos(kxx) with a0 being
the initial amplitude of the perturbed ion density. Note that the perturbation is not
depending on time, thus this is an initial value problem.
In order to compare our simulation results of the linear ELD of LHWs with the
corresponding linear theories, the real frequency and the linear ELD rate as a function
of the wave number, k?e are plotted as Fig. 3.2 for magnetic ?eld B0x = 0.0659 and
B0y = 0.999, which corresponds to a ?xed k?/k = 0.066, while k?e varies from 0.15 to
0.30. The temperature ratio Ti/Te = 1.0, and the initial wave amplitude E0 = 0.001.
In Fig. 3.2, the top plot is the real frequency and the bottom one is the linear ELD
rate, and the error bars are estimated based on the uncertainties in FFT introduced
to the frequency space due to discrete data points in the ?nite time sequence. It is
seen from Fig. 3.2 that the real frequency has a good match with the theoretical
predictions from Eq. 1.21, which is shown by the solid curve in the ?gure. The
linear ELD rate, ?L/?i as a function of k?e also agrees very well with the theoretical
predications given as Eq. 1.37. For k?e < 0.19, we expect very low linear electron
Landau damping rate, i.e., ?L/?i ? 0. The damping rate grows nearly exponentially
with respect to k?e. The results of Fig. 3.2 show that the GeFi model is really good
for the simulation of LHWs.
3.4 Nonlinear Landau Damping of Lower Hybrid Waves
As discussed above, we can increase the perturbed electric ?eld, E0, to obtain the
nonlinear ELD of LHWs. The nonlinear ELD of LHWs is discussed in this section
with lots of interesting properties. In Section 1.3.2, we have discussed that both
electrons and ions can interact with the LHWs. The electrons are usually treated
44
Figure 3.2: Diamonds show the real frequency and the linear ELD rate as a function
of the wave number for LHWs withB0x = 0.0659 andB0y = 0.999, which corresponds
to a ?xed k?/k = 0.066, and Ti/Te = 1.0. The solid lines are based on the analytical
theory.
as magnetized, their gyro-motions cannot be neglected, and they interact resonantly
with the waves when the condition, ? = k?ve?, is satis?ed. However, the ions are
always treated as unmagnetized, so their gyromotions may be neglected and they
interact resonantly with the waves when the condition, ? = k?vi, is satis?ed. Here,
ve and vi are the velocities of an electron and ion, respectively. Here, in this section,
both electron and ion Landau damping of LHWs will be studied. The long time
evolution of the propagation of LHWs and the particles trapping in the nonlinear
Landau damping will be presented in details. Among other runs, three cases are
discussed in details. In case 1, only electrons are resonant with LHWs. Only ions are
resonant with the wave in case 2. Both electrons and ions are resonant in case 3.
45
3.4.1 Case 1, Electron Landau Damping of LHWs
Wave pro?le in case 1
In case 1, k?/k = 0.066, the ion-to-electron temperature ratio Ti/Te = 1.0, and
the initial perturbed wave amplitude E0 = 0.1, the wave number k?e = 0.2255, the
real frequency is thus calculated, ? = 103.41?i = 3.03?LH. The resonant condition
for electrons are calculated, ?/(k?Vte) = 3.78, thus the electrons at the tail part of
the velocity distribution can satisfy the electron resonant condition, ? = k?ve?. On
the other hand, the ion resonant condition is calculated, ?/(kVti) = 10.69 ? 1 (with
Vti =
?
?Ti/mi being the ion thermal speed), which results in that almost none of
the ions can satisfy the resonant condition, ? = k?vi.
Figure 3.3: Time evolution of the spatial Fourier mode of the electric ?eld (in the
natural logarithm scale) for case 1 with k?/k = 0.066 and Ti/Te = 1.0, and k?e =
0.2255, E0 = 0.1.
Fig. 3.3 shows the time evolution of the spatial Fourier modes of the electric ?eld
(in the natural logarithm scale) in case 1. Compared with Fig. 3.1, which is purely
46
linear ELD, the LHW pro?le of case 1 is seen to go through an electron nonlinear
Landau damping, in which the wave amplitude oscillates and slowly decreases. After
a few oscillations, the amplitude settles down to a steady oscillation, which is very
above the noise level. In essence, As discussed in section 1.1.2, O?Neil?s[2] theory has
already indicated that large amplitude waves trap the resonant electron population,
and thus the Landau damping is dominated by the nonlinear e?ects. The nonlinear
damping rate can be measured as ?/?i = 0.07 based on the red dash line in the
?gure. As the electrons bounce back and forth inside the potential well of the wave,
the wave amplitude oscillates with nearly the particle trapping time scale. The phase
mixing of the trapped electrons causes the wave amplitude to reach an asymptotic
constant, forming a stable equilibrium with an undamped nonlinear plasma wave,
i.e., the BGK equilibrium.[4]
Particle motion, electron trapping in case 1
The particles? trapping in the wave electric ?eld potential plays an important
role in the nonlinear ELD of LHWs. In order to see the evolution of the electrons
trapped in the electron Landau damping of LHWs, it is shown in Fig. 3.4 the contour
plots of the electron distributions in the (x,ve?) phase space and the corresponding
average parallel velocity distribution functions at times ?it = 0,0.22,2.41 and 8.78
obtained in case 1, in which only electrons are resonant with the waves. The electron
distribution functions are plotted with the black solid curves, and the ion distribution
functions are shown with the red curves. Note that for the purpose of presentation,
here, the ion velocities are normalized to the ion thermal speed Vti. It can be seen
from Fig. 3.4a that the electron distribution gradually forms a vortex structure in the
phase space at ?it = 0.22, roughly around the parallel (to the magnetic ?eld) phase
velocity ver/Vte = ?/(k?Vte) = 3.785 of the LHW. Note that a similar structure can
be found at the corresponding negative velocities since the wave propagates in both
47
the positive and negative directions, although only the ve? > 0 space is shown in the
contour plots. Correspondingly, in Fig. 3.4b for ?it = 0.22, the black solid line at
ve? > 0 side shows that the electron parallel velocity distribution function deviates
strongly from the initial Maxwellian distribution in the region around the resonance
speed ver/Vte = 3.785, as marked by the black vertical dashed line. On the other
hand, nothing happens to the ion velocity distribution function, as shown by the red
solid line in the same plot, con?rming that no ions resonate with the LHW.
As time increases, more electrons are trapped, and the vortex structure as well
as the pair of ?plateaus? in the distribution function ?nally reach a steady state,
corresponding to the ?nal BGK equilibrium. The steady structures are propagating
with the wave along the x direction. Therefore, the number of electrons that can be
trapped in the waves is ?nite for any arbitrarily long time.
3.4.2 Case 2, Ion Landau Damping of LHWs
Wave pro?le in case 2
In case 2, k?/k = 0.001, the ion-to-electron temperature ratio Ti/Te = 1.0, the
initial perturbed wave amplitude E0 = 0.1, and the wave number k?e = 0.2255. The
real frequency is then calculated, ? = 37.24?i = 1.09?LH. The resonance condition
for electrons yields to be ?/(k?Vte) = 90 ? 1, and for ions, the condition becomes
?/(kVti) = 3.85. Thus none of the electrons can satisfy the resonant condition,
? = k?ve?, but ion at the tail part of the distribution can satisfy the resonant condition
? = k?vi. Therefore, only ions can be resonant with the LHWs in this case.
Fig. 3.5 shows the time evolution of the spatial Fourier modes of the electric
?eld, in a natural logarithm scale, in case 2, in which only the ion Landau damping
is present. It is seen from the ?gure, the wave amplitude decreases linearly with time
on the logarithm to a level that is much smaller than the saturation level for case 1
that is dominated by the electron nonlinear Landau damping, as seen in Fig.3.3. The
48
Figure 3.4: (a) Contour plots of the electron distributions in the particle phase space
(x,ve?) and (b) the corresponding parallel velocity distribution functions averaged
overx, shown with the black solid curves, at times ?it = 0,0.22,2.41 and 8.78 obtained
from Case 1. All the distribution functions are plotted in the logarithm scales The red
solid lines show the ion distribution functions, and the black dashed lines show the
electron resonant phase velocity ver = 3.785Vte based on the theoretical prediction.
The ion velocities are normalized to the ion thermal speed.
49
Figure 3.5: Time evolution of the spatial Fourier mode of the electric ?eld (in the
natural logarithm scale) for case 2 with k?/k = 0.001, Ti/Te = 1.0, k?e = 0.2255, and
E0 = 0.1.
simulation thus indicates that ions interact with LHWs mainly through the linear
Landau damping. After the wave has reached a weak (but still above the noise) level,
the wave amplitude is seen to oscillate on a long time scale with period ? 2?/?i.
Such long-period pattern is due to that on the time scale much longer than the LHW
period, ions are still magnetized, trapped, and gyrating with frequency ? ? ?i.
Therefore, the Landau damping dynamics is of a broad/hybrid frequency range of
both ?LH and ?i. In addition, the wave amplitude shows a sudden jump, followed
by a sudden damp, at time intervals of ?t ? 2?/?i. This sudden jump and damp
phenomenon is due to the phase bunching of ions, and will be investigated in section
3.6.
50
Particle motion in case 2
Here, we present the particle dynamics in case 2, in which only ions interact
directly with the LHW through the Landau damping. Fig. 3.6 shows the resulting
contour plots of the electron phase-space distribution in case 2. The ion phase space
evolution is shown in Fig. 3.6b, and the corresponding electron (black solid lines)
and ion (red solid lines) distribution functions are depicted in Fig. 3.6c. Again, the
ion velocities are normalized to the ion thermal speed.
Since electrons do not resonantly interact with the LHW, their velocity distribu-
tion is nearly unchanged. The phase-space structure of ion distribution, however, is
signi?cantly di?erent from that of the electrons in the ELD. The vortex structure does
not form in the ion phase space, as seen for ?it = 0.11 and 1.21, in Fig. 3.6b. The
ions are weakly trapped (i.e., trapped in the ion gyro-period time scale, much longer
than that of the LHWs), and the nonlinear e?ects are absent in their interactions
with the LHW. As a result, the ions are continuously heated, as seen at ?it = 4.40.
In Fig. 3.6c, the tail of the ion velocity distribution function expands due to the
increase of the ion thermal speed.
3.4.3 case 3, Electron And Ion Landau Damping of LHWs
Wave pro?le in case 3
In case 3, k?/k = 0.0404, ion-to-electron temperature ratio Ti/Te = 4.0, the
initial perturbed electric ?eld E0 = 0.1, and the wave number k?e = 0.2255. The real
frequency in this case? = 68.24?i = 2.05?LH. Thus, the electron resonance condition
?/(k?Vte) = 4.08, and the ion resonance condition ?/(kVti) = 3.53. Therefore, both
electrons and ions can resonate with LHWs.
The time evolution of the spatial Fourier mode of the electric ?eld in this case
is plotted in Fig. 3.7. In case 3 with both the electron and ion Landau damping,
51
Figure 3.6: Time (normalized to 1/?i) evolution of the electron phase-space distri-
bution in case 2: contour plots of the (a) electron and (b) ion phase-space distribu-
tions and (c) the corresponding electron (black solid lines) and ion (red solid lines)
distribution functions. The ion velocities are normalized to the ion thermal speed.
The black (red) dash line shows the theoretical electron (ion) resonant phase-velocity
vir = 3.264Vti.
52
Figure 3.7: Time evolution of the spatial Fourier mode of the electric ?eld (in the
natural logarithm scale) for case 3 with k?/k = 0.0404, Ti/Te = 4.0, k?e = 0.2255,
and E0 = 0.1.
the overall reduction of the wave amplitude is dominated by the linear ion Landau
damping, whereas the electron Landau e?ects are also seen in the envelop of the
nonlinear oscillations from times ?it = 0.0 to 5.2, as seen in the ?gure. The wave
saturation level is small, similar to that in case 2 shown in Fig. 3.5.
Particle motions in case 3
Fig. 3.8 shows the resulting electron and ion phase-space distributions in case 3
with the same format as Fig. 3.6. Indeed, in such case with both the electron and ion
Landau damping, wave-particle resonance is found for both electrons and ions around
the predicated phase velocities ver = 4.08Vte and vir = 3.53Vti, respectively, as indi-
cated by the vertical dashed lines in Fig. 3.8c. Correspondingly, plateaus are present
in both electron and ion velocity distribution functions around the corresponding
resonant phase velocities.
53
Although the phase-space structure of ions is similar to that of electrons at the
early ?it = 0.22, the vortex structure does not form. The overall ion phase space
contour plots (column b) show structures very similar to those in Fig. 3.6 for case 2.
Again, the ion velocity distribution function curve expands tailward of the resonant
velocity (vir/Vti = 3.53), and ion heating is observed. For the electrons, the resulting
structures are similar to those in Fig. 3.4 for case 1, except that the phase-space vortex
structure becomes more ambiguous in the later time. The reason for this di?erence is
that as the ions keep absorbing energy from the waves, the wave amplitude becomes
smaller and smaller, and thus the trapping of electrons becomes weaker. The plateaus
in the particle distributions, however, do not disappear eventually.
3.5 Nonlinear Damping Rates and Saturation of LHWs
According to the discussions above, when the initial perturbed electric ?eld E0
is small enough, the nonlinear e?ects of the ELD of LHWs are absent, and the waves
will decay into the noise level linearly in logarithm. As E0 gets larger, the nonlinear
ELD e?ects will be involved. The wave will saturate in a BGK equilibrium, and the
damping of the wave is weak, in the nonlinear ELD. Moreover, the ion Landau damp-
ing mainly goes through a linear damping in the resonance with LHWs. Because of
these di?erent characteristics of the electron or/and ion Landau damping for di?erent
initial perturbed electric ?eld, it is interesting for us to research in the damping rate
as a function of E0 for di?erent cases. Meanwhile, for di?erent cases, the saturation
levels of the electric ?eld are quite di?erent. Thus, the saturated electric ?eld, Es, is
also of interest.
3.5.1 Landau damping rate of LHWs
By ?xing k?e = 0.2255 and Ti/Te = 1, the damping rate, ?, as a function of the
initial wave amplitude E0 for various k?/k is shown in Fig. 3.9. For k?/k = 0.0904,
54
Figure 3.8: Resulting particle distributions obtained from case 3 in the same format
as Fig. 3.6.
55
0.066, and 0.0404, all of which are similar to case 1 that only electrons are resonant
with the LHW, there exists a transition from a strong decay at smaller amplitudes to
a weak decay at larger amplitudes. In such cases, three distinct regimes are observed:
(i) a small amplitude, fast damping regime, where only the linear Landau damping
occurs; (ii) a nonlinear damping regime, in which the damping rate decreases as
the amplitude increases, and (iii) a large amplitude, weak decay regime, where the
damping rate saturates to be close to zero. The above results can be understood as
the following. As the electric ?eld becomes larger, the nonlinear e?ects of electron
Landau damping becomes strong, which prevent the wave energy from decaying into
the particles kinetic energy. The damping of the wave is therefore weaker.
Figure 3.9: Damping rate, ?/?i, as a function of the initial wave amplitude E0 in the
cases with k?e = 0.2255 and Ti/Te = 1, for k||/k = 0.0904, 0.066, 0.0404, and 0.001.
On the other hand, for k?/k = 0.001, similar to case 2 with only ion Landau
damping, the damping rate is found to change only slightly with the initial wave
amplitude. The three distinct regimes are not clearly identi?ed, and little nonlinear
56
Landau damping e?ects are observed. Note that the black dot for E0 = 0.1 in Figure
3 corresponds to case 2.
Fig. 3.10 shows the damping rate vs. E0 for k?e = 0.2255 and k?/k = 0.0404,
with Ti/Te = 4.0, 1.0, and 0.33. The electron temperature (and thus the electron
thermal speed) is ?xed in all these cases, while the ion temperature varies with Ti/Te.
In the case with Ti/Te = 1.0, most of the particles that are resonant with the LHWs
are electrons, while few ions are involved in the resonant interaction. The resulting
curve, therefore, shows a trend with the three distinct regimes. For Ti/Te = 0.33,
the ions are colder and thus possess a smaller thermal speed. As a result, there are
almost no ions involved in the interaction. The damping rate ? in this case is nearly
identical to that for Ti/Te = 1.0.
For Ti/Te = 4.0, similar to case 3 in which both the electrons and ions are
resonant with the LHW, the damping rate decreases a bit as the initial amplitude E0
increases, but never reaches a level near zero. These results, again, indicate that the
nonlinear Landau damping e?ects are dominant in the ELD, whereas the ion Landau
damping is of the linear characteristics. As Ti/Te increases, the reduction of the wave
amplitude becomes stronger due to the linear nature of the Landau damping as more
ions participate in the wave-particle resonance.
3.5.2 Saturated Electric Field, Es
Fig. 3.11 shows the saturation level of electric ?eld, Es/E0, as a function ofTi/Te
for cases with an initial E0 = 0.1 and k?/k = 0.066, 0.0404, and 0.001. The electron
temperature is again ?xed. In the cases withk?/k = 0.066, only electrons are resonant
with the LHWs, with?/(k?Vte) = 3.78. The wave decay is dominated by the nonlinear
ELD, and the wave saturation level is found to be nearly constant, at Es/E0?60%
within the plotted temperature range. In the cases with a decreased k?/k = 0.0404,
the LHW frequency is also reduced, but the ratio?/(k?Vte) is increased to 4.08. When
57
Figure 3.10: Damping rate vs. E0 in the cases with k?e = 0.2255 and k?/k = 0.0404,
for Ti/Te = 4.0, 1.0 and 0.33.
Ti/Te?2.0, onlyelectronsareresonantwiththewaves. ThesaturationlevelEs isagain
constant, but at a higher number of 75% due to the larger ?/(k?Vte) and thus smaller
portion of resonant particles. WhenTi/Te>2.0, ions are also resonant with the waves,
which leads to a decreasedEs because of the linear e?ects of the ion Landau damping.
The saturation level Es/E0, however, is ?nite (? 15%), because of the contribution
from the nonlinear Landau damping associated with electrons and trapped ions in the
long time scales. In the cases with k?/k = 0.001, ?/(k?Vte) = 90, and thus electrons
cannot be resonant with the waves. When Ti/Te = 0.25, ?/(kVti) = 7.7, ions also
cannot be resonant with the waves. The level Es/E0 is thus nearly 100% in the
absence of both the electron and ion Landau damping. When Ti/Te>0.25, the ions
resonance condition is satis?ed. The saturation level Es/E0 decreases sharply to a
very small number, near the noise level, due to the purely linear ion Landau damping
e?ects.
58
Figure 3.11: Saturated electric ?eld Es/E0 vs. Ti/Te in the cases with k?e = 0.2255
and k?/k = 0.066, 0.0404,0.001
3.6 Calculation of Driven Current by LHWs
The lower hybrid waves can resonate with both electrons and ions through Lan-
dau damping mechanism, and thus are able to heat the particles and to generate
electric currents in the plasmas. Momentum and energy can be exchanged between
waves and particles obeying the resonance condition, ??k?v, with ? being the fre-
quency, k the wave vector and v the velocity of the particles. We will only consider
the driven currents in the parallel direction to the magnetic ?eld, thus the resonance
condition is ? = k?v?. In this section, we will ?rst review the theory of current drive
in collision plasmas[53, 50, 51, 52]. Then, the results of the driven currents from
the Landau damping of LHWs from our GeFi simulations will be given. However, in
our simulations, the problems are treated as initial value problems and the plasma is
collisionless. Our results are the driven currents from the Landau damping of LHWs
in the collisionless plasmas.
59
3.6.1 Fast Electrons Drive
The theory of driving fast electrons at the tail part of the velocity distribution
function (fast electrons) will be discussed. Although it may be easier to push slow
electrons, it may actually be more e?ective to push fast electrons. In practice, inject-
ing waves with faster parallel phase velocities would be used to deposit momentum in
faster resonant electrons. Although our simulations are in the collisionless plasmas,
our discussions on the current drive here are in the plasmas with collisions in order
to provide a complete concept of current drive.
The Coulomb collision cross section becomes smaller as relative speed between
the colliding particles increases. Thus fast, superthermal electrons collide less often
than slower, thermal electrons. This is because that the average relative speed be-
tween superthermal electrons and most other electrons and ions is much greater than
the relative speed between thermal electrons and most other electrons and ions. In
fact, the ratio of these speeds is roughly v/Vte, where v is the superthermal electron
velocity. Although it may be energetically expensive to accelerate fast electrons in the
?rst place, this energy deposition need to occur less often. But, the advantage is that
currents lasts longer when carried by relatively less collisional electrons. The power
requirements to sustain a given current against collisions can be small. Assume that
the velocity v of an electron is randomized by collisions in a momentum destruction
time of 1/?(v). An increment energy input ?? then produces an incremental current
?j that persists for time 1/?. The parallel momentum absorbed by this electron is
m?v?; the incremental current carried by this electron is ?j = q?v?; and the incre-
mental increase in the electron kinetic energy is ?? = mv??v?. Thus, we could ?nd
the following relationship
?j = ?? qmv
?
. (3.1)
60
The power requirement to refresh this current at time intervals of 1/? is
Pd = ???. (3.2)
Combining Eqs. 4.1 and 4.2 and assuming that the only current is the drive current,
J = ?j, we have the steady-state e?ciency
J
Pd =
q
mv??(v). (3.3)
Apparently, the e?ciency (current per power dissipated) is maximized when the ex-
pression v??(v) is minimized. It is usually identi?ed by Fisch for the utilization of
lower hybrid wave in the limit, v? ? Vte. In this limit, a high e?ciency can be
acquired, since J/Pd ?v2?, with v? being large.
To illustrate the e?ect on the electron distribution function f caused by the
injection of high phase velocity waves, we present the results of a numerical calculation
of our electron Landau damping study of LHWs in Fig. 4.1, which is the plot of
parallel velocity distribution as a function of the parallel velocity (v par is v?/Vte). It
is seen in the ?gure that the plateau (solid line) deviates from the initial distribution
(dashed line). In the current drive problems, the injected waves propagate in only
one direction, thus this plateau will exist in only one side of the velocity distribution
function. This deviation in the distribution forms in the resonant region. Due to the
plateau, the distribution function is asymmetric, indicating the presence of current.
The asymmetry, it turns out, is large enough to signify very large currents.
3.6.2 Electric Currents obtained from GeFi Simulation of LHWs
Lower hybrid wave is a preferable source to generate currents in lots of plasma
devices, i.e., Tokamak. As it can be seen in the GeFi simulation results discussed
in Chapter 3, the parallel velocity distribution function is deviated from the initial
61
Figure 3.12: The plateau in the parallel velocity distribution function in the electron
Landau damping of LHWs (V par is the parallel velocity, v?, normalized in Vte)
Maxwellian distribution in the electron Landau damping of LHWs. Although our
studies are in the collisionless plasmas, we can still calculate the generated current
due to the deviation of the distribution. However, there is not the dissipation or
current drive e?ciency concepts in our study.
In our simulation, the waves are allowed to propagate in both the positive and
negative directions symmetrically, and thus the wave-particle interaction occurs in
both directions. The deviated distribution function, therefore, is nearly symmetric.
In the lower-hybrid drive in fusion plasmas, however, the waves will propagate only in
one direction. The plateau will only exist in one side of the distribution, which is thus
asymmetric. And the asymmetric deviation in the distribution f(ve?) can generate
net parallel current. Here, we estimate the electron parallel current J? generated in
the ve? direction.
62
Let f = fm + ef, where fm = 1/(?2?Vte)exp[?v2?/(2V2te)] is the background
parallel Maxwellian distribution function and ef is the deviation of the distribution
function from the Maxwellian one. Fig. 4.1 shows us one of the typical deviated
distribution functions in our simulations. The current can be calculated as
J? = ?
?
ev? ef(v?)dv?. (3.4)
Since only the plateau (resonance region) of the deviated distribution function con-
tributes to the current, our calculation of J? is only for the plateau region.
Weperformagainasimulationusingtheparametersincase1, whichisintroduced
in the last chapter. In this new run we obtain the electrons? velocities and parallel
velocity distribution information to calculate the parallel current by Eq. 4.18. Fig. 4.2
shows the time evolution ofJ? obtained in case 1, in which only electrons are resonant
with the wave. Here, the current is normalized to en0Vte, and the time is expressed
in units of 1/?i. It is seen that the current reaches a maximum value within a few
wave periods and then keeps oscillating, due to the oscillating wave-particle energies
in the nonlinear ELD.
In order to examine the e?ects of the wave amplitude on the generation of the
currents, Fig. 4.3 presents the resulting J? in the logarithmic scale as a function of
the initial wave amplitude E0. The three curves correspond to k?/k = 0.066 and
Ti/Te = 1.0, k?/k = 0.0404 and Ti/Te = 4.0, and k?/k = 0.0404 and Ti/Te = 1.0,
respectively. In all these cases, k?e = 0.2255. And all the other parameters remain
the same as that given in last chapter, section 3.1. The currents are calculated at
later times when the ?nal BGK state is reached. For k?/k = 0.0404 and Ti/Te = 4.0,
although the linear ion Landau damping is dominant on top of the ELD, the driven
currents show a trend that is similar to that for k?/k = 0.066 and Ti/Te = 1.0 with
ELD only. For all three curves in Fig. 10, the current grows very fast at smaller
63
Figure 3.13: Time evolution of J? (normalized to en0Vte) obtained from case 1.
E0. The growth then slows down as E0 increases. The reason is that when the wave
amplitude is small enough, the linear ELD is dominant, and most of the wave energy
is converted to the particle energy. When E0 is large, the nonlinear ELD e?ects limit
the decay of the wave energy into the particles kinetic energy, and Thus the growth
of the current slows down. For k?/k = 0.066 and Ti/Te = 1.0, and k?/k = 0.0404 and
Ti/Te = 1.0, both dominant with ELD, larger amplitudes of currents are generated
by larger ratio of k?/k (or a larger ?). For the two curves corresponding to the same
ratio of k?/k, the one with a larger temperature ratio Ti/Te possesses ion Landau
damping. The ion Landau damping leads to smaller saturated wave amplitudes, and
thus generates smaller steady-state currents in the electron-wave particle interaction.
Note that the plasmas here are collisionless, and we are unable to calculate the
energy dissipation and the current drive e?ciency. It is aimed to research the electron
Landau damping of LHWs in collisionless plasmas in this thesis. It is seen that there
64
Figure 3.14: Resulting parallel currents as a function of the initial LHW amplitude
E0.
is very close relations between the nonlinear electron Landau damping e?ects with
the current generations. Moreover, the ions? Landau damping can play important
role in the Landau damping of LHWs and thus a?ect the current drive mechanisms.
Further work can about the current drive issue by LHWs be done if we can put the
collisions in the plasmas and do current drive problems with consistent perturbing
LHWs.
3.7 Ions Phase Bunching
In either ?g. 3.6 or ?g. 3.7, it can be seen that the wave amplitude shows a
sudden jump, followed by a sudden damping, at time intervals of ?t? 2?/?i, which
is equal to the ion gyromotion period. The sudden jump and damp is relevant to ions
gyromotion, ions phase bunching in the velocity phase space and ions distribution at
low velocities.
65
Figure 3.15: Scatter plots of ions distribution for selected particles in velocity phase
space (Vix, Viz) at ?it = 0.0, 0.8075, 1.6231, 6.0275, 6.3537, 6.8431.
66
Figure 3.16: Ions velocity distribution function in terms ofVix at ?it = 5.7012, 6.0275,
6.1906, 6.8431.
Figure 3.15 is scatter plots of ions distribution for a group of ions randomly
selected at the beginning of simulation in velocity phase space (Vix,Viz) at ?it = 0.0,
0.8075, 1.6231, 6.0275, 6.3537, 6.8431. It can be seen from the ?gure that the ions
are uniformly distributed in the phase space at ?it = 0.0. At ?it = 0.8075 and
?it = 1.6231, when the wave is damping, ions are highly bunched at a few regions in
the phase space. At ?it = 6.0275, when the amplitude of the wave is small and just
before the sudden jump, ions become uniformly distributed again. During the sudden
jump ?it = 6.3537 and after the judden jump ?it = 6.8431, ions are highly bunched
again at two major regions. Some small bunching structures can also be observed in
those two plots.
67
Figure 3.16 shows plots of the ion velocity distribution function in terms of Vix
at ?it = 5.7012, 6.0275, 6.1906, 6.8431. From the ?gure, evolution of the ion velocity
distribution function from a time before the sudden jump to a time after the sudden
jump can be seen. Initialy at ?it = 5.7012, which is before the sudden jump, the
distribution function is Maxwellian. From plots at ?it = 6.0275 and ?it = 6.1906,
both of which are during the sudden jump, it is seen that the distribution function
gradually changes. The number of ions at low positive velocities increases, while the
number of ions at low negative velocities decreases. At ?it = 6.8431, which is after
the sudden jump, the distribution recovers to Maxwellian.
3.8 Summary
In this chapter, Simulation results of Landau damping of LHWs from the novel
Gyrokinetic electron and Fully kinetic ion (GeFi) code model are presented, including
linear and nonlinear propagations of LHWs and calculations of driven electric current
in the process. The main conclusions are summarized in the following.
(1) From the simulations of the linear electron Landau damping of LHWs, it
is seen the waves decay into the noise level linearly in logarithm. Both the real
frequency and the linear electron Landau damping rate from the simulations show
excellent agreement with that from the analytical theory of electron Landau damping
of LHWs.
(2) Both electrons and ions can resonantly interact with the LHWs. The electrons
are magnetized, and their resonance condition follows ? = k?ve?. On the other hand,
the ions are highly unmagnetized in the magnetic ?eld of LHWs, and their resonance
condition is determined by ? = k?vi.
(3) Trapped electrons are observed correspondingly in the nonlinear electron
Landau damping, as predicated by the nonlinear theory. In the long time nonlinear
68
evolution of the LHWs associated with the electron Landau damping, the amplitude
of the wave becomes oscillatory asymptotically, reaching a ?nal BGK equilibrium.
(4) The ion Landau damping, on the other hand, is dominated by the linear
physics in the LHW time scales, with nearly no trapped ions in the wave-particle
interaction. In the case with solely ion resonance or the case in which there exist
both the ion and the electron resonance, the wave amplitude is signi?cantly reduced
by the ion Landau damping. On the long time scales, however, the ions are still
weakly trapped. Behaviors of magnetized ions appear, with a frequency ? ?i.
(5) As the initial wave amplitude increases, a transition occurs in the electron
Landau damping from a strong linear decay of LHWs to a weak decay, which is
dominated by the nonlinear physics.
(6) Generation of the parallel currents through the ELD from our GeFi simula-
tions is discussed. While the presence of the ion Landau damping results in a smaller
J? than that with solely electron Landau damping, similar trends are observed in the
dependence of the currents on the initial wave amplitudes for cases with or without
the ion Landau damping. The current increases quickly with the initial wave ampli-
tude when the amplitude is small, but slows down when the initial wave amplitude
is large.
69
Chapter 4
GeFi Simulation of Electron-ion Hybrid instability
Kinetic simulation investigation of an instability in a magnetized plasma with a
localized electron cross-?eld ?ow is performed by utilizing GeFi model.
4.1 Introduction
Study of the dynamics that governs the release of free energy associated with
sheared ?ows was given considerable attention in both hydrodynamics[56, 57] and
plasma physics[58, 59], since macroscopic ?ows are commonly encountered in various
plasmas, such as plasmas in tokamak devices, Earth?s magnetotail, plasma sheet
boundary layer, laser-produced plasmas, and Earth?s magnetopause. The shear-
driven instabilities have signi?cant e?ects on particles, momentum, and energy trans-
port. For instance, in tokamak devices, the transition from a low (L) mode to a high
(H) mode of energy and particle con?nement is thought to be excited by the sheared
poloidal ?ows[60, 61]. In the Earth?s magnetopause and plasma sheet boundary layer,
variety of wave activities, that are responsible for the broadband electrostatic noise
observed by satellites, are associated with the induced sheared electron cross-?eld
?ow due to the presence of steep density gradients[62, 68]. In a word, the instabilities
excited by the sheared ?ow play important roles in the wealth of physics activities.
The sheared ?ows can excite di?erent instabilities given di?erent conditions,
such as the low frequency and long wavelength Kelvin-Helmholtz (KH) modes[69].
KH mode can be sustained by a transverse velocity shear for L > ?i, where L is
the velocity shear scale length and ?i is the ion gyroradius. Electron-ion Hybrid
(EIH) instability[70, 71] is another shear-driven instability that can be sustained by
70
the velocity shear ?e < L ? ?i with ?e being the electron gyroradius. Di?erent
from Kelvin Helmholtz mode, EIH instability is a short wavelength (ky?i ? 1 but
ky?e ? 1) and high frequency (?r ??LH, where ?LH is the lower hybrid frequency)
mode.
The linear theory of EIH instability is reviewed in the ?rst place in slab geometry
in uniform plasmas and nonuniform plasmas with kzL = 0, and uniform plasmas
with a ?nite kz ?= 0. Here in the slab geometry, static magnetic ?eld is at z axis,
and electron shear ?ow as a function of x is put at y axis. GeFi kinetic simulations
are then performed in slab geometry and cylindrical geometry with either uniform
plasmas or nonuniform plasmas, with kzL = 0 and kzL ?= 0. In the cylindrical
geometry, static magnetic ?eld is also at z direction, while electron shear ?ow as a
function of radial position r is in ? (poloidal) direction. Results are compared with
the linear theory in a slab geometry. Realistic experimental parameters in Auburn
Linear EXperiment for Instability Studies (ALEXIS) device are adopted in the linear
GeFi kinetic simulations, and the results are compared with the experiments as well.
Nonlinear GeFi kinetic simulations in a slab geometry and uniform plasmas are also
studied, and show that the EIH instability mode ?nally evolves from a short wave
length (kx?i ? 12) mode to a long wave length (kx?i ? 3) mode with frequency??LH
in the nonlinear stage.
4.2 Linear theory of EIH instability
Linear theories of EIH instability in magnetized plasmas in slab geometry with
a sheared electron ?ow channel are reviewed in this section, including general phys-
ical model, dispersion equations in uniform plasmas and nonuniform plasmas with
kzL = 0, and unifrom plasmas with kzL?= 0, and numerical methods for solving the
dispersion equations (shooting code method).
71
4.2.1 General physical model
Consider a schematic con?guration, which consists a magnetized plasma slab in
which a localized electric ?eld is present.[70] The external static magnetic ?eld is
chosen to be directed along z axis and the electric ?eld is in x direction, that is,
B0 = B0(x)ez,
E0 = E0(x)ex, (4.1)
where ex and ez are the unit vectors inxandz directions, respectively. The amplitude
of the electric ?eld E0(x) is assumed to be localized over a region with a width L,
which satis?es ?e < L ? ?i. This con?guration would result in a E0 ? B0 ?ow
velocity in the y direction, whose spatial extent is smaller than ?i. As it is seen, all
the physical quantities vary only along thexaxis. Electron and ion densities can vary
spatially in the whole region of interest, with their scales of variation to be chosen of
the same order of magnitude as that of the electron E0 ?B0 ?ow.
In the case of perpendicular (to static magnetic ?eld) ?ows, an unstable mode
develops whose frequency and growth rate are on the order of lower hybrid wave fre-
quency, ?LH when k? = 0. The perpendicular wave length of the instability mode is
assumed to be much larger than electron Larmor radius, i.e., k??e ? 1. In this limit,
the ?uid description consisting of mass conservation and momentum balance equa-
tions can be adopted. Ions can be assumed to be unmagnetized, since the frequency
and growth rate of the instability mode are larger than ion cyclotron frequency, ?i.
To consider an instability mode with frequency on the order of?LH, cold ?uid plasma
approximation can also be used to determine ion responses. Perpendicular wave num-
ber is also assumed k??i ? 1.
72
The ?uid equations are linearized according to the following normal mode rep-
resentation for the plasma density N and ?ow velocity U :
N = n (x) + bn (x)exp[?i(?t?kyy?k?z)],
U = V (x) +v (x)exp[?i(?t?kyy?k?z)], (4.2)
whereky is the wave number along they axis, and? is the complex angular frequency
of the instability mode. n (x) is the equilibrium plasma density, and V (x) denote
the equilibrium sheared ?ow velocity. ? = i, e denotes the species of the particles.
bn (x) and v (x) are the density and velocity modi?cations due to the perturbation
in the system. To be speci?c, the equilibrium ?ow velocity is chosen as follows:
V (x) = VE(x)ey +Vd (x)ez. (4.3)
Where ey is the unit vector in the y direction; Vd (x) represents the ?ow velocity par-
allel to the equilibrium magnetic ?eld for species ?; and VE(x) denotes the amplitude
of the E0 ?B0 ?ow velocity:
VE(x) = ?cE0(x)B
0(x)
. (4.4)
Assuming the mode is electrostatic, we can obtain the electric ?eld
E = E0(x)ex ??{?(x,?)exp[?i(?t?kyy?k?z)]}, (4.5)
where ?(x,?) is the complex electrostatic potential.
73
From the Poisson?s equation, one can obtain the equilibrium electron density ne
in terms of the equilibrium ion density ni and the equilibrium electric ?eld:
ene(x) = ?
i
Zieni(x)? 14?dE0(x)dx . (4.6)
Here, Zi is the charge number of ions.
The perturbed ?ow velocity v (x) in terms of the electrostatic potential ?(x,?)
can be expressed in the following.
v x = ?( iq m
 D 
)(? d?(x,?)dx ?ky? (x)?(x,?)), (4.7)
v y = ?( q m
 D 
)(? (x)? (x)d?(x,?)dx ?ky? ?(x,?)), (4.8)
v z = ( q m
 ? 
)k??(x,?)?(iV
?
d 
? )v x, (4.9)
where ? (x) denotes the cyclotron frequency of species ? = i, e. The quantities ? ,
D , and ? (x) are obtained by the following expressions.
? (x,?) = ??kyVE(x)?k?Vd (x),
D (x,?) = ?2 (x,?)?? (x)?2 (x),
? (x) = 1 + V
?
E(x)
? (x). (4.10)
The perturbed plasma density bn (x) can be found from the mass conservation
equation, in terms of the particle velocities:
i? bn (x) = ddx[n (x)v x] +in (x)(kyv y +k?v z). (4.11)
74
By linearizing the Poisson?s equation there results in the di?erential dispersion equa-
tion for the instability mode.
d
dx(A(x,?)
d?(x,?)
dx )?q(x,?)?(x,?) = 0, (4.12)
where
q(x,?) = k2yB(x,?) +k2?C(x,?)?kydE(x,?)dx ,
A(x,?) = 1??
 
?2p (x)
D (x,?),
B(x,?) = 1??
 
?2p (x)
D (x,?)(1?
? (x)V?E(x)
?2 (x,?) ),
C(x,?) = 1??
 
?2p (x)
?2 (x,?),
E(x,?) = ?
 
?2p (x)? (x)
D (x,?)? (x,?), (4.13)
and ?p (x) denotes the plasma frequency of species ?. To assume that the ions are
unmagnetized, we can take the limit that ?i = 0 when evaluating Eq. 4.13.
Eq. 4.12 is the wave di?erential dispersion equation in a general form. We will
simplify this equation in the uniform plasma assuming that the ions are stationary in
all of our simpli?cations in the next sections.
Since there is a localized electric ?eld with a shear scale length L>?e, the elec-
tron velocity distribution needs to be modi?ed. The electric ?eld causes a cross-?eld
electron ?ow in they direction. And because the electric ?eld is nonuniform inxdirec-
tion, the distribution function deviates from being a simple drifting Maxwellian[71].
Expand the distribution function in terms of the small parameter ?e/L leads to the
following expression
Fe = ne(Xg)?
?(x)(?Vte)3=2
exp(?w
2
? +v
2
z
V2te ), (4.14)
75
where w2? = v2x + w2y/?(x), wy = vy ?VE(x), VE(x) = ?E(x)/B0 is the sheared
electron cross-?ow, Xg = x+vy/?e is the guiding center, Vte is the electron thermal
speed, and ?(x) = 1 + V?E/?e. This expression is valid up to the second order of
?e/L. The electron density is determined from the quasineutrality condition ne(x) =
ni(x)?(1/e)E?, where e is the electron charge.
4.2.2 EIH instability dispersion equations
The governing wave equation Eq. 4.12 introduced in the last section can be
simpli?ed in a uniform density plasma, given that the ions are stationary, the parallel
wave number k? ? 0, and the electron ?ow is weakly sheared, i.e., that V?E/?e ? 1.
To consider the linear dispersion relation for electrostatic potential ?(x) of lower
hybrid waves, and assuming a ?utelike perturbation, Eq. 4.12 can be expressed
approximately in the following:[71]
d
dx[A(x)
d?(x)
dx ]?k
2
yA(x)?(x) = ?
2( ky?e
??kyVE)Se?(x), (4.15)
where ? is the complex angular frequency of the mode and ky is the wave number in
the y direction.
A(x) = (1 +?2)(1??2LH/?2),
Se = (lnne)? ?V??E/?e, (4.16)
with ? = ?pe/?e and ?2LH = ?2pi?2e/(?2pe + ?2e).
Fortheinstabilitymodeintheuniformplasmas, Eq. 4.15canthenbesimpli?ed:[69]
( d
2
dx2 ?k
2
y +F(?)
kyV??E(x)
??kyVE(x))?(x) = 0, (4.17)
76
where F(?) = ?2/{(1 + ?2)[1 ? (?LH/?)2]}. As it is seen from this equation, the
second derivative of the dc electric ?eld is essential for driving this instability.
For the instability mode in a nonuniform plasma with a density gradient sheared
in the length Ln, one can obtain the dispersion relation equation from Eqs. 4.15 and
4.16.
( d
2
dx2 ?k
2
y +F(?)
ky(V??E ??e/Ln)
??ky?e )?(x) = 0. (4.18)
Ganguli et al.[73] derived the EIH instability dispersion relation for a uniform
plasma with a ?nite kzL?= 0 in slab geometry,
[ d
2
dx2 ?k
2
y ?k
2
z +F(?)(
kyV??E(x)
??kyVE(x) +
k2z?2e
??kyVE(x))]?(x) = 0 (4.19)
4.2.3 Numerical methods for solving the dispersion equation
We have so far discussed about the dispersion relation equation of EIH instability
in both uniform and nonuniform plasmas. In this subsection, a numerical method[70,
72] will be introduced to solve the dispersion equation to obtain the electrostatic
potential ?(x,?) in the region of the ?ow channel and the complex frequency of the
instability mode. The frequency is decomposed into a real frequency and an imaginary
frequency (growth rate):
? = ?R +i?, (4.20)
where ?R denotes the real frequency of the EIH mode, and ? gives the growth rate
(if positive) or damping rate (if negative).
77
The condition, that the electrostatic potential ?(x,?) is damped away from the
region in which the electron ?ow is localized, is imposed to obtain an accurate nu-
merical solution. The region of interest is divided into four intervals:
?(x,?) =
?
????
???
????
????
???
???
?
vr(x,?), x>xr
?(x,?), xm <x<xr
B?(x,?), xl <x<xm
Bvl(x,?), x<xl,
where vr and vl are the WKB solutions that are damped in the intervals x > xr
and x < xl, respectively; ?r and ?l are the numerically obtained solutions in the
intervals xm <x<xr and xl <x<xm, respectively; B is a constant. The function
?l is obtained by integrating from the point xl toward the point xm using the initial
conditions that at the point x = xl, ?l and its spatial derivative are equal to vl
and its spatial derivative, respectively. Similarly, the function ?r is obtained by
integrating from the point xr toward the point xm using the initial conditions that at
the pointx = xr,?r and its spatial derivative are equal tovr and its spatial derivative,
respectively.
The electrostatic potential?(x,?) and its spatial derivative should be continuous
at the point xm. This condition can be used to obtain the complex frequency ? and
the constant B, both of which yield to the following results:
B = ?r(xm,?)?
l(xm,?)
,
F(xm,?) = ddxln(?r(xm,?)?
l(xm,?)
) = 0. (4.21)
78
Note that for a randomly chosen ?, it doesn?t follow that F(xm,?) vanishes
unconditionally, since for a given ?, ?r(x,?) and ?l(x,?) satisfy speci?c initial con-
ditions at the points xr and xl, respectively. Thus, their spatial derivatives are such
that at the point xm, Eq. 4.21 is not automatically satis?ed.
The numerical procedure can be described as follows. Firstly, the density pro?le,
magnetic ?eld, and ?ow velocity are chosen along with values for the quantities xl,
xm and xr, and the initial guess of the frequency ?. The solution to Eq. 4.21 can be
obtained by using the Muller?s method[74]. This results in the desired value of the
frequency and growth rate of the instability. The structure of the electrostatic po-
tential can be obtained for the acquired ?. It is found that the electrostatic potential
and the value of ? are nearly independent of xm as long as it is chosen near the peak
of the localized shear ?ow. In the same sense, the solution is nearly independently
of xr, and xl as long as they lie further than a few skin depth layers away from the
region of localization of the electrostatic potential ?(x,?).
4.3 Theoretical numerical solutions of the EIH instability
Based on the above numerical scheme, a shooting code is developed to solve the
EIH dispersion equations to obtain the eigen functions and eigen values. Results in
both uniform and nonuniform plasmas will be presented.
4.3.1 Numerical Investigation of EIH Instability in Uniform Plasmas
Based on the numerical method discussed above, a shooting code is developed
to solve the dispersion equation in the uniform plasma Eq. 4.17. From the equation,
it is seen that the eigenfunction (?(x,?)) and eigenvalue (?) should depend on three
quantities: ? = ?pe/?e, kyL and ?1 = V0E/L?e. V0E is de?ned in the following way,
VE(x) = V0Ef(x) and f(x) = sech2(x).
79
Fig. 4.1 shows the plot of the typical eigenfunction electrostatic potential?(x,?)
of the EIH instability in the real space. The case used here has ? = 1.0, ?1 = 0.3 and
kyL = 0.5. The eigenvalue frequency can be obtained along with the eigenfunction
as ? = (4.35 + 2.84i)?LH for this case. Note that the black line represents the real
part of ?(x,?), and the red line shows the imaginary part.
Figure 4.1: A typical eigenfunction electrostatic potential (?(x,?)) of EIH instability.
Here ? = 1.0, ?1 = 0.3 and kyL = 0.5. For this case, the frequency ? = (4.35 +
2.84i)?LH. The black line shows the real part of the eigenfunction, and the red line
represents the imaginary part of the eigenfunction.
As it can be seen, the EIH instability in uniform plasmas depends on the three
quantities ?, ?1 and kyL, it is of interest to investigate how the eigenvalue frequency
? changes as one of the three quantities varies.
Fig. 4.2 shows the real frequency and growth rate (imaginary frequency) as a
function of kyL for ? = 0.5, 1.0 and 5.0 with ?1 = 0.3 and mi/me = 1836. It is seen
from the ?gure that both real and imaginary frequency are on the order of ?LH. The
real frequency increases monotonously as the kyL increases, while the growth rate
increases and decreases.
80
Figure 4.2: The real parts (dashed lines) and imaginary parts (solid lines) of the
eigenvalues (?) of Eq. 4.17 as a function of kyL for ? = 0.5 (red lines), 1.0 (black
lines) and 5.0 (blue lines). Here mi/me = 1836 and ?1 = 0.3.
Fig. 4.3 is the plot of real frequency and growth rate versus the quantity ?1 =
V0E/?eL for cases with ? = 0.5 and 1.0. Here, kyL = 0.5 and mi/me = 1836. The
?gure shows that both the real frequency and growth rate increase monotonously and
almost linearly as the quantity ?1 increases.
4.3.2 Numerical investigation of EIH instability in nonuniform plasma
The same numerical method is used to solve the dispersion equation (Eq. 4.18)
of the EIH instability in the nonuniform plasma with a density gradient Ln. Here, in
order to assist our further research with ALEXIS experimental device, the ion-electron
mass ratiomi/me = 73440 used in the device will be adopted in my current numerical
investigations. The new parameter Ln introduced in the dispersion equation Eq. 4.18
will be evaluated together with L by the ratio L/Ln.
Fig. 4.4 shows the plot of the eigenvalue (the complex frequency ?) as a function
of kyL for the case with ? = 1.0, ?1 = 0.1 and L/Ln = 0.5. It is seen from the ?gure
81
Figure 4.3: The real parts (dashed lines) and imaginary parts (solid lines) of the
eigenvalues (?) of Eq. 4.17 as a function of ?1 for ? = 0.5 (red lines), 1.0 (black
lines). Here mi/me = 1836 and kyL = 0.5.
Figure 4.4: The real parts (black solid line) and imaginary parts (red dashed line)
of the eigenvalue (?, ?) of Eq. 4.18 as a function of kyL for ? = 1.0, ?1 = 0.1 and
L/Ln = 0.5
82
Figure 4.5: The real parts (black solid lines) and imaginary parts (red dashed line)
of the eigenvalue (?, ?) of Eq. 4.18 as a function of L/Ln for ? = 1.0, kyL = 0.4 and
?1 = 0.1, L is ?xed as 5.0?e.
that both the real frequency (black solid line) and the growth rate (red dashed line)
vary slightly as kyL changes. In the nonuniform plasma, the frequency and growth
rate of the EIH instability is not sensitive to the values of kyL in our selected cases.
The eigenvalue of Eq. 4.18 as a function of L/Ln is plotted in Fig. 4.5 for the
case with ? = 1.0, ?1 = 0.1 and kyL = 0.5. The black solid line is the real frequency
and the red dashed line is the growth rate. It can be seen from the ?gure that as
the ratio quantity L/Ln increases, both the real frequency and growth rate decreases.
Since as L/Ln approaches zero, the nonuniform cases actually become uniform cases,
we can expect both the frequency and growth rate approach the frequencies in the
corresponding uniform cases.
83
4.4 Linear GeFi simulation of the EIH instability
Linear simulation results of EIH instabilities from linear GeFi model in the elec-
trostatic limit are presented and compared with the results from the numerical solu-
tions of the dispersion equations. The results include kinetic simulations in slab geom-
etry with uniform plasmas and nonuniform plasmas (kzL = 0), and kinetic simulations
in cylindrical geometry with uniform plasmas and nonuniform plasmas (kzL = 0), and
kinetic simulations in slab geometry and uniform plasmas with kzL?= 0. Real experi-
mental parameters in the ALEXIS device are adopted in the kinetic simulations, and
results are compared with the experiments as well.
4.4.1 Simulation model
Initial setup of the GeFi particle simulation in both slab geometry and cylindri-
cal geometry are introduced here along with the particle distribution function and
boundary condition.
In the slab geometry GeFi particle simulation, an external static magnetic ?eld
B0 is in the z direction with unit magnitude, and an electric ?eld as a function
of x, E(x) is in the x direction. Thus the E ? B drift ?ow will be in y direction
and dependent on x. 2-D simulations are performed in x?y plane, thus k? = 0.
The boundary condition in y direction is set to be periodic, whereas the boundary
condition in x direction is a conductor boundary. Assuming that the domain length
in x direction is Lx, and the shear scale length is denoted as L, the electric ?eld takes
this form E(x) = ?1L?esech2((x?Lx/2)/L). Ion-electron mass ratio mi/me = 1836
and light-Alfv`en speed ratio c/VA = 42.8, thus the parameter ? = ?pe/?e = 1.
The quantity ? can be varied by changing c/VA and mi/me, since the relation ? =
(c/VA)/
?
mi/me can be derived.
When the GeFi particle simulation of EIH instability is extended to a cylindrical
geometry (r,?,z) and uniform plasmas. An external static magnetic ?eld B0 is in the
84
z direction with unit magnitude, and an electric ?eld as a function of radial position
r, E(r) is in the radial direction. Thus the E ? B drift motion is in ? direction
and dependent on r. 2-D simulations are performed in r ?? plane, thus k? = 0.
The boundary condition in ? direction is set to be periodic, whereas the boundary
condition in r direction is a conductor boundary. Assuming that the domain length
in r direction is Lr, and the shear scale length is denoted as L, the electric ?eld takes
this form E(r) = ?1L?esech2((r?Lr/2)/L). Ion-electron mass ratio mi/me = 1836
and light-Alfv`en speed ratio c/VA = 42.8 yield ? = 1.0.
Velocity distribution in the GeFi particle simulation for electrons takes the form
of drift Maxwellian distribution,
f(v) = ( me2??T)3=2 exp(?me|v?vd|
2
2?T ), (4.22)
where vd is the E?B drift velocity. While for ions, the velocity distribution is normal
Maxwellian.
4.4.2 EIH instability in a slab geometry and uniform plasma with kz = 0
The case with ? = 1.0, kyL = 0.5 and ?1 = 0.3 is performed in the GeFi particle
simulation in slab geometry and uniform plasmas. In this case, we would like to
compare the eigenfunction and eigenvalue from simulation with that from numerically
solving the dispersion equation. Fig. 4.6 shows the plot of the eigenfunction from
both the simulation and numerical solution. The solid black (red) line shows the
real (imaginary) eigenfunction from the numerical solution, whereas the dashed lines
are from the simulations. As it can be seen from the ?gure, both the real and
imaginary eigenfunction from our GeFi simulation can match well with that from
the numerically solving the dispersion relation. The eigenvalue (frequency) from the
numerical solution is ?t = (4.34,2.84)?LH = (131.54,86.17)?i, and the eigenvalue
85
Figure 4.6: Plot of the eigenfunction of the perturbed electrostatic potential. Solid
black (red) line is the real (imaginary) eigenfunction from numerically solving the
dispersion equation. Dashed black (red) line is the real (imaginary) eigenfunction
from the GeFi simulation.
from our simulation can be analyzed by FFT method to yield ?s = (4.24,2.69)?LH =
(128.30,81.40)?i. The eigenvalue from the GeFi simulation is close to that from the
theory as well. Therefore, the conclusion, that our GeFi simulation of EIH instability
can agree with the theory very well, can be drawn.
Fig. 4.7 is the contour plot of perturbed physical quantities ni (ion density), ne
(electron density), ? (electrostatic potential), E (electric ?eld), Ui (ion ?ow velocity)
and Ue (electron ?ow velocity) in the (x,y) real space for a typical case in the linear
GeFi kinetic simulation in the electrostatic limit. The contour structure of Uey shows
the sheared electron ?ow in y direction at the center of x domain. Structures of ni,
ne, ?, Ex and Ey clearly demonstrate the EIH instability excited in the regime of the
sheared electron ?ow.
Fig. 4.8 is the plot of real frequencies (?) and the growth rates (?) normalized
by lower hybrid frequency (?LH), as a function of kyL, which is the wave vector in
86
Figure 4.7: Contour plot of perturbed physical quantities in (x,y) real space.
87
Figure 4.8: Plot of real frequency (?) and growth rate (?) normalized by lower hybrid
frequency (?LH), as a function of kyL. Dashed black line is the real frequency and
solid red line is the growth rate, both of which are from the theory. The black dots
(red dots) are real frequency (growth rate) from the GeFi simulation
y direction and normalized by the shear scale length L. Here, ? is hold to be 1, and
?1 = 0.1. It can be seen from the ?gure that the dots from either the real frequencies
or the growth rates are on the theoretical lines, indicating that the GeFi kinetic
simulations have excellent agreement with the linear theoretical predications.
Fig. 4.9 is similar to Fig. 4.8. For the same purpose, here we hold ? = 1.0
and kyL = 0.5, but vary the values of ?1 to study the dependence of EIH instability
on the shear ?ow strength. It can also be seen from this ?gure, the real frequencies
(black dots) and growth rates (red dots) from GeFi kinetic simulations have excellent
agreement with the eigenvalues of the dispersion equation.
It can be then summarized that in uniform plasmas, EIH instability can be
excited perfectly in the linear GeFi kinetic simulations for a wealth of parameters. The
example eigenfunction plot (Fig. 4.6) and the plot of real frequencies and growth rates
as functions of kyL (Fig. 4.8) and ?1 (Fig. 4.9) demonstrate very good agreement
88
Figure 4.9: Plot of real frequency (?) and growth rate (?) normalized to lower hybrid
frequency (?LH), as a function of ?1. Dashed black line is the real frequency and
solid red line is the growth rate, both of which are from the theory. The black dots
(red dots) are real frequency (growth rate) from the GeFi simulation
with the eigen functions and eigen values of the dispersion di?erential equation in
slab geometry and uniform plasmas, Eq. 4.17.
4.4.3 EIH instability in a slab geometry and nonuniform plasma with
kz = 0
GeFi kinetic simulations of EIH instability in a slab geometry and nonuniform
plasma with a density gradient Ln are presented and compared with the eigen values
of the corresponding dispersion di?erential equation, Eq. 4.18.
Setup of the sheared electric ?eld is the same as before, while the sheared density
pro?le is embedded in the system as a function of x
n(x) = n0 ? ?n2 [1 + tanh(x?
Lx
2
Ln )]. (4.23)
89
Here Lx is the domain length in x direction, n0 is a constant as the lowest density,
?n is the di?erence of the highest density and the lowest density.
For the purpose of further investigation in the experiments, the ALEXIS device
parameters are adopted in the current simulations. External background magnetic
?eld b0 = 127.5(Gauss), background density n0 = 3.86 ? 108(cm?3), peak value of
external electric ?eldE0 = 20.3(Volt/cm), shear length of electric ?eldL = 0.35(cm),
shear length of density Ln = 0.85(cm), electron temperature Te = 5.33(eV), ion
to electron temperature ratio Ti/Te = 0.001, ion to electron mass ratio mi/me =
73440 (Argon plasma), light speed to Alfv?en speed ratio c/VA=42.8, the beta value
of electron ?e = 0.004.
Fig. 4.10 is the plot of real frequencies and growth rates as a function of kyL
from both simulations and numerical solutions of the dispersion equation. The ?gure
shows that GeFi simulation results have very good agreement with the theoretical
predications. In the same format, in Fig. 4.11 real frequencies and growth rates are
plotted as a function of the ratio L/Ln with L = 0.35cm and kyL = 0.22 being ?xed.
From the ?gure, it can be seen that GeFi simulations have good agreement with the
theory.
4.4.4 EIH instability in cylindrical geometry and uniform plasma with
kz = 0
The investigation of EIH instability is extended to a cylindrical geometry (r,?,z)
and uniform plasma by utilizing GeFi particle simulation. Results are presented in
the following.
Fig. 4.12 is a contour plot of perturbed physical quantities for a typical EIH
simulation in a similar format to Fig. 4.7. In the contour structure in Uey, the
sheared electron ?ow is seen to be in ? direction and dependent on r. From the
contour plots of ne, ?, Ex and Ey, It can be seen that the EIH instability structures
90
Figure 4.10: Plot of real frequency and growth rate as a function of kyL.
Figure 4.11: Plot of real frequency and growth rate as a function of the ratio L/Ln
with L = 0.35cm and kyL = 0.22 being ?xed.
91
Figure 4.12: Contour plot of perturbed physical quantities in (r,?) real space
92
Figure 4.13: Plot of real frequency (?) and growth rate (?) normalized to lower hybrid
frequency (?LH), as a function of k L.
are formed near the center of the radial direction. In ? direction, the mode number
is clearly shown to be ?ve as expected.
It is of interest for us to compare the simulation results in cylindrical geometry
with the theory, although the EIH theory is in a slab geometry. Fig. 4.13 is a plot
of real frequency and growth rate as a function of k L for cases with ? = 1.0 and
?1 = 0.1. It can be seen from the ?gure that real frequencies have excellent agreement
with the dispersion relation, while growth rates are in the same trend as that from
the dispersion theory. The di?erence of the growth rates between simulations and
theory may be due to the cylindrical e?ects in our simulations, whereas the theoretical
calculation has assumed a slab geometry.
93
Figure 4.14: Real frequency and growth rate as a function of k L
4.4.5 EIH instability in cylindrical geometry and nonuniform plasma with
kz = 0
EIH instability in cylindrical geometry and inhomogeneous plasma with a density
gradient Ln in radial direction is studied. The setup of the sheared electric ?eld is
the same as before, while the sheared density pro?le is embedded in the system as a
function of r
n(r) = n0 ? ?n2 [1 + tanh(r?
Lr
2
Ln )]. (4.24)
Here Lr is the domain length in r direction, n0 is a constant as the lowest density, ?n
is the di?erence of the highest density and the lowest density.
Figs. 4.14 and 4.15 are plots of real frequency and growth rate, from measure-
mentsofGeFikineticsimulations, asfunctionsofk L(? = 1.0,?1 = 0.1,L/Ln = 0.41)
and ?1 (? = 1.0,k L = 0.5,L/Ln = 0.41), respectively. It is seen from ?g. 4.14, the
frequency increases as k L increases, while the growth rate increases ?rst to reach a
94
Figure 4.15: Real frequency and growth rate as a function of ?1
peak value, and decrease afterwards. However, both the real frequency and growth
rate increase along with ?1, as shown in ?g. 4.15.
4.4.6 EIH instability with kz ?= 0
It is also of interest to study ?nite k? e?ects on the excitation of EIH instability.
Here, a uniform plasma in cylindrical geometry is used. An external static magnetic
?eld B0 is in the z direction with unit magnitude, and an electric ?eld as a function
of radial positions r, E(r) is in the radial direction. Thus the E?B drift motion is in
? direction and dependent on r. 3-D simulations are performed in r,?,z directions,
thus kz (k?) is not zero. The boundary conditions in ? and z directions are set to be
periodic, whereas the boundary condition in r direction is a conductor boundary.
Realistic parameters in ALEXIS device are adopted in GeFi kinetic simulations:
external magnetic ?eld B0 = 127.5(Gauss), electron temperature Te = 5.33(eV),
background density n0 = 6.0?108(cm?3), peak value of external electric ?eld E0 =
95
Figure 4.16: Growth rates vs kzL.
7.0(V/cm), shear scale length of electric ?eld L = 0.4(cm) and wave number in ?
direction (normalized to L) k L = 0.5. Wave vector in z direction (normalized to L)
kzL varies from 0.0008 to 0.007.
Fig. 4.16 is a plot of growth rates as a function of kzL from theoretical pred-
ications and GeFi kinetic simulations, both of which indicate that the growth rate
decreases as the wave vector kzL increases, and ?nally the growth rate is zero as
kzL becomes greater than certain values. GeFi simulation result is in the same trend
as theory, whereas the theory predicates the growth rate dcreases much faster than
the GeFi simulations. From the theory, critical value of kzL for the growth rate to
become zero is about 0.003, while GeFi kinetic simulation determines this value to
be about 0.006.
Withthesameparametersasgivenabove, anexperimentisperformedinALEXIS
device, and the frequency is measured in the uniform plasma region, as shown in Fig.
4.17. From the ?gure, it is seen that there is not any frequency observed except the
96
Figure 4.17: A measurement from the ALEXIS experiment, in which no instability is
found. The spikes are the noise in the circuit.
noise spikes in the circuit, indicating no instability modes are excited. The reason is
that the length of ALEXIS device isLz = 1.7(m), and thuskzLcan be estimated from
this length to be about 0.0148, which is much greater than either the critical value
from theory 0.003 or that from GeFi simulations 0.006. The result from ALEXIS
experiment is consistent with that from GeFi kinetic simulations.
For the values ofkzLthat can excite EIH instability, it is interesting to investigate
the e?ects of ?nitekzLon the frequencies and growth rates for di?erentk L, as shown
in Fig. 4.18. It can be seen from the ?gure that both frequencies and growth rates
from a ?nite kzL = 0.005 are lower than that from kzL = 0.000. The di?erences
between kzL = 0.005 and kzL = 0.000 in either frequency or growth rate are larger at
smaller k L region and smaller at larger k L region. It can be concluded that ?nite
kzL e?ects decrease either the frequency or the growth rate of EIH instability, and
the e?ects are more signi?cant with small k L than that with larger k L.
4.5 Nonlinear evolution of EIH instability
Nonlinear GeFi particle simulations are then performed to investigate the nonlin-
ear evolution of EIH instability in a slab geometry and uniform plasma with kzL = 0.
97
Figure 4.18: Frequencies and growth rates vs k L for kzL = 0.005 and kzL = 0.000.
Figure 4.19: Contour plot of electrostatic potential in (?,x) space in a nonlinear
simulation.
98
Figure 4.20: Contour plot of electrostatic potential in (?,x) space in a linear simula-
tion.
Figs. 4.19 and 4.20 are contour plots of electrostatic potentials in (?,x) space
for the case with ? = 1.0,?1 = 0.3,kyL = 0.5 in a nonlinear simulation and linear
simulation, respectively. The normalization in x axis is ion gyro radius ?i, while
the normalization in y axis is ion gyro frequency ?i. It is seen from Fig. 4.19
that the nonlinear frequency is about 1.0?LH (?LH = 30.26?i in this case), and the
linear frequency from Fig. 4.20 yields 4.24?LH. The linear frequency structure is at
the center of x domain, where the shear ?ow locates, while the nonlinear frequency
structure is not at the center of x domain. Figure 4.21 shows the frequency FFT
analysis for a nonlinear simulation with ? = 1.0,?1 = 0.1,kyL = 0.7. The peak of
the plot indicates the nonlinear frequency, which is 1.64?LH. The red dashed lines
indicate the linear frequency 2.25?LH. From both cases, the nonlinear EIH frequency
is smaller than the linear frquency.
99
Figure 4.21: Plot of electrostatic potential as a function of ?.
Figure 4.22: Contour plots of physical quantities in real space from a nonlinear sim-
ulation with ? = 1.0,?1 = 0.3,kyL = 0.35 at ?it = 0.0599.
100
Figure 4.23: Contour plots of physical quantities in real space from a nonlinear sim-
ulation with ? = 1.0,?1 = 0.3,kyL = 0.35 at ?it = 0.8497.
Figures 4.22 and 4.23 are contour plots of physical quantities in real space from
a nonlinear simulation with ? = 1.0,?1 = 0.3,kyL = 0.35 at ?it = 0.0599 and
?it = 0.8497, respectively. From the comparison between the two ?gures, it can be
seen that at earlier time ?it = 0.0599, the structures of the EIH instability indicate
a shorter wave length (k?i ? 12), while at later time ?it = 0.8497 the structures
become wider in x direction, indicating a longer wave length (k?i ? 3).
Figure 4.24 is a plot of EIH frequency from nonlinear GeFi particle simulations
as a function of kyL for cases with ? = 1.0,?1 = 0.1. It can be seen from the
?gure that when the wave number kyL is less than about 0.65, the nonlinear EIH
frequency changes slightly, and is around 1.0?LH. When kyL is greater than 0.65, as
kyL increases, the frequency increases and approaches the corresponding linear EIH
frequency.
101
Figure 4.24: Nonlinear EIH frequency as a function ofkyLfor cases with? = 1.0,?1 =
0.1. Black solid line is nonlinear GeFi particle simulation results, and red dashed line
is for linear theory.
4.6 Summary
In this chapter, the electron-ion hybrid (EIH) instability is reviewed in theory
and investigated in our GeFi simulations. With a sheared electric ?eld present in
either uniform or nonuniform plasma, the EIH instability mode can be generated.
The dispersion equation for EIH mode in uniform plasma is introduced as Eq. 4.17,
and the dispersion relation in the nonuniform plasma with a density gradient sheared
in the length Ln is given in Eq. 4.18. A numerical method is introduced to solve
both dispersion equations to obtain the complex eigen function and eigen values. By
this numerical method, a shooting code is programmed in Fortran to numerically
calculate the eigen functions and values of the dispersions. Various parameters are
adopted to obtain the numerical solutions to the dispersion equations, and the results
are presented in section 4.1.4.
102
GeFi model is then used to simulate the EIH instability. The following conclu-
sions can be drawn from GeFi particl simulations.
(1) Results of EIH instability in a uniform plasma and slab geometry from linear
GeFi particle simulations have very good agreement with the linear theory.
(2) In the cylindrical geometry, frequencies from linear GeFi particle simulations
have good agreement with the linear theory for a slab geometry, while the growth
rates are lower than the theoretical values, but in the same trend.
(3) In the uniform plasma in a cylindrical geometry, GeFi particle simulations
indicate that the presence of a ?nite kz can signi?cantly change the threshold of EIH
instability. The results are consistent with ALEXIS experiments.
(4) The nonlinear GeFi particle simulations of the EIH instability in a slab ge-
ometry and uniform plasma shows that in the nonlinear stage, the frequency is lower
than that in the linear simulations. Structures in the real space show that the insta-
bility evolves from a short wavelength mode (k?i ? 12) to a longer wavelength mode
(k?i ? 3). The nonlinear EIH instability deserves a further study.
103
Chapter 5
Nonlinear Parametric Decay Instability of Lower Hybrid Wave
Parametric Instability (PI), a multiple wave modes interaction process, is in-
vestigated in the nonlinear propagation of elestrostatic LHWs by utilization of GeFi
kinetic particle simulation.
5.1 Introduction
Parametric instability of LHWs has been extensively analyzed due to great inter-
est in rf heating of plasmas and current drive in fusion devices, such as tokamaks. Un-
like Landau damping of LHWs, which is a wave-particle interaction process between
LHWs and electrons and/or ions, PI is a wave-wave interaction process among mul-
tiple wave modes. Referring to the literatures on the theory of PI[63, 21, 19, 65, 64],
a pump wave (?0,k0), with frequency ?0 near lower hybrid wave frequency ?LH and
wave vector k0, decays into two high frequency sidebands (?1;2,k1;2), driven by a low
frequency mode (?,k) from various decay channels. Here ?1;2 = ???0,k1;2 = k?k0
are the selection rules determined by momentum and energy conservation of the cou-
pled modes (subscripts 1, 2 refer to the lower and upper sidebands, respectively).
The theory of PI of electrostatic LHWs has been studied in great details. PI of
LHWshasalsobeenobservedinsomelowerhybridcurrentdrive(LHCD)experiments[67,
66]. It has been demonstrated that PI plays important roles in LHCD with spectral
broadening of LHWs in the experiments. However, the experimental study of PI still
has much more to accomplish in the future work. Few kinetic simulation literatures
have been published on PI of LHWs. Computational simulation is one of the most
104
advanced tools to study complicate nonlinear physics. Thus it is of great interest to
investigate PI of LHWs by taking bene?ts of GeFi kinetic particle simulations.
5.2 GeFi Simulation of Parametric Instability of LHWs
Basic idea of the PI simulation is to drive a LHW with (?0,k0) (pump wave)
in a homogeneous plasma in a slab geometry, and observe the time evolution of the
pump LHW. The pump wave decays, and some other wave modes in the system can
be excited and grow up. One of the growing wave modes should be identi?ed as the
low frequency mode with (?,k), and one of the other growing wave modes could be
identi?ed as either the lower sideband (?1,k1) or the upper sideband (?2,k2), with
(?1 = ???0,k1 = k?k0) and (?2 = ?+?0,k2 = k+k0) being the selection rules.
Note that in our previous simulations on the Landau damping of LHWs, the initial
value problem was solved, and thus the LHW was driven only once at the beginning
of the simulation. Here, in the simulations of PI of LHWs, the pump LHW is driven
constantly all through the simulation time.
Simulations are performed in a uniform plasma with an external static mag-
netic ?eld in x direction. Simulation model is complicated since the PI requires 3
dimensional (3-D) simulations in the space and various decay channels can exist. The
pump LHW is placed to propagate in x and z directions, that is k0 = k0xex +k0zez
with k0y = 0. But for the low frequency wave mode or the sidebands, they usually
have ?nite wave vectors in x,y,z directions, k = kxex + kyey + kzez and k1;2 =
k1;2xex +k1;2yey +k1;2zez. Thus 3-D simulations in the space are required. Moreover,
the pump wave is allowed to decay into all kinds of possible wave modes, thus reso-
lutions (grids number) in all three directions should be large enough to exclude the
numerical inaccuracy. The large grid numbers in 3 directions require huge amount of
simulation time, and data ?les? sizes are very large as well. It is di?cult to analyze
enough data ?les to plot time evolutions and frequencies of wave modes.
105
Various decay channels allow not only one decay process to occur in one single
simulation. This makes it even more di?cult to analyze di?erent wave modes in the
system and to identify the sidebands and the low frequency mode. Originally, it is
expected to ?nd either the lower or the upper sideband and the low frequency mode
of the pump wave. However, the sidebands may not be able to be observed, since they
can be subjected to another di?erent PI process, following another decay channel to
dump away before growing up. This kind of similar processes could possibly occur
many times. Finally, several wave modes may grow up, however we may not be able
to identify any of them to be the sidebands.
As discussed in the previous chapters, LHWs can decay due to direct wave-
particle Landau interactions. In order to exclude Landau damping from damping
sources, it is necessary to select a pump wave, which cannot satisfy the Landau
damping condition. This is to make sure that the damping source is only the PI
process. According to the discussions in the above paragraph, although it may be
unable to identify the exact corresponding sidebands and low frequency mode in PI,
the conclusion, that PI process is the only mechanism for the decay of the pump
LHW, can then be drawn.
Mode number m instead of wave number k is used to describe di?erent wave
modes. Let the simulation domain length belfor a 1-D simulation, thenm = k/(2?/l)
is the equation to convert wave number to mode number in simulations. Therefore,
for a pump wave with (?0,m0), and a low frequency mode with (?,m), the lower
sideband and upper sideband should be (?1,m1 = m?m0) and (?2,m2 = m+m0)
as the selection rules. Since GeFi kinetic simulations of PI of LHWs are 3-D, the
mode number has three components m = (mx,my,mz).
106
Figure 5.1: Time evolution of the Fourier electric ?eld of the pump wave in a single
mode simulation.
5.2.1 Single Mode simulation of Pump Wave
As discussed above, it is necessary to exclude the direct wave-particle Landau
damping of the pump LHW from damping sources. A single mode simulation is that
all the other wave modes except the pump wave mode in the simulation system are
?ltered, as in the simulations of Landau damping of LHWs.
The selected pump LHW has wave vector k0?e = (0.004,0.0,0.08746), thus
k0??e = 0.004 and k0??e = 0.08746, where ?e is the electron gyro-radius. From Eq.
1.27, the frequency of this pump wave is calculated to be ?0 = 66.77?i = 2.207?LH.
Ion-electron mass ratio mi/me = 1836, ion-electron temperature ratio Ti/Te = 1.0,
light-Alfv`en speed ratio c/VA = 42.8, external static magnetic ?eld?s magnitude is
set to be 1.0, electron thermal speed Vthe = 1.0. By these parameters, we have
?pi = 42.8?i, ?pe = 1836?i, ?e = 1836?i, and thus ?LH = 30.26?i.
Figure 5.1 shows the plot of time evolution of Fourier electric ?eld of the pump
wave in a single mode simulation. As it is seen from the ?gure that the amplitude
107
of the pump wave grows up gradually due to the constantly driving source. After a
certain time, the amplitude of the wave saturates down to a constant value due to
the nonlinear e?ects. The pump wave does not decay after the saturation. In other
words, there are no or few wave-particle Landau interactions in the propagation of
pump wave. Thus this pump LHW is good for PI simulation to exclude the Landau
damping from the damping sources. Real frequency of the pump wave from GeFi
kinetic simulation is measured to be ?0s = 67.88?i = 2.24?LH, which has a good
agreement with the theoretical predication.
5.2.2 GeFi simulation of PI of LHWs
Case 1
Pump LHW with k0?e = (0.004,0.0,0.08746) is driven constantly in plasmas.
Mode number of the pump wave is set to be m0 = (4,0,4). In y direction, a wave
vector with ky?e = 0.09 is set to have mode number my = 4. In the simulation
system, all wave modes with mode numbers larger than (8,8,8) in any of the three
dimensions are ?ltered. In other words, Only wave modes with |mx|? 8 and |my|? 8
and |mz| ? 8 are allowed to exist in the simulation system. Since in the simulation,
waves can propagate in both positive and negative directions, negative wave modes
are also present. Therefore, there exist totally 17 ? 17 ? 17 = 4913 (zero modes
included) wave modes in the simulation system.
To be di?erent from the simulations of Landau damping of LHWs, where the
initial value problems were solved, here pump LHWs are driven constantly. It is seen
from Fig. 5.1, it takes a while for the pump wave to saturate. Thus the strategy
of the simulation here is the following. Before the saturation of the pump wave,
a single mode simulation is performed. After the pump wave is saturated, other
wave modes with mode numbers less than (8,8,8) are then allowed in the simulation
system. Taking this pump LHW as an example, only the pump wave mode exists in
108
Figure 5.2: Contour plot of the electric ?eld in the (kx,kz) space with the mode
number in y direction my being 0.
the system before ?it = 0.5, and when ?it ? 0.5 all the wave modes smaller than
(8,8,8) are allowed to exist.
In order to ?nd out growing wave modes with the largest amplitudes, contour
plots of electric ?eld in k space are necessary. Since there are three dimensions, mode
numbers in y direction my are ?xed, and contour plots in (kx,kz) space are made.
Wave modes propagating in positive and negative directions are symmetric, thus it
is only necessary to plot the positive my and all (mx,mz) pairs.
Figure 5.2 shows the contour plot of electric ?eld in (kx,kz) space with my = 0.
Vertical (horizontal) dashed lines indicate the mode numbers of mx (mz). For ex-
ample, in the x direction, from the middle kx = 0 (mx = 0) to the right, mx =
1,2,???,7,8 corresponding to the 8 vertical dashed lines, and to the left, mx =
?1,?2,???,?7,?8. It can be seen from the ?gure that there exist two strong dark
points that are symmetric, indicating presence of the pump wave with mode numbers
(4,0,4) (or the symmetric (?4,0,?4), here and through this chapter, we will not
109
Figure 5.3: Contour plot of the electric ?eld in the (kx,kz) space with the mode
number in y direction my being 1. Wave modes with (0,1,0), (4,1,4) and (?4,1,?4)
are found to have the largest amplitudes of electric ?eld.
mention the symmetric mode numbers again). Red color in the ?gure indicates the
presence of other wave modes with much smaller amplitude of electric ?eld. In the
same way, other wave modes with large amplitude of electric ?eld can be found out
by making similar plots for my = 1,2,???,7,8. By doing this, wave modes (0,1,0)
and (4,1,4) are found to have the largest amplitudes of electric ?eld, as shown in Fig
5.3.
Among all contour plots of electric ?eld in (kx,kz) space with di?erent ky and
time, ?g 5.3 is found to have two strongest dark points, indicating that both wave
modes (0,1,0) and (4,1,4) have the largest amplitude of electric ?eld (except the
pump wave). If both modes (0,1,0) and (4,1,4) are considered as the lower hybrid
wave eigen modes, their eigen frequencies can be calculated from the dispersion equa-
tion (Eq. 1.27). The wave mode (0,1,0) with wave vector k1?e = (0,0.0225,0) prop-
agates perpendicular to the magnetic ?eld, and has frequency ? = ?LH ? 30.26?i.
110
Figure 5.4: Time evolution of Fourier electric ?eld of the 3 wave modes: pump wave
(4,0,4) (black line), wave mode (0,1,0) (red line) and wave mode (4,1,4) (green
line).
The wave mode (4,1,4) with wave vector k?e = (0.004,0.0225,0.08746) is calculated
to have frequency ? = 2.145?LH ? 64.92?i. Note that the frequencies calculated
here are LHW eigen frequencies.
After wave modes with the largest growth amplitude are identi?ed, time evolution
of the wave modes are plotted, as shown in Fig. 5.4. The pump wave (4,0,4) is
depicted as black line; the wave mode (0,1,0) is represented by red line; green line
is the wave mode (4,1,4). In the ?gure, it can be seen that the pump wave decays
after the other wave modes are released in the system. The wave modes (0,1,0) and
(4,1,4) grow up ?rst, and then decay along with the pump wave. The whole process
happens very fast within several lower hybrid wave periods. Apparently, wave vectors
of the three wave modes can match up to be identi?ed as parametric decay instability,
that is k1 = k?k0. However, match of the three frequencies is also necessary to be
investigated.
111
Figure 5.5: Frequencies of pump wave (4,0,4), wave mode (0,1,0) and wave mode
(4,1,4) from top to bottom.
Fig. 5.5 is plot of frequencies of pump wave (4,0,4), wave mode (0,1,0) and
wave mode (4,1,4) from top to bottom. The x axis is frequency, and the y axis
is frequency power. From the plot of wave mode (4,0,4), which is the pump wave,
frequency can be measured to be ?0 ? 68.44 = 2.262?LH, close to the theoretical
predication. Frequency of wave mode (0,1,0) can be estimated from the peak of the
middle plot, ?1 ? 30.94 = 1.02?LH, which is very close to its LHW eigen frequency.
Wave mode (0,1,0) is identi?ed to be a lower hybrid wave since its frequency analysis
indicates that it possesses lower hybrid wave eigen frequency. In the bottom plot of
Fig. 5.5, the pulse is not obvious, however the measurement of the largest peak yields
a frequency ? = 101.25?i = 3.346?LH. It can be seen that frequencies of the 3 wave
modes can match up as well, ?1 ????0.
From analysis in k space and frequencies, the following conclusion can be drawn.
A constantly driven pump wave with k0?e = (0.004,0,0.08746) and frequency ?0 ?
2.24?LH decays in to a lower hybrid wave mode (0,1,0) with k1?e = (0,0.0225,0) and
112
frequency?1 ? 1.02?LH, and a wave mode (4,1,4) withk?e = (0.004,0.0225,0.08746)
and frequency ? ? 3.346?LH. This decay process is identi?ed as a parametric de-
cay instability from the aspect of selection rules of wave vectors k1 = k ? k0 and
frequencies ?1 = ???0.
Inthe bottomplot of Fig. 5.5, the largest peak ismeasured as the frequency of the
wave mode?, while at the left of the peak, there is a broad frequency spectrum in the
range [36,101]?i. Explanation of this smaller frequency is the following. In this broad
spectrum, there might be a frequency as the lower hybrid wave eigen frequency of wave
mode (4,1,4), which is about 64.92?i. In the parametric instability simulation and
with an ideal condition, the eigen frequency of mode (4,1,4) should not be observed
in the analysis. However, due to limited simulation time and limited datas used in
FFT and the fact that wave mode (4,1,4) grows up and damps away within several
lower hybrid wave periods, the decay wave mode frequency (? = 3.346?LH) cannot
obtain so large a power to completely dominate the eigen frequency of the wave mode
in the analysis. In the broad spectrum there is another frequency which should be
about 36.94?i. This frequency is also from parametric instability process with the
same mode (0,1,0), but the negative frequency of this mode. It can be seen in the
frequency plot of mode (0,1,0) in Fig. 5.5, the wave mode has a negative frequency
about ??1 = ?30.26?i. By the selection rules of frequency, ??1 = ???0, wave mode
(4,1,4) should have a frequency ? = 36.94?i. Thus, wave mode (4,1,4) might have 3
frequencies 36.94?i, 64.92?i and 101.25?i. Since the time range of analysis is short,
less than 1/?i, the resolution of frequency analysis is about 10?i. Thus it is not
di?cult to understand why the frequency pulses in the frequency spectrum of wave
mode (4,1,4) are not well formed.
Thedecaychannelofthiscaseisalsoidenti?ed. Inthissimulationcase,k?Vte/? =
0.13 < 0.3, thus it falls in ?Case A? discussed in Secthion 1.4.2. And ?0 = 2.26?LH >
2?LH, thus the decay channel is identi?ed as that a pump wave decays into two lower
113
Figure 5.6: Time evolution of Fourier electric ?eld of the 3 wave modes: pump wave
(4,0,4) (black line), wave mode (0,1,0) (red line) and wave mode (?4,1,?4) (green
line).
hybrid waves. From Eq. 1.73, the growth rate in the absence of damping can be
estimated to obtain ?0t/?i ? 27.71. Growth rate in the simulation case is estimated
to be ?0s/?i ? 33.61, which has a good agreement with the theory estimation.
Besides wave modes (0,1,0) and (4,1,4), it is found that there also exist some
other wave modes with very large amplitude of electric ?eld in the simulation sys-
tem. Take a look at Fig. 5.3 again, there is another highlight black spot at mode
(?4,1,?4). Note that (?4,1,?4) is not the symetric mode of (4,1,4) (its symmetric
mode is (?4,?1,?4)). Time evolution of Fourier electric ?eld of wave modes (4,0,4),
(0,1,0) and (?4,1,?4) are plotted in Fig. 5.6 and their respective frequency spec-
trum is plotted in Fig. 5.7. Denote (?0,k0), (?1,k1) and (?,k) as wave modes (4,0,4),
(0,1,0) and (-4,1,-4), respectively. The selection rules of frequency and wave vectors
are ? = ?1 ??0 and k = k1 ?k0. Rest of the analysis is similar to that above.
114
Figure 5.7: Frequencies of pump wave (4,0,4), wave mode (0,1,0) and wave mode
(?4,1,?4) from top to bottom.
Analysis of wave modes in k space also shows some other wave modes with a
little smaller but still large amplitudes of electric ?eld. Fig. 5.8 shows the contour
plot in (kx,kz) space with mode number my = 2. Three highlight spot are found
in the ?gure at wave modes (0,2,0), (4,2,4) and (?4,2,?4). Time evolution and
frequency spectrum of these 3 wave modes with the pump wave are plotted in Fig.
5.9 and Fig. 5.10, respectively. Analysis and discussion of these modes by selection
rules in frequencies and wave vectors are the same as that above.
Case 2
In case 2, pump LHW has wave vector k0?e = (0.008,0.000,0.08746), thus the
eigen frequency of this wave is?0 = 122.34?i = 4.045?LH from the dispersion relation
of LHWs. The pump wave is also set up at wave mode (4,0,4). All the other wave
modes are ?ltered until t?i ? 0.5, when the pump wave is about saturated, wave
115
Figure 5.8: Contour plot of the electric ?eld in (kx,kz) space with mode number in
y direction my being 2. Wave modes with (0,2,0), (4,2,4) and (?4,2,?4) are found
to have the largest amplitudes of electric ?eld.
Figure 5.9: Time evolution of Fourier electric ?eld of the 4 wave modes: pump wave
(4,0,4) (black line), wave mode (0,2,0) (red line), wave mode (4,2,4) (yellow line)
and wave mode (?4,2,?4) (green line).
116
Figure 5.10: Frequencies of pump wave (4,0,4), wave mode (0,2,0), wave mode
(4,2,4) and wave mode (?4,2,?4) from top to bottom.
modes with mode number less than (8,8,8) are then allowed to exist in the simulation
system.
Analysis of Fourier electric ?eld in (kx,kz) space shows that wave modes (0,3,1),
(4,3,5) and (?4,3,?3) have large amplitude as the highlight spots in Fig. 5.11. Since
electric ?elds of wave modes are oscillating, they may not be observed at a time due to
small value at the trough of the oscillations, but they do exist in the simulation system
and can be observed some other times when they reach the crest of the oscillations.
Figs. 5.12 and 5.13 are contour plots in (kx,kz) space at ?it = 0.58824 and ?it =
0.54466, respectively. In the ?gures, wave modes (1,3,0), (?1,3,0), (5,3,4) and
(?5,3,?4) are found to have large amplitudes.
Time evolutions and frequency spectrums of these wave modes are then plotted.
Figs. 5.14 and 5.15 are plots of time evolution and frequency spectrum of modes
(4,0,4), (0,3,1), (4,3,5) and (?4,3,?3). Figs. 5.16 and 5.17 are plots of time
evolution and frequency of modes (4,0,4), (1,3,0) and 5,3,4. Figs. 5.18 and 5.19 are
117
Figure 5.11: Contour plot of the electric ?eld in (kx,kz) space with mode number in
y direction my = 3. Wave modes with (0,3,1), (4,3,5) and (?4,3,?3) are found to
have the largest amplitudes of electric ?eld.
Figure 5.12: Contour plot of the electric ?eld in (kx,kz) space with mode number in
y direction my = 3 at ?it = 0.58824. Wave modes (1,3,0) and (?1,3,0) are found
to have large amplitudes of electric ?eld.
118
Figure 5.13: Contour plot of the electric ?eld in (kx,kz) space with mode number in
y direction my = 3 at ?it = 0.54466. Wave modes with (5,3,3) and (?5,3,?3) are
found to have large amplitudes of electric ?eld.
Figure 5.14: Time evolution of Fourier electric ?eld of the 4 wave modes: pump wave
(4,0,4) (black line), wave mode (0,3,1) (red line), wave mode (4,3,5) (yellow line)
and wave mode (?4,3,?3) (green line).
119
Figure 5.15: Frequencies of pump wave (4,0,4), wave mode (0,3,1), wave mode
(4,3,5) and wave mode (?4,3,?3) from top to bottom.
similar plots for modes (4,0,4), (?1,3,0) and?5,3,?4. Analysis of these wave modes
by the selection rules of wave vectors and frequencies, for the purpose of momentum
and energy conservation, are similar to that in case 1. Some wave modes, such as
(4,3,5), (5,3,4) and (?5,3,?4), have small frequency power. However in their time
evolution plots, it can be seen that their peak values of electric ?eld amplitude are
as large as other wave modes with much larger frequency power. The explanation for
this is that although wave modes (4,3,5), (5,3,4) and (?5,3,?4) have large peak
values of electric ?eld amplitude, they grow up and damp away very fast, within
about 0.1/?i. The short time existence of those wave modes causes low frequency
power in their frequency spectrums.
5.2.3 Particles distribution in PI simulation
Although wave-wave interaction plays the most important role in the PI process
of LHWs, it is also of great interest to study electons and ions distributions in the
120
Figure 5.16: Time evolution of Fourier electric ?eld of 3 wave modes: pump wave
(4,0,4) (black line), wave mode (1,3,0) (red line), and wave mode (5,3,4) (green
line).
Figure 5.17: Time evolution of Fourier electric ?eld of 3 wave modes: pump wave
(4,0,4) (black line), wave mode (?1,3,0) (red line), and wave mode (?5,3,?4)
(green line).
121
Figure 5.18: Frequencies of pump wave (4,0,4), wave mode (1,3,0) and wave mode
(5,3,4) from top to bottom.
Figure 5.19: Frequencies of pump wave (4,0,4), wave mode (?1,3,0) and wave mode
(?5,3,?4) from top to bottom.
122
Figure 5.20: Time evolution of electron and ion distributions in Case 1 for ?it =
0.0218 (top), ?it = 0.5229 (middle), ?it = 0.6972 (bottom). Column (a) is electron
distribution in phase space (Ve?,x). (b) is ion distribution in phase space (Vi,x). (c)
are plots of electron (solid lines) and ion (red lines) distribution functions.
process. Wave-particle interactions change velocity distributions of particles, as pre-
sented in the Landau damping of LHWs. The change of the velocity distributions
of particles may in turn modify the characteristics of wave modes in the simulation
system. Thus it is necessary and interesting to investigate velocity distributions of
particles along with properties of wave modes.
Both electron and ion distributions in Case 1 are investigated in phase space
and their respective velocity distribution functions at di?erent time in Fig. 5.20.
In the ?gure, column (a) and column (b) are contour plots of electron distribution
in phase space (Ve?,x) and contour plots of ion distribution in phase space (Vi,x),
respectively. Ion (solid lines) and electron (red lines) velocity distribution functions
are plotted in column (c). Electron velocity is normalized to electron thermal speed
Vthe, and ion velocity is normalized to ion thermal speed Vthi. From top to bottom,
?it = 0.0218, ?it = 0.5229 and ?it = 0.6972. Note that in Case 1, ?it = 0.5 is
123
the time when the simulation system is released, and other wave modes besides the
pump wave are allowed to exist in the system. It can be seen from column (a) that
electrons are heated at time ?it = 0.5229 and ?it = 0.6972. Lots of electrons move
to higher velocity region, especially electrons with negative velocities. In column (b),
it is seen that the blue structure does not expand to higher velocity region, thus ions
are not heated. In velocity distribution functions in column (c), electron distribution
function, depicted as black lines, expands at the tail, especially in negative velocity
part, while ion distribution function is unchanged as seen from red curves. This
change in electron parallel velocity distribution may in turn a?ect properties of the
wave modes. For example, the expected low frequency modes ?,k in either Case 1
or Case 2 have large frequencies rather than low frequency (? ?i). Large frequencies
of the expected low frequency modes may be excited by the alteration of electron
distribution.
5.3 Summary
In this chapter, parametric instability of LHWs is studied by utilization of GeFi
kinetic simulation. A constantly driven source in a uniform plasma is used to gen-
erate a pump LHW with wave vector k0 and frequency ?0. Single mode simulation
shows that the pump LHW saturates at a constant amplitude after some time. In
PI simulations, the pump LHW is allowed to grow to saturation as a single mode
before other wave modes, with mode numbers smaller than (8,8,8), are allowed to
exist in simulation system. In either Case 1 or Case 2, the pump LHW is found to
decay into a lower sideband lower hybrid wave (k1,?1) and an expected low frequency
wave mode (k,?). The three waves satisfy the selection rule in both wave vectors
and frequencies, k1 = k?k0 and ?1 = ???0, or k = k1 ?k0 and ? = ?1 ??0,
for the purpose of momentum and energy conservation. Results are concluded in the
following.
124
(1) PI process is observed in the nonlinear propagaton of LHWs in a magnetized
uniform plasma in GeFi kinetic simulations. The process is found to be very fast
within several lower hybrid wave periods.
(2) More than one PI process usually occur simultaneously. Di?erent decay
channels are possible to excite wave modes in di?erent PI processes as shown in both
Case 1 and Case 2.
(3) Electron is heated and the electron velocity distribution is modi?ed signi?-
cantly, while ion is not heated and the ion velocity distribution is unchanged. The
alteration of electron parallel velocity distribution has signi?cant e?ects on low fre-
quency wave modes.
5.4 Future work
In order to understand the parametric instability of LHWs in details, more work
needs to be accomplished. Since LHWs can interact with both electrons and ions,
it is necessary to decouple the electron nonlinear e?ects and ion nonlinear e?ects to
understand their respective contributions to the parametric instability. Speci?cally,
since the electron is magnetized, the nonlinear coupling in such a system is dominated
by a 3-D physics through the E ? B drift velocity. The ions, on the other hand,
is nearly unmagnetized. Its dynamics may be predominatly in a 2-D fashion. By
systematically singling out the electron and ion nonlinear dynamics in both 2-D and
3-D simulations, we could be able to limit the decay channels to understand the 3-D
nonlinear physics. GeFi particle simulations with linear electrons and noninear ions,
as well as nonlinear electrons and linear ions are planned to perform in the future
work. Based on these analysis, more detailed simulation of both nonlinear electrons
and nonlinear ions can be performed, and the physics of the parametric instability of
LHWs can be understood.
125
Chapter 6
Summary
In this thesis, propagation of LHWs as the main topic is investigated in great
details in two aspects: wave-particle interactions (Landau damping) and wave-wave
interactions (parametric instability) of LHWs. As another application of our GeFi
model, the excitation of a shear ?ow driven instability, electron-ion hybrid instability,
is also investigated and discussed.
In chapter 1, physics of Landau damping and parametric instability of lower hy-
brid waves are reviewed. The theory of Landau damping, including both linear and
nonlinear Landau damping rates of Langmuir waves, is introduced. Main characteris-
tics of nonlinear Landau damping, such as particles trapping and long time evolution
into a ?nal BGK equilibrium, etc, are discussed by referring to some recent published
literatures. For lower hybrid wave, the electromagnetic dispersion relation as well as
the electrostatic dispersion relation in both cold and warm plasmas is introduced in
great details. The coe?cient of linear electron Landau damping of LHWs is derived
given the condition that the damping rate is much smaller than the real frequency.
The important property of Landau damping of LHWs, that unlike most other wave
modes LHWs can interact with both electrons and ions directly, is ?nally discussed.
In the last section, theory of PI is reviewed including both electron and ion nonlinear
response e?ects. The nonlinear dispersion relation equation of parametric instability
of LHWs is given as Eq. 1.70. Based on the theory and dispersion equation, various
decay channels are categorized according to the quantity k?ve/?. For various possible
decay channels, growth rates are also estimated.
126
Original novel GeFi plasma simulation model is introduced in chapter 2. This
new model, in which the electrons are treated as GK particles and ions are treated as
FK particles, is particularly applicable to problems in which wave modes ranging from
magnetosonic and Alfv?en waves to lower-hybrid/whistler waves need to be handled
on an equal footing. To utilize this code, the simulated physical processes should
be dominated by wave frequencies ? < ?e, and wave numbers k? < k?. With
fast electron gyromotion and Langmuir oscillations removed from the dynamics, the
GeFi model can readily employ realistic mi/me mass ratio. As discussed above, the
GeFi model has been improved and modi?ed to allow the existence of modes with
wavelengths on the same scale of the background nonuniformity.
In chapter 3, simulation results of Landau damping of LHWs from our novel
Gyrokinetic electron and Fully kinetic ion (GeFi) model, including linear and nonlin-
ear propagations of LHWs and calculations of driven electric currents in the process,
are presented. Main conclusions are summarized in the following. From the simu-
lations of the linear electron Landau damping of LHWs, it is seen the waves decay
into the noise level linearly in logarithm. Both the real frequency and the linear elec-
tron Landau damping rate from the simulations show excellent agreement with that
from the analytical theory of electron Landau damping of LHWs. Both electrons and
ions can resonantly interact with the LHWs. The electrons are magnetized, and their
resonance condition follows? = k?ve?. On the other hand, the ions are highly unmag-
netized in the magnetic ?eld of LHWs, and their resonance condition is determined by
? = k?vi. Trapped electrons are observed correspondingly in the nonlinear electron
Landau damping, as predicated by the nonlinear theory. In the long time nonlinear
evolution of the LHWs associated with the electron Landau damping, the amplitude
of the wave becomes oscillatory asymptotically, reaching a ?nal BGK equilibrium.
The ion Landau damping, on the other hand, is dominated by the linear physics in
the LHW time scales, with nearly no trapped ions in the wave-particle interaction.
127
In the case with solely ion resonance or the case in which there exist both the ion
and the electron resonance, the wave amplitude is signi?cantly reduced by the ion
Landau damping. On the long time scales, however, the ions are still weakly trapped.
Behaviors of magnetized ions appear, with a frequency ? ?i. As the initial wave am-
plitude increases, a transition occurs in the electron Landau damping from a strong
linear decay of LHWs to a weak decay, which is dominated by the nonlinear physics.
Generation of the parallel currents through the ELD from our GeFi simulations is
discussed. While the presence of the ion Landau damping results in a smaller J? than
that with solely electron Landau damping, similar trends are observed in the depen-
dence of the currents on the initial wave amplitudes for cases with or without the ion
Landau damping. The current increases quickly with the initial wave amplitude when
the amplitude is small, but slows down when the initial wave amplitude is large.
Investigation of electron-ion hybrid (EIH) instability is presented in chapter 4.
With a sheared electron ?ow in either uniform or nonuniform plasma, the EIH insta-
bility mode can be excited. The dispersion equation for EIH mode in uniform plasma
is introduced as Eq. 4.17, and the dispersion relation in nonuniform plasma with a
density gradient sheared in the length Ln is given in Eq. 4.18. A numerical method
is introduced to solve both dispersion equations to obtain the complex eigen func-
tions and eigen values. By this numerical method, a shooting code is programmed
in Fortran to numerically calculate the eigen functions and values of the dispersion
equations. Various parameters are adopted to obtain the numerical solutions to the
dispersion equations, and the results are presented in section 4.1.4. GeFi particle
simulation is then used to simulate the EIH instability. Results of EIH instability in
a uniform plasma and slab geometry from GeFi particle simulations have very good
agreement with the theory. In the cylindrical geometry, frequencies from GeFi parti-
cle simulations have good agreement with the theory for a slab geometry, while the
growth rates are in the same trend. In the uniform plasma in a cylindrical geometry,
128
with a ?nite kz, GeFi particle simulations indicate that the presence of a ?nite kz can
signi?cantly change the threshold of EIH instability. The results are consistent with
ALEXIS experiments. In the nonlinear GeFi particle simulations of EIH instability
in a slab geometry and uniform plasma, the frequency is lower than that in linear
simulations. Structures in the real space show that the instability evolves from a short
wavelength mode (k?i ? 12) to a longer wavelength mode (k?i ? 3). The nonlinear
EIH instability deserves a further study.
Parametric instability of LHWs is discussed in chapter 5. Our GeFi kinetic
scheme is used to simulate the process of the parametric decay of LHWs. A constantly
driving source in the plasma is used to generate a pump lower hybrid wave with
selected wave vector k0 and frequency ?0. The single mode simulation shows that
the pump wave saturates at a constant amplitude after a certain time. The pump
lower hybrid wave is allowed to grow to the saturation level as a single mode before the
other wave modes with mode numbers smaller than (8,8,8) are released in the system.
After the other wave modes are allowed to exist in the system, the pump LHW is
found to decay into a lower sideband lower hybrid wave (k1,?1) and a wave mode
(k,?). The three wave modes follow the selection rules of energy and momentum,
k1 = k?k0 and?1 = ???0. The growth rates estimated from the simulation is close
to that estimated from the theory. The main conclusions can be drawn as following.
PI process is observed in the nonlinear propagaton of LHWs in a magnetized uniform
plasma in GeFi kinetic simulations. The process is found to be very fast within several
lower hybrid wave periods. More than one PI process usually occur simultaneously.
Di?erent decay channels are possible to excite wave modes in di?erent PI processes
as shown in both Case 1 and Case 2. Electron is heated and the electron velocity
distribution is modi?ed signi?cantly, while ion is not heated and the ion velocity
distribution is unchanged. The alteration of electron parallel velocity distribution
has signi?cant e?ects on low frequency wave modes.
129
Bibliography
[1] L. D. Landau, J. Phys.(Moscow), 10, 25 (1946).
[2] T. M. O?Neil, Phys. Fluid, 8, 2255 (1965).
[3] N. J. Fisch, Theory of RF Current-Drive, Reviews of Modern Physics, 59, 175
(1987).
[4] I.B. Bernstein, J.M. Greene, and M.D. Kruskal, Phys. Rev. 108, 546 (1957)
[5] G. Manfredi, Phys. Rev. Lett., 79, 2815 (1997).
[6] M. Brunetti, F. Califano, and F. Pegoraro, Phys. Rev E, 62, 4109 (2000).
[7] J. R. Danielson, F. Anderegg, and C. F. Driscoll, Phys. Rev. Lett., 92, 245003
(2004).
[8] J.H. Malmberg and C.B. Wharton, Phys. Rev. Lett. 6, 184 (1964).
[9] R.R.J. Gagn?e and M.M. Shoucri, J. Comput. Phys. 24, 445 (1977).
[10] D.R. Nicholson, Introduction to Plasma Theory (Wiley, New York, 1983).
[11] K. Miyamoto, Plasma Physics for Nuclear Fusion (MIT Press, cambridge, MA,
1987).
[12] G. Brodin, Phys. Rev. Lett. 7, 78, 1997
[13] C.B. Wharton, J.H. Malmberg, and T.M. O?Neil, Phys. Fluids, 8, 17, 1968
[14] J.H. Malmberg and C.B. Wharton, Phys. Rev. Lett. 19, 775 (1967).
[15] N.J. Fisch, Theory of RF current-Drive, Reviews of Modern Physics, 59, 175
(1987).
[16] T. H. Stix, Waves in Plasmas, (AIP, New York, 1992).
[17] P. Bonoli, IEEE Trans. Plasma Sci. 12, 95, (1984).
[18] M. Brambilla, Plasma. Phys. 18, 669, (1976).
[19] M. Porkolab, Phys. Fluids, 17, 1432 (1974).
[20] M. Porkolab, Phys. Fluids 20, 2058, (1977).
130
[21] V. K. Tripathi, C. S. Liu, and C. Grebogi, Phys. Fluids 22, 301, (1979).
[22] P. M. Bellan and M. Porkolab, Phys. Fluids 19, 995, (1976).
[23] M. Porkolab, et al., Phys. Rev. Lett., 53, 1229, (1984).
[24] J. J. Schuss, Phys. Fluids 18, 1178, (1975).
[25] R. E. Bell, S. Bernabei, A. Cavallo, et al., Phys. Rev. Lett. 60, 1294, (1988).
[26] T.H. Stix, Phys. Rev. Lett. 15, 878 (1965)
[27] A.L. Verdon, Iver H. Cairns, D.B. Melrose, and P.A. Robinson, Phys. Plasmas,
16, 052105 (2009)
[28] Y.A. Omer?chenko, R.Z. Sagdeev, V.D. Shapiro, and V.I. Shevchenko, Sov. J.
Plasma Phys. 15, 427 (1989)
[29] R. Bingham, J.M. Dawson, and V.D. Shapiro, J. Plasma Phys. 68, 161 (2002)
[30] J.B. McBride, E. Ott, J.P. Boris, and J.H. Orens, Phys. Fluids 15, 2367 (1972)
[31] J.F. Drake, J.D. Huba, and N.T. Gladd, phys. Fluids 26, 2247 (1983)
[32] Y.A. Omel?chenko, R.Z. Sagdeev, V.D. Shapiro, and V.I. Shevchenko, Sov. J.
Plasma Phys. 15, 427 (1989)
[33] I.H. Cairns, Publ. -Astron. Soc. Aust. 18, 336 (2001)
[34] B.F. McMillan and I.H. Cains, Phys. Plasmas 14, 012103 (2007)
[35] Y. Lin, X.Y. Wang, Z. Lin and L. Chen, Plasma Phys. Controled. Fusion, 47,
657-669 (2005)
[36] Y. Lin, x. Y. Wang, L. Chen, X. Lu and W. Kong, Plasmas Phys. Control.
Fusion, 53, 054013 (2011)
[37] Hoshino M, J. Geophys. Res. 92, 7368 (1987)
[38] P.L. Pritchett, J. Geophys. Res. 106, 3783 (2001)
[39] H. J. Cai and Lee. L.C. Phys. Plasmas, 4, 590 (1997)
[40] M.E. Mandt, E.R. Denton and J.F. Drake, Geophys. Res. Lett. 21, 73 (1994)
[41] Y. Lin and D. W. Swift, J. Geophys. Res. 101, 19, 859 (1996)
[42] M. Nakamura and M. Scholer, J. Geophys. Res. 105, 23179 (2000)
[43] Frieman E. A. and Chen L., Phys. Fluids,25, 502 (1982)
[44] Brizard A., J. Plasma Phys. 41, 541 (1989)
131
[45] S. Wang, Phys. Rev. E 64, 056404, (2001)
[46] Chen, L., Waves and Instabilities of Plasmas, World Scienti c, 1987, Chap3
[47] L. Qi and S. Wang, Phys. Plasmas 16, 062504, (2009)
[48] Lee W.W., J. Comput. Phys., 72, 243 (1987)
[49] Hahm T.S., Lee W.W. and Brizard A., Phys. Fluids, 31, 1940 (1988)
[50] N. J. Fisch, Theory of RF Current-Drive, Reviews of Modern Physics, 59, 175
(1987).
[51] Porkolab, M., 1985a, in Wave Heating and Current Drive in Plasmas, edited by
V. L. Granatstein and P.L. Colestock (Gordon and Breach, New York), p. 219.
[52] Charles F. F. Karney and Nathaniel J. Fisch, Phys. Fluids 28, 116 (1985)
[53] C. F. F. Karney, N. J. Fisch, and F. C. Jobes, Phys. Rev. A 32, 2554C2556
(1985)
[54] Landau, L., Phys. Z. Sowjetunion, 10, 154 (1936)
[55] Vedenov, A.A., in Reviews of PLasma Physics, edited by M.A. Leontovich (Con-
sultants Bureau, New York), Vol. 3, p. 229 (1967)
[56] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability (Dover, New
York, 1981), Chap. XI.
[57] P.G. Drazin and L.N. Howard, in Advances in Applied Mechanics (Academic,
New York, 1966), Vol. 9, Chap. 1.
[58] A.B. Mikhailovskii, Theory of Plasma Instabilities (Consultants Bureau, New
York, 1974), Vol. 1, Chap. 1.
[59] A.B. Mikhailovskii, Theory of Plasma Instabilities (Consultants Bureau, New
York, 1974), Vol. 2, Chap. 1 and 8.
[60] K.C. Shaing and E.C. Crume, Jr., Phys. Rev. Lett. 63, 2369 (1989).
[61] R.J. Groebner, K.H. Burrell, and R.P. Seraydarian, Phys. Rev. Lett. 64, 3015
(1990).
[62] P. Song, R.C. Elphic, C.T. Russell, J.T. Gosling, and C.A. Cattell, J. Geophys.
Res. 95, 6375 (1990).
[63] V.K. Tripathi, C. Grebogi, and C.S. Liu, Physics of Fluids, Vol. 20, No. 9 (1977)
[64] L. Chen, R.L. Berger, Nuclear Fusion 17, 4 (1977)
[65] R.L. Berger and F.W. Perkins, Physics of Fluids, 19, 3 (1976)
132
[66] R.Cesario, et al Nucl. Fusion 46 (2006) 462-476
[67] R. Cesario, L. Amicucci, D. De Areangelis, M. Ferrari, F. Napoli, E. Pullara,
Journal of Physics: Conference Series 260, 012008 (2010)
[68] H. Romero, G. Ganguli, P. Palmadesso, and P.B. Dusenbery, Geophys. Res. Lett.
17, 2313 (1990).
[69] G. Ganguli, Y.C. Lee, P.J. Palmadesso, Phys. Fluids, 31, 4 (1988)
[70] H.Romero, G. Ganguli, Y.C. Lee, P.J. Palmadesso, Phys. Fluids B 4, 7 (1992).
[71] G. Ganguli, Y.C. Lee, and P.J. Palmadesso, Phys. Fluids, 31, 10 (1988)
[72] H. Romero, G. Ganguli, and Y.C. Lee, Phys. Rev. Lett. 69, 24 (1992)
[73] G. Ganguli, M. J. Keskinen, H. Romero, R. Heelis, T. Moore, and C. Pollock, J.
Geophys, Res, 99, 8873 (1994).
[74] K.E. Atkinson, An Introduction to Numerical Analysis (Wiley, New York, 1978),
Sec. 2.10.
133