Whistler mode chorus is one of the most important emissions affecting the energization of the radiation belts. Recent laboratory experiments that inject energetic electron beams into a cold plasma have revealed several spectral features in the nonlinear evolution of these instabilities that have also been observed in high-time resolution in situ wave-form data. These features include (1) a sub-element structure which consists of an amplitude modulation on time-scales slower than the bounce time, (2) closely spaced discrete frequency hopping that results in a faster apparent frequency chirp rate, (3) fast frequency changes near the sub-element boundaries, and (4) harmonic generation. In this paper, we develop a finite dimensional self-consistent Hamiltonian model for the evolution of the resonant beam of electrons. We analyze a single wave case and demonstrate that the instability occurs due to a Krein collision, which manifests as a coupling between a negative and positive energy mode. This analysis revealed that the nonlinear evolution of the spectrally stable fixed-points of the self-consistent Hamiltonian develop a sub-packet structure similar to that of space observations. We then analyze the case of two whistler waves to show that the model reproduces the nonlinear harmonic generation and leads to a hypothesis for the closely spaced frequency hopping observed in laboratory experiments and space data.

Whistler mode chorus waves are some of the most spectacular electromagnetic waves in the earth's magnetosphere. They come in large amplitude bursts called elements, where the frequency is observed to change as a function of time, and their amplitude can reach up to δB/B0.01.1 They occur in the outer radiation belts (often in the dawn sector),2 and are believed to be driven by a temperature anisotropy created by the velocity dependent azimuthal drift of energetic electrons.3,4 They can also be triggered by VLF transmitters or lightning generated whistlers5 and are thus known as “triggered emissions.” Chorus is believed to play an important role in the dynamics of the Earth's outer radiation belts primarily through the resonance between the Doppler shifted wave frequency and the cyclotron harmonics. They can pitch angle scatter electrons into the loss cone6 where they are lost, but more importantly they can energize electrons to relativistic energies.7,8 A theory that could explain the saturated amplitude of whistler mode chorus as well as the complex frequency-vs-time dynamics is important for understanding how these acceleration and loss processes occur and ultimately for space weather applications.

Plasmas with larger temperatures in the direction perpendicular to the magnetic field rather than parallel, T>T||, are unstable to whistler mode waves,9 but the linear instability is insufficient to explain the features of chorus such as the time-dependent frequency. It has long been recognized that resonant particles play a definitive role in the frequency chirping and nonlinear amplification observed.5,10 Recently, analysis of high time-resolution electromagnetic waveform data from NASA's Van Allen Probe mission has revealed oscillations in the frequency and growth rate that are consistent with the bouncing of electrons trapped in the electromagnetic wave structure.11 These trapping oscillations are direct evidence of the resonant nature of instability. Another striking feature and clue to the physical processes responsible for the nonlinear evolution of the wave-dynamics of chorus is its sub-element structure, which is defined by a slow-time scale modulation of the amplitude of the chorus element.12,13 Recently it was also demonstrated that a discrete change in the frequency occurs across the sub-element boundary.11 

Laboratory experiments in the Naval Research Laboratory's (NRL) Space Physics Simulation Chamber (SPSC) have been performed in which an energetic electron beam is launched into a cold argon plasma.14 Similar experiments have been carried out at University of California at Los Angeles (UCLA) in the Large Area Plasma Device.15 In both cases, large amplitude electromagnetic whistlers with frequencies that change coherently with time were observed. The waves were consistent with resonant interaction with energetic electrons. In the case of the NRL experiments, the wave-electron interaction most likely occurred over a region with little to no magnetic field inhomogeneity (UCLA experiments have begun to look at the role of inhomogeneity). One of the striking features of the wave dynamics is the sub-element structure where the amplitude of the wave-packets has a slow time-scale modulation, and a frequency hop, i.e., a discrete change in frequency, coincident with the boundary of the sub-elements. In this paper, we keep these features in mind, which were observed in both space and laboratory nonlinear whistler waves, when analyzing the nonlinear model equations.

The nonlinear self-consistent interaction of a beam of particles with a single wave has been studied by many authors, beginning with O'Neil, Winfrey, and Malmberg16,17 who developed a model for the interaction of an electrostatic wave with an electron beam. This model was shown to be Hamiltonian, and a reduced model was developed from it called the rotating bar model.18 These kinds of single wave models have been extended to other kinds of waves, including Alfven waves19 and whistler waves.20,21 For the case of whistler waves, the application was to explain chorus. The difference between these past works and the present work is in the analysis of the model and in the inclusion of multiple wave modes. These previous single whistler-wave models, also did not use the Hamiltonian structure of these equations, nor emphasize the role of the conservation of both energy and momentum. Thus it was not previously recognized that the growth of the waves comes about by total momentum conservation, nor was a Hamiltonian spectral analysis possible which has revealed a new class of nonlinear solution. These types of models have also been extended to include the interaction between multiple wave modes22 and different wave modes.23 

In Section II, we construct a self-consistent Hamiltonian model for the nonlinear interaction of waves and energetic near resonant electrons. We construct a general Hamiltonian that includes the action angle variables for waves as well as canonical variables for the particles. In Section III, we analyze a restricted version of the model that includes only a single parallel propagating whistler wave. Previous single whistler-wave models did not use the Hamiltonian structure of these equations, nor emphasize the role of conservation of both energy and momentum. In Section III, A we find and analyze the motion around the fixed point of the Hamiltonian, we demonstrate that the whistler instability occurs via a Krein collision.24 In Section III B, we compare this analysis to the traditional linear plasma theory. Motivated by this analysis, we demonstrate that the nonlinear evolution of the linearly unstable but spectrally stable wave modes is qualitatively different than the previously studied solutions. These new solutions exhibit a modulation of the amplitude of the wave-packet that is consistent with the sub-element structure observed both in laboratory experiments and space data. In Section IV, we present the laboratory data that exhibits frequency hopping of chirping whistlers, and we analyze a version of the model that includes multiple parallel propagating whistlers. We demonstrate that nonlinear induced scattering25 is included in this model. We examine, numerically the evolution of a pair of modes and demonstrate that a discrete frequency hopping can occur. In Section V, we draw some conclusions based on the present work.

The model we wish to explore is described by a uniform magnetic field along the z direction and a background cold plasma with an infinite helical electron beam. We choose the periodic boundary conditions in all three directions, such that in reality we are considering the grid (in x and y) of electron beams of infinite extent in the z direction.

Before beginning, we note that all of the following can be derived from a variational principle, however, the constructive derivation followed here offers some insight into the approximations and assumptions being made.

Since we are interested in using canonical Hamiltonian methods for the electron dynamics, we work with the electrodynamic potentials for the wave equations. We fix the gauge where Φ, the electrostatic potential, is zero so that, B=×A and E=(1/c)A/t and Maxwell's equations may be combined to form

2A+(·A)+1c22At2=4πcJbg+4πcJb,
(1)

where we have broken up the current into a background component Jbg which we will treat in a linearized manner, and the beam component Jb which we will treat nonlinearly. Next we introduce

A=lAl(t)eikl·xiωlt+c.c.,
(2)

where c.c is the complex conjugate, Al(t) has a slow time dependence, and l is a mode index. Note that at this point, we could have included a slow spatial dependence as well, but we are restricting our attention in this study to homogeneous systems. Then we use the linear susceptibility tensor χs for species s to describe the response of background current to the field

Jbg=lωl24πsχs·Aleikl·xωlt+c.c.,
(3)

at this point, we are ignoring the nonlinear response of the background plasma, and thus we are ignoring important effects such as nonlinear induced scattering and 3-wave processes, etc., which are due to the background electrons and have been observed experimentally,14 but to keep this first model somewhat tractable, we ignore these effects here. With the linear susceptibility, the wave equation can be written as

lD(ωl,kl)·Aleikl·xiωlt+c.c.=4πcJb,
(4)

where D(ωl,kl) is the linear dispersion tensor which in the cold plasma limit may be written as

D(ω,k)=(ω2c2Sky2kz2iω2c2H+kxkykxkziω2c2H+kykxω2c2Skx2kz2kykzkzkxkykzω2c2Pkx2ky2),
(5)

where

S=1sωps2ω2Ωs2,H=sωps2Ωsω(ω2Ωs2)P=1sωps2ω2,
(6)

and where Ωs=qsB/(msc) is the signed particle gyrofrequency for species s and ωps2=4πnsqs2/ms is the plasma frequency. We note that there is an arbitrary sign convention used here which we have taken to be the conventional one. We then take the Fourier transform in space to arrive at

D(ωl,kl)·Al(t)=4πc1VVd3xeikl·x+iωltJb,
(7)

where V is the volume of the periodic box that we consider and where we note that in a finite periodic box, it is necessary to choose V such that the modes considered have wavelengths equal to integer multiples of the box size considered to isolate the term on the left hand side in the usual manner. We then find a zero eigenvalue of the dispersion matrix and with the corresponding complex unit magnitude eigenvector, el we write Al=Alel. Next, we expand the dispersion matrix around the normal mode frequency, and operate from the left with el. Then we note that Eq. (7) should be viewed as an operator equation18 where ω is an operator so that the left hand side of Eq. (7) becomes

el·D(ωl,kl)·Al(t)el·(D(ωl0,kl)+iωD(ωl0,kl)t)·Alel,
(8)

and using the fact that el·D(ωl0,kl)·el(t)=0, we arrive at the evolution equation for the mode amplitude

Ȧl=i4πc(ωD)1VVd3xeikl·x+iωltel·Jb,
(9)

which gives the change in the amplitude and phase of mode l as a function of time due to the nonlinear current of the beam. Next we introduce the wave action and phase of mode l by defining, Al=βJlexp(iθl), where θl and Jl are real variables, and β is a real constant for later convenience. Furthermore, we write the complex polarization vector as el=ρ̂lexp(iφl), where ρ̂l and φl are real, and this notation is meant to be thought of component-wise. With these definitions, the real vector potential can now be written as

A=lβ2Jlρ̂lcos(kl·xωltθl+φl).
(10)

We describe the dynamics of the electrons in the beam by using a canonical non-relativistic Hamiltonian formulation, in guiding center coordinates, where we have taken the constant background magnetic field to be described by A0=(By)ey, so that the single particle Hamiltonian may be written as

H=ΩePϕ+pz22m+2mΩePϕΩeB[ Aycos(ϕ)+Axsin(ϕ) ]+ΩeBpzAz+12mΩe2B2| A |2,
(11)

where Ωe>0 and the cartesian position variables may be expressed in the canonical guiding center coordinates as

x=PY2|Ωe|+2Pϕm|Ωe|cos(ϕ)y=Y2PϕmΩesin(ϕ).
(12)

We note that in this work, we focus on the generation and saturation of chorus-like emissions, such that the relativistic equations are not necessary. Laboratory experiments of chorus emissions are successful with 10's of keV electrons. Next we use our definition of the vector potential (Eq. (10)) to write this as

H=ΩePϕ+pz22mlβ2mΩePϕJlΩeB[ ρlycos(Ψl+φlyϕ)+ρlycos(Ψl+φly+ϕ)ρlxsin(Ψl+φlxϕ)+ρlxsin(Ψl+φly+ϕ) ]2ΩeBpzlJlρ̂lzcos(Ψl+φlz)+mΩe2B2l,lβ2JlJl[ ρ̂lxρ̂lx(cos(ΨlΨl+φlxφlx)+cos(Ψl+Ψl+φlx+φlx))+ρ̂lyρ̂ly(cos(ΨlΨl+φlyφly)+cos(Ψl+Ψl+φly+φly))+ρ̂lzρ̂lz(cos(ΨlΨl+φlzφlz)+cos(Ψl+Ψl+φlz+φlz)) ],
(13)

where

Ψl=k·xωltθl.
(14)

Next, we need the expression for the cartesian velocities (vx,vy,vz) in terms of the new canonical variables

vx=1m(2m|Ωe|Pϕsin(ϕ)β2|e|clJlρ̂xcos(Ψl+φlx))vy=1m(2m|ΩePϕ|cos(ϕ)β2|e|clJlρ̂ycos(Ψl+φly))vz=1m(pzβ2|e|clJlρ̂zcos(Ψl+φlz)).
(15)

In addition to the linear current of the background cold plasma, we consider the current due to beam particles. Which, we can write as

Jb=enbVMi=0Md3pd3qv(q,p,t)δ(xr(q,p,t))δ(qqi(t))δ(ppi(t)),
(16)

where enbV/M is the beam charge density per particle inside the volume V, and where (qi(t),pi(t)) is the position in phase space of the ith beam particle. There are 9 delta functions and 6 integrals, which leave three delta functions. Here we write

Jb=enbVMi=0Mv(ϕi(t),Pϕi(t),z)×δ(x+PYi(t)2|Ωe|2Pϕi(t)m|Ωe|cos(ϕi(t)))×δ(yYi(t)+2Pϕi(t)mΩesin(ϕi(t)))×δ(zzi(t)).
(17)

Using the expression for the beam current density in the equation giving the evolution of the complex mode amplitude (Eq. (9)), we arrive at

Ȧl=i4πenbc(ωD)1Mi=0MeiΨliθlel·v,
(18)

where

el·v=(ρ̂leiφl)T·(2|Ωe|Pϕmsin(ϕ)2|Ωe|Pϕmcos(ϕ)pz)+ΩeB(ρ̂leiφl)T·A,
(19)

putting in our representation for Al and equating real and imaginary parts of Eq. (18), and after some straightforward algebra we can write the evolution equations for Jl and θl as

J̇l=4πenbcβ(ωD)1Mi=0M{ 2ΩePϕJlmρ̂lx[ cos(ψl+φlxϕ)cos(ψl+φlx+ϕ) ]+2ΩePϕJlmρ̂ly[ sin(ψl+φlyϕ)+sin(ψl+φly+ϕ) ]+2pzJlρ̂lzsin(ψl+φlz)2βΩeBlJlJl[ ρ̂lxρ̂lxsin(ψl+φlxψlφlx)+ρ̂lxρ̂lxsin(ψl+φlx+ψl+φlx)+ρ̂lyρ̂lysin(ψl+φlyψlφly)+ρ̂lyρ̂lysin(ψl+φly+ψl+φly)+ρ̂lzρ̂lzsin(ψl+φlzψlφlz)+ρ̂lxρ̂lxsin(ψl+φlz+ψl+φlz) ] }
(20)

and

θ̇l=4πenbcβ(ωD)1Mi=0M{ 122ΩePϕmJlρ̂lx[ sin(ψl+φlxϕ)+sin(ψl+φlx+ϕ) ]+122ΩePϕmJlρ̂ly[ cos(ψl+φlyϕ)+cos(ψl+φly+ϕ) ]+pzρ̂lzJlcos(ψl+φlz)βΩeBlJlJl[ ρ̂lxρ̂lxcos(ψl+φlxψlφlx)+ρ̂lxρ̂lxcos(ψl+φlx+ψl+φlx)+ρ̂lyρ̂lycos(ψl+φlyψlφly)+ρ̂lyρ̂lycos(ψl+φly+ψl+φly)+ρ̂lzρ̂lzcos(ψl+φlzψlφlz)+ρ̂lxρ̂lxcos(ψl+φlz+ψl+φlz) ] }.
(21)

We note that the wave equations come in a form that can be seen to come from the single particle Hamiltonian, except for the constant factor outfront. Now we can see that the purpose of the constant factor β was to make the variables, θl,Jl, canonically conjugate variables. If we take the M-particle, N-wave Hamiltonian to be the sum of the single particle Hamiltonians indexed by the particle label then we find

H=i=0M{ ΩePϕi+pzi22mlβl2mΩePϕiJlΩeB[ ρlycos(Ψli+φlyϕ)+ρlycos(Ψl+φly+ϕi)ρlxsin(Ψli+φlxϕi)+ρlxsin(Ψli+φlx+ϕi) ]2ΩeBpzilJlρ̂lzcos(Ψli+φlz)+mΩe2B2l,lβlβlJlJl[ ρ̂lxρ̂lx(cos(ΨliΨli+φlxφlx)+cos(Ψli+Ψli+φlx+φlx))+ρ̂lyρ̂ly(cos(ΨliΨli+φlyφly)+cos(Ψli+Ψli+φly+φly))+ρ̂lzρ̂lz(cos(ΨliΨl+φlzφlz)+cos(Ψli+Ψli+φlz+φlz)) ] },
(22)

where we note that ψli is not a canonical coordinate of phase space, but is rather a definition of a conglomeration of variables. Then we see that the equations, θ̇l=H/Jl and J̇l=H/θl will give the appropriate wave dynamics if we choose

βl=(4πnb(ωDl)1M)1/2.
(23)

Analyzing this Hamiltonian in general may be of limited use, however, there are reasons that we find it useful, namely: (1) as a starting point for analyzing simpler cases as we will do in Sections III and IV, (2) as a means to write a general purpose code to analyze beam-wave interactions. It should also be pointed out here that the consideration of oblique whistler wave modes is only more difficult because the single wave case does not have a 2-dimensional phase space such that the electron dynamics may be easily visualized. Once multiple waves of any kind are included, the phase space becomes inherently more difficult to visualize and analyze.

Here, we consider the case of a single parallel propagating whistler in a background cold plasma interacting with the electron beam. This case has been analyzed before,20,21 but we feel that in light of new data, a second study of these equations to highlight certain features may be of use. In this case, we use the dispersion tensor in Eq. (5) with kx=ky=0. The determinant of the dispersion tensor gives the dispersion relation

(ω2c2Sk2)2ω4c4H2=0,
(24)

which can be solved for k2 as

k2=ω2c2(S±H),
(25)

where the + sign gives the right hand circularly polarized whistler dispersion relation, whose polarization vector is given by

êl=12(i10),ρ̂l=(12120),φ̂l=(π200),
(26)

so that

D(ω,k)êl·D(ω,k)·el̂=ω2c2(S+H)+k2,
(27)

where S + H = R and

R=1sωps2ω(ω+Ωs).
(28)

With two species and in the limit ωΩi and with me/mi1, the dispersion relation can be written as

Ωe2c2D(ω¯,k¯)ω¯2(1+δ2ω¯(1ω¯))2k¯2,
(29)

where ω¯=ω/Ωe,k¯=ck/Ωe, and δ2=ωpe2/Ωe2. The self-consistent Hamiltonian reduces significantly to

H=i=0M{ ΩePϕi+pzi22m2βmΩePϕiJΩeBcos(kziωtθϕi)+mΩe2B2β2J }.
(30)

There are two things about this equation to note. First, the coupling between the fields and the particles occurs only in one term. Second, the arguments of the trigonometric function can be simplified by introducing a single new particle variable, ψi=kziωktϕi. To take advantage of this, we introduce the type-2 mixed-variable generating function that takes us from (t;zi,ϕi;Pϕi,pzi) to (t;Ψi,ξi;PΨi,Ii) via

F=i=0M(kziωtϕi)PΨi+ziIi,
(31)

so that the transformation equations are

ψi=PψiF2=(kziωtϕi)ξi=IiF2=zipzi=ziF2=kPψi+IiPϕi=ϕiF2=Pψi,
(32)

note that the generating function is time-dependent, thus the new Hamiltonian becomes H=H+tF2. The new Hamiltonian becomes

H=i=1M{ (Ωeω)PΨi+12m(kPψi+Ii)22βmΩePΨiJΩeBcos(Ψiθ)+mΩe2B2β2J }.
(33)

This Hamiltonian has M + 1 degrees of freedom, Ii are constants of the motion. We can make this dimensionless

h=i=1M{ (1ω¯)pψi+12(k¯pψi+Ii)22αpψiJMcos(ψiθ)+α2JM },
(34)

where time is measured in Ωe, energy with mc2, particle momentum Pψi and wave action with mc2/Ωe, and the new constant α=δ2(nb/n0)/(Ωe/c2ω¯D). This Hamiltonian has two immediately apparent constants of motion which can be seen from the equations of motion, the total energy h, and the total momentum

P=J+ipψi.
(35)

From momentum conservation, we see that as the beam particles lose pψi, the wave amplitude grows, by adjusting J, thus momentum conservation controls the growth rate.

In Figure 1, we show the evolution of the normalized wave amplitude of a parallel propagating whistler driven by an electron beam of energy 40 keV, 45° pitch angle, and a density that was 1% of the cold plasma density. The resonant frequency was about ω=0.225Ωe, and the parallel wave vector was, kc/Ωe2.769. The background plasma parameters were chosen typical of chorus so that δ5.1. We used 20 time-steps per cyclotron period and the explicit Runge-Kutta method of order (4)5 due to Dormand and Prince to achieve an energy and momentum conservation near the machine precision (1013%). We initiated the particles around the fixed point (as discussed in Section III A) with ψi0=2πi/M where M = 1000 was the number of particles, and added a perturbation of the fixed point from the unstable eigenvector found from the linear analysis. We see that the mode initially grows at the linear rate (plotted in green in the figure with an initial offset so that it can be distinguished), and the mode saturates at a finite amplitude.

FIG. 1.

Evolution of normalized wave action for parallel propagating whistler driven by an electron beam of energy 40 keV, 45° pitch angle, and a density that was 1% of the cold plasma density.

FIG. 1.

Evolution of normalized wave action for parallel propagating whistler driven by an electron beam of energy 40 keV, 45° pitch angle, and a density that was 1% of the cold plasma density.

Close modal

In Figure 2, we show the evolution of the resonant particles at four different snapshots in time. In these figures, we plot the phase space position of each particle as a blue dot, and we draw contours of the single particle Hamiltonian using the wave amplitude and phase at that given instance. At t=90.0Ωe, the resonant island in phase space has opened up, and the particles are clearly moving to smaller momentum values, and thus by conservation of total momentum, the wave is growing. At t=135.0Ωe, the instability is near the saturation, and we see that a good fraction of particles have been trapped in the island structure and are again gaining momentum as they orbit around the island center near pψ=0.04. Conservation of momentum providing the mechanism of saturating the instability. In the third panel, at t=180.0Ωe, the particles are oscillating around the fixed points while low-density tendrils of particles are spread throughout the phase space. At this point, one can see that the rotating bar model18 may be a fair approximation for times just after the nonlinear saturation. However, in the final panel, at t=2500.0Ωe many of the particles are exhibiting apparently chaotic orbits, with two separate coherent structures observed. The rotating bar has broken up into two clumps by this time. In Figure 3, we show the evolution of the instantaneous frequency (θ̇ and the instantaneous growth rate J̇/J. We find that in this case, there is an overall shift in the frequency due to the beam particles, but there is no secular change in the frequency that would correspond to either laboratory or space observations of the chorus phenomenon.

FIG. 2.

Evolution of resonant particles in phase space at four different times corresponding to the wave evolution depicted in Figure 1.

FIG. 2.

Evolution of resonant particles in phase space at four different times corresponding to the wave evolution depicted in Figure 1.

Close modal
FIG. 3.

Evolution of instantaneous frequency (top) and growth rate (bottom) as a function of time.

FIG. 3.

Evolution of instantaneous frequency (top) and growth rate (bottom) as a function of time.

Close modal

To perform an analysis of the fixed points, we can avoid singularities caused by zero wave-amplitude by making a canonical transformation of the Hamiltonian of Eq. (34) using X=2Jcos(θ),Y=2Jsin(θ), and X is a momentum variable and Y a position variable so that the Hamiltonian becomes

h=i=0M{ (1ω¯)pψi+12(k¯pψi+Ii)22αpψi2MXcos(ψi)2αpψi2MYsin(ψi)+α212N(X2+Y2) }.
(36)

Then we look for a fixed point of the Hamiltonian with zero wave amplitude, which leads to a condition on pψi, from ψ̇i=pψih=0, such that

0=1ω¯+k¯(k¯pψ+Ii)=Ωω+kvz,
(37)

which means that if all of the particles are in exact Doppler shifted cyclotron resonance with the wave and the wave amplitude is zero X=Y=0, then we are at a fixed point of the Hamiltonian for any values of ψi that are evenly distributed (as is typically done). We then form a 2M + 2 dimensional vector z=(δψi,δpψi,δY,δY)T and calculate the linearized motion around the fixed point, which may be written in the matrix form as z.=Lż, where L is a block matrix of the form

L=(0k¯21[ α12Mpψisin(ψi) ][ α12Mpψicos(ψi) ]00[ 2αpψi2Mcos(ψi) ][ 2αpψi2Msin(ψi) ](2pψi2Mαsin(ψi))(αcos(ψi)2Mpψi)0α2(2pψi2Mαcos(ψi))(αsin(ψi)2Mpψi)α20),
(38)

where 0 is the M × M zero matrix, 1 is the M × M identity matrix, and brackets, e.g., [xi] refers to M × 1 column vectors and parenthesis, e.g., (xi), to 1 × M row vectors. Since this matrix comes from a Hamiltonian system, its eigenvalues σj have certain properties, namely, that the eigenvalues come in real pairs ±σj, imaginary pairs ±iσj, zero, or complex quadruplets ±a±ib. The solution to the linearized motion can be written as a sum jcjξjexp(σjt), where ξj is the associated eigenvectors. For a Hamiltonian that depends smoothly on some parameter, as we have here, there are only two ways for a system to go unstable, (1) a saddle node bifurcation where a pair of imaginary eigenvalues merge at zero and then split onto the real axis, or (2) a Krein bifurcation where a pair of eigenvalues (of opposite signature, as explained below) collide and split into quadruplets in the complex plane.

In Figure 4, we plot the eigenvalues of the matrix given by Eq. (38) for a beam density of 1% of the background, a beam energy of 11 keV, and background plasma parameters corresponding to δ=5.11, and for two different values of the pitch angle: 17.0° which has no unstable solutions (all eigenvalues are on the imaginary axis), and 18.0° which has both stable and unstable solutions. A Krein collision occurs between these parameter values. The motion of the eigenvalues as we vary the pitch is depicted by the arrows. The eigenvalues are plotted using a color corresponding to the sign of the signature sj which was computed by

sj=sign(iω(ξj,ξj*)),
(39)

where ω is the symplectic 2-form. Red corresponds to a positive energy and Blue to a negative energy. We find that the instability here arises through Krein collisions, which in this case is the coupling of a positive and negative energy wave. Recently, this mechanism has been extended to the complex Hamiltonian structure of the two-stream instability.26 Such an extension would be necessary to consider the case of a warm electron beam.

FIG. 4.

Eigenvalues of linearized motion around zero amplitude fixed point for parallel whistler case. Positive signature (energy) corresponds to the red dots, and negative signature corresponds to the blue dots. The motion as we increase the pitch angle from 17.0° to 18.0° is depicted by the arrows.

FIG. 4.

Eigenvalues of linearized motion around zero amplitude fixed point for parallel whistler case. Positive signature (energy) corresponds to the red dots, and negative signature corresponds to the blue dots. The motion as we increase the pitch angle from 17.0° to 18.0° is depicted by the arrows.

Close modal

Spectral stability does not imply that the system under investigation is linearly unstable. Indeed, the traditional plasma physics literature has already worked out the dispersion relation for an electron beam interacting with a parallel propagating whistler.27 The dispersion relation is

ω¯2(1+δ2ω¯(1ω¯))=nbnδ2(ω¯k¯V||/c(ω¯k¯V||/c1)+V2c(ω¯2k¯2)(ω¯k¯V||/c1)2),
(40)

where ω¯,k¯, and δ are defined before. V||,V are the parallel and perpendicular beam velocity.

In Figure 5, we plot the linear growth rate Im(ω¯) as computed from Eq. (40) as a function of k¯ for two cases. Both cases consider a 50 keV beam propagating at a pitch angle of 17°, through a plasma with δ=5.115. The first case, with a beam density of nb/n0=0.015, corresponding to a spectrally stable case is shown in the top panel. The second case, with a beam density of nb/n0=0.01, corresponds to a spectrally unstable case (after a Krein collision). The resonant k¯=2.04 is indicated by a vertical grey line. As dictated by the nonlinear model and the requirements of the fixed point, this k¯ and corresponding ω¯ are determined by simultaneously satisfying the background dispersion relation without the beam and the wave-particle resonance condition. Thus, we see that spectrally stable solutions are solutions where the ω and k determined by the resonance condition and the unmodified dispersion relation are linearly stable. This can only happen when the beam density is large enough to modify the real frequency enough.

FIG. 5.

Linear growth rate as a function of k¯ for case that is spectrally stable (top) and unstable (bottom), where the resonant k is indicated by the vertical grey line.

FIG. 5.

Linear growth rate as a function of k¯ for case that is spectrally stable (top) and unstable (bottom), where the resonant k is indicated by the vertical grey line.

Close modal

Next, we investigate the nonlinear evolution of a spectrally stable but linearly unstable solution. To do this, we choose a 11 keV beam propagating at a pitch angle of 17.975° in a plasma with δ=5.115 with a beam density of nb/n0=0.01. We choose a ω¯=0.31407694 and k¯=3.47541495 such that the linear dispersion relation (without the beam) is satisfied as well as the Doppler shifted resonance condition. We compute the eigenvalues of the matrix L of Eq. (38) and confirm that all the eigenvalues lie on the real axis and is thus spectrally stable. We then perturb the beam particles by 4%. This corresponds to shifting the beam energy and pitch angle slightly so that ω¯ and k¯ are not quite in resonance with the beam but close. With the perturbed beam energy and pitch angle, we can calculate the linear growth rate according to Eq. (40) and determine the linear growth rate to be Im(ω¯) = 0.00271396. In Figure 6, we plot the evolution of the wave action for this case on a log scale. We find that the wave proceeds to grow according to the linear theory. In green, we plot the evolution that we would expect based on the linear growth rate. At some amplitude, the wave becomes damped. In red, we plot the negative growth rate, which from Figure 5, we see is also a solution to the linear dispersion relation. The mode goes through this cycle three more times before the behavior changes qualitatively. We note that these parameters were chosen to make this plot have these specific features that we find interesting. If we choose too large of a growth rate, the transition to a qualitatively different behavior happens earlier in time without the apparently reversible behavior. Whereas too small of a growth rate and the reversible behavior happens for too long a time. In this simulation, the total error in energy conservation is 10−4% and similarly for momentum conservation.

FIG. 6.

Evolution of wave-action for non-resonant mode that is spectrally stable, but linearly unstable.

FIG. 6.

Evolution of wave-action for non-resonant mode that is spectrally stable, but linearly unstable.

Close modal

In Figure 7, we show what a time series of the wave amplitude looks like near the transition to the irreversible behavior. We find that the amplitude is modulated on a slow time scale. This modulation does not happen for the case that we studied that was spectrally unstable and shown in Figure 1, however, this modulation is reminiscent of the sub-element structure that has been observed in both laboratory experiments14 and in space observations.11–13 In Figure 8, we plot a time series of one component of the magnetic field during a chorus emission. This was recorded by Van Allen probe A on October 14, 2012 during a geomagnetic storm. The boundaries of the sub-elements are indicated with vertical grey lines. This model result suggests a new hypothesis to be explored, that the modulation of the amplitude is due to slightly non-resonant particles. This may arise due to the wave parameters slipping away from the resonant particles as it propagates in an inhomogeneous magnetic field.

FIG. 7.

Reproduction of time series as would be observed in laboratory indicating the modulation of the amplitude on a slow time-scale.

FIG. 7.

Reproduction of time series as would be observed in laboratory indicating the modulation of the amplitude on a slow time-scale.

Close modal
FIG. 8.

One component of the magnetic field of a chorus emission from high-time resolution waveform data obtained by the EMFISIS instrument on board the Van Allen Probe A. Sub-element boundaries are indicated by grey line.

FIG. 8.

One component of the magnetic field of a chorus emission from high-time resolution waveform data obtained by the EMFISIS instrument on board the Van Allen Probe A. Sub-element boundaries are indicated by grey line.

Close modal

In Figures 9 and 10, we show the evolution of the instantaneous frequency (θ̇ and the instantaneous growth rate J̇/J over different time intervals. Figure 9 shows the time period corresponding the reversible growth and decay. We see that the period of change in frequency is twice that of the period of the linear growth rate. We also see that there are large oscillations in the frequency and growth rate near where the amplitude of the wave is smallest. Figure 10 shows the time period where the more complicated motion is taking place. There are very large amplitude changes in the frequency near the points where the wave-amplitude becomes small. This large change in the frequencies is also seen near the sub-element boundaries in the space data of chorus.11 

FIG. 9.

(Top) Change in real frequency and (bottom) instantaneous growth rate for case where mode is spectrally stable but linearly unstable.

FIG. 9.

(Top) Change in real frequency and (bottom) instantaneous growth rate for case where mode is spectrally stable but linearly unstable.

Close modal
FIG. 10.

(Top) Change in real frequency and (bottom) instantaneous growth rate for case where mode is spectrally stable but linearly unstable. Plots show the time period near transition to irreversible behavior.

FIG. 10.

(Top) Change in real frequency and (bottom) instantaneous growth rate for case where mode is spectrally stable but linearly unstable. Plots show the time period near transition to irreversible behavior.

Close modal

In Figure 11, we show the phase space, similar to Figure 2, at four different times. What is missed in this plot (but can be seen in movies) is that the island rapidly rotates around ψ, due to the beam being slightly off resonance. At t=14500Ωe, we can see that there are two portions of the beam that are being distorted in a circular pattern, unlike the spectrally unstable case where there was only one rotating bar like structure. By t=15000Ωe, these bar like structures have detached from each other, and by t=25000Ωe, the particles appear as two clumps, one near the island, and one near the x-point, with some particles being drawn into chaotic orbits. These figures suggest that a two-rotating bar model may explain the sub-packet structure.

FIG. 11.

Particle phase space at four different times for the case that is spectrally stable, but linearly unstable.

FIG. 11.

Particle phase space at four different times for the case that is spectrally stable, but linearly unstable.

Close modal

Ongoing experiments, investigating the details of triggered emissions are being conducted at the NRL in the SPSC.14 The experiments are conducted in the afterglow plasma from a pulsed helicon source into an expanding magnetic field geometry going from uniform 420 G at the source to uniform 20 G field in the experimental region. An energetic electron beam is used propagating towards the source region so that the interaction between the whistler and the electron beam is believed to be occurring in a region of homogeneous magnetic field. Relevant parameters include: electron beam energy Eb=1.8 keV, beam pitch angle 40°, ratio of beam to background plasma density nb/n0=0.050.13, and ratio of plasma to cyclotron frequency δ=1.759.2, where the ranges are due to the decaying background plasma.

The beam-plasma instability generates a broadband spectrum of whistler waves, and as the ratio of beam to background density increases, the amplitude of the broadband whistlers increases. When the ratio is large enough, a subset of the generated whistler waves can interact with the beam and triggers a large amplitude chirping whistler wave, a portion of which is shown in Figure 12. Figure 12 highlights that the chirping wave is composed of slowly chirping elements that are generated in frequency hops, where the ensemble looks like a fast chirping whistler wave. The dashed lines illustrate the linear chirp rate that is maintained over the many frequency hops.

FIG. 12.

Frequency vs. time of laboratory experiment revealing discrete hops of frequency.

FIG. 12.

Frequency vs. time of laboratory experiment revealing discrete hops of frequency.

Close modal

Motivated by the laboratory observation of multiple whistlers, we consider the case of multiple interacting parallel propagating whistler waves in a background cold plasma with a resonant electron beam. The polarization vector and dispersion relation are the same as in Section III, the only difference being that we retain the sum over modes. The self-consistent Hamiltonian of Eq. (22) reduces to

H=i=0M{ ΩePϕi+pzi22ml2βlmΩePϕiJlΩeBcos(klziωltθlϕi)+mΩe2B2l,lβlβlJlJlcos((klkl)zi(ωlωl)t(θlθl)) }.
(41)

Note that the presence of at least two wavelengths with a gyrophase dependence ϕi prevents transforming into the moving frame and reducing the problem to 1 degree of freedom per particle without further approximation. For now, we simplify the analytical manipulations by considering only two waves and use a type two generating function the same as the single wave case, Eq. (31), except that the k and ω in the generating function must be indexed to a particular wave, l = 1. The Hamiltonian is still time-dependent so we must transform the Hamiltonian once more, using another type-2 generating function of the form F2=(ω2ω1)τ+θ2)I2, which removes the time dependence from the arguments of the cosines. Then we introduce the dimensionless variables similar to what was done for the single wave case to arrive at the dimensionless Hamiltonian

h=i=0M{ (1ω¯)pψi+12(k¯pψ+pξ)2+α12MJ12+α22+(ω2ω1)MJ222α1pψiJ1Mcos(ψiθ1)2α2pψiJ2Mcos(Δkξiθ2+Ψi)+α1α2J1J2Mcos(Δkξi(θ2θ1)) },
(42)

where Δk=k2k1. In addition to h, examination of the equations of motion reveals that there are two additional constants of motion

P1=J1+J2+ipψP2=ΔkJ2+ipξ,
(43)

representing the conservation of two components of the total momentum. It is interesting to point out, that by considering the two waves here, the interaction term (the last term in the Hamiltonian) becomes secular at the nonlinear Landau beat-wave resonance condition when α22Δω

ΔωΔk=k1pψ+pξv||,
(44)

where Δω=ω2ω1. If we considered a beam of particles at zero temperature, this Hamiltonian would allow one to study the beat-wave trapping dynamics and the saturation of nonlinear induced-scattering. This will be addressed in a future article.

In space observations of chorus, chirping harmonic emissions are often observed. In Figure 13, we show the short-time Fourier transform (using the FFT with a square window) of high-time resolution wave data from the EMFISIS instrument on board the Van Allen Probe A satellite. The figure shows a main chorus element that starts at 0.12 s at a frequency of about 600 Hz. Near 0.15 s a peak near 1200 Hz appears. The magnitude of the FFT coefficient corresponding to the Harmonic is about 15% of the magnitude of the coefficient of the primary wave.

FIG. 13.

Short-time Fourier transform of one component of the magnetic field from the EMFISIS instrument aboard NASA's Van Allen Probe. Harmonic structure can be seen.

FIG. 13.

Short-time Fourier transform of one component of the magnetic field from the EMFISIS instrument aboard NASA's Van Allen Probe. Harmonic structure can be seen.

Close modal

To model the harmonic generation seen in space observations, we choose the background plasma to have δ=5.115, and resonant electrons at 65 keV with a pitch angle of 25°, and with a beam density to background density nb/n0=1×104. The unstable wave that resonates with this beam has a frequency of ω¯1=0.123 and a wavelength of k¯1=1.92. We consider a second wave with k¯2=2k¯1 and a frequency ω¯=0.358 that satisfies the dispersion relation. With these parameters, the dimensionless constants α1=8.74×103 and α2=6.38×103. In this case, the simulation domain only needs to be as long as a single wavelength of the longest wavelength mode. In Figure 14, we plot the wave action coordinates J1 (blue) and J2 (green) as a function of time. Energy and both total momenta were conserved to be within 10−6%. The resonant mode grows linearly from the start in agreement with the spectral theory. The harmonic does not grow until the first wave has reached sufficient amplitude (J103), at which point it grows rapidly and then saturates at a low amplitude. In the final nonlinear state, the harmonic has a wave amplitude that is about 10% of the primary wave, consistent with the space observations.

FIG. 14.

J vs. time of primary unstable wave (blue) and harmonic (green).

FIG. 14.

J vs. time of primary unstable wave (blue) and harmonic (green).

Close modal

We next consider the nonlinear evolution of a resonant wave mode (one whose unperturbed dispersion relation is such that the Doppler shifted resonance condition is satisfied) along with a non-resonant wave mode. For this case, we choose δ=5.115, and resonant electrons at 11 keV with a pitch angle of 20.5°, and with a beam density of nb/n0=0.01. The resonant mode has a frequency of ω¯=0.318 and a wave-vector of k¯=3.508, while the non-resonant mode has a frequency of ω¯=0.402 and a wave-vector of k¯=4.21. Both modes are linearly unstable. In Figure 15, we plot the evolution of the wave-action variables J1 and J2. The non-resonant mode in this case has a larger linear growth rate as can be seen with the solid-green curve than the resonant mode (with solid-blue). We confirm that both modes are linearly unstable by individually turning each mode off (by setting the corresponding α to zero). These curves are plotted with a dotted line of the same color. The figure shows that when either mode is allowed to evolve by itself, the amplitude the mode reaches is about the same. However, when the modes are allowed to grow simultaneously, we see that the fastest growing mode reaches the saturated state first, but does not prevent the evolution of the resonant mode from proceeding. As the resonant mode grows, the non-resonant mode loses wave-action and eventually becomes subdominant.

FIG. 15.

Evolution of wave action variables. Blue curves represent the evolution of the resonant mode. Green curves represent the evolution of the non-resonant mode. Dotted curves correspond to the evolution of that mode in the absence of the other.

FIG. 15.

Evolution of wave action variables. Blue curves represent the evolution of the resonant mode. Green curves represent the evolution of the non-resonant mode. Dotted curves correspond to the evolution of that mode in the absence of the other.

Close modal

In Figure 16, we calculate the short-time FFT of the simulated time-series of the evolution of the two modes. In this figure, the evolution appears as a discrete change in the frequency. In the laboratory data presented in Figure 12, a discrete downshift in frequency is observed. Those modes are believed to be oblique whistlers and thus the quantitative comparison of the simulation results cannot be done at this time. What is interesting about these results, is that the resonant mode is the dominant nonlinear mode, and not the mode with the largest growth rate, as is often assumed. These calculations lead to a new hypothesis, to be investigated in the future, that the discrete changes in frequency are due to the natural evolution of non-resonant modes with larger growth rates in a quantized system where only a discrete set of k is allowed.

FIG. 16.

Short-time FFT of simulated data of the evolution of two closely spaced modes.

FIG. 16.

Short-time FFT of simulated data of the evolution of two closely spaced modes.

Close modal

In this paper, we have constructed a general model that captures the nonlinear self-consistent dynamics of energetic electrons in a finite number of waves supported by a background plasma. We have analyzed two classes of models, a single wave model consisting of a parallel propagating whistler wave interacting with an energetic electron beam, and a model with two parallel propagating whistlers. Models like these offer a number of advantages. They are simple to implement numerically, consisting only of a series of ordinary differential equations. They are efficient to run. Simulations can be run without a lot of computing resources, making the dynamics easier to diagnose and understand, and enabling long-time evolution. In addition, the compact nature of the equations will, in the future, allow efficient analytical perturbation methods, such as Lie transform perturbation theory, to be applied to explain the observed numerical dynamics.

There are a number of ways to extend these types of models and yet retain the advantages. One would be to include the effect of inhomogeneity (as was done for electrostatic waves with DC electric fields22). Another way to extend these types of models, is to include nonlinearities of the background plasma, either perturbatively or not. This may offer an efficient way to investigate the effects of trapped particles on weak-turbulence processes such as the saturation of nonlinear Landau damping by particles trapped in the beat-potential28 or the parametric decay of a whistler wave into another whistler wave and an electron acoustic wave that then forms a “time-domain structure.”29 In space it is not necessarily the case that the electron distribution has such sharp beam-like features as has been modeled here. Thus, for space applications, it is important to investigate whether a more broad electron distribution can give rise to similar phase space structures that can be responsible for the observed amplitude modulations in space. One of the real strengths of this kind of model is the “taming” of the singular nature of resonant particle interactions.

In this work, we were motivated to explain recent laboratory and space observations of the nonlinear evolution of whistler waves. We used a simple model of a beam of electrons interacting with a single parallel propagating whistler wave to investigate the trapping motion. We demonstrated that the instability arises due to a Krein collision between a negative and positive energy mode. This analysis turned out to be more than academic. The linear (as opposed to spectral) theory revealed that modes are unstable in a previously unexplored parameter regime. The nonlinear evolution of these slightly off-resonant modes was qualitatively different than the evolution of the typical on-resonant evolution. The simulations revealed a competition between a growing and damping solution that eventually led to a nonlinear state that carried a slow amplitude modulation on a time-scale longer than the bounce-time, with rapid and large changes in the frequency near the minima of the amplitude modulations. This behavior is similar to observations of chirping whistlers in real plasmas. And yet, it was found that parallel propagating whistlers do not have a secular change in the frequency as is typically observed in whistlers in real plasmas. It remains to be seen, what the physical mechanism is that causes the observed frequency chirp of chorus. The most conventional hypothesis is that the inhomogeneous magnetic field is responsible for the chirp, and the frequency adjusts to maintain resonance, but the mechanism for how this happens has not been explained.

Another feature that is clearly seen in laboratory experiments is a discrete frequency hop, i.e., a discontinuous jump in the frequency. Similar features have been observed in the late-time evolution of chorus observations in space.11 We have attempted to model this as the interaction between two closely spaced frequencies that exchange energy at a finite amplitude. This hypothesis raises other interesting theoretical questions. Can the nonlinear generation of secondary waves (closely spaced or not) be responsible for saturating the instability, by either removing wave energy from the resonant source of instability drive, or by causing the electrons to become detrapped and disrupting the nonlinear current that drives the nonlinear instability? These are fascinating questions that are not only important for radiation belt dynamics. The same physics, with the same implications for energetic particle transport, of the trapping of energetic particles and chirping waves is important for tokamak plasmas as well,30,31 and speaks to the universal behavior of nonlinear wave-particle interactions.

We acknowledge the hard work of the entire Van Allen Probes team on the instrument design, testing, calibration, operations, and open dissemination of science data. We would like to explicitly thank Craig Kletzing and George Hospodarsky for insightful discussions about data from the EMFISIS instrument. This work was partially supported by the Naval Research Laboratory base program and by the National Aeronautics and Space Administration under Grant No. NNH15AZ90I.

1.
L. B.
Wilson
,
C. A.
Cattell
,
P. J.
Kellogg
,
J. R.
Wygant
,
K.
Goetz
,
A.
Breneman
, and
K.
Kersten
,
Geophys. Res. Lett.
38
,
L17107
(
2011
).
2.
W.
Li
,
R. M.
Thorne
,
V.
Angelopoulos
,
J.
Bortnik
,
C. M.
Cully
,
B.
Ni
,
O.
LeContel
,
A.
Roux
,
U.
Auster
, and
W.
Magnes
,
Geophys. Res. Lett.
36
,
L09104
, doi: (
2009
).
3.
A. B.
Collier
and
A. R. W.
Hughes
,
Adv. Space Res.
34
,
1819
(
2004
).
4.
K.
Min
,
J.
Lee
, and
K.
Keika
,
J. Geophys. Res.
115
,
A00I02
(
2010
).
5.
Y.
Omura
,
D.
Nunn
,
H.
Matsumoto
, and
M. J.
Rycroft
,
J. Atmos. Terr. Phys.
53
,
351
(
1991
).
6.
B. T.
Tsurutani
,
O. P.
Verkhoglyadova
,
G. S.
Lakhina
, and
S.
Yagitani
,
J. Geophys. Res.
114
,
A03207
(
2009
).
7.
J. C.
Foster
,
P. J.
Erickson
,
D. N.
Baker
,
S. G.
Claudepierre
,
C. A.
Kletzing
,
W.
Kurth
,
G. D.
Reeves
,
S. A.
Thaller
,
H. E.
Spence
,
Y. Y.
Shprits
, and
J. R.
Wygant
,
Geophys. Res. Lett.
41
,
20
, doi: (
2014
).
8.
D. L.
Turner
,
V.
Angelopoulos
,
W.
Li
,
M. D.
Hartinger
,
M.
Usanova
,
I. R.
Mann
,
J.
Bortnik
, and
Y.
Shprits
,
J. Geophys. Res.: Space Phys.
118
,
2196
, doi: (
2013
).
9.
C. F.
Kennel
and
H. E.
Petschek
,
J. Geophys. Res.
71
,
1
, doi: (
1966
).
10.
K. B.
Dysthe
,
J. Geophys. Res.
76
,
6915
, doi: (
2012
).
11.
C.
Crabtree
,
E. M.
Tejero
,
G.
Ganguli
,
G. B.
Hospodarsky
, and
C. A.
Kletzing
, “
Bayesian Spectral Analysis of Chorus Sub-Elements from the Van Allen Probes
,”
J. Geophys. Res.
(unpublished).
12.
O.
Santolík
,
D. A.
Gurnett
,
J. S.
Pickett
,
M.
Parrot
, and
N. C.
Wehrlin
,
Geophys. Res. Lett.
31
,
L02801
, doi: (
2004
).
13.
O.
Santolík
,
C. A.
Kletzing
,
W. S.
Kurth
,
G. B.
Hospodarsky
, and
S. R.
Bounds
,
Geophys. Res. Lett.
41
,
293
, doi: (
2014
).
14.
E. M.
Tejero
,
C.
Crabtree
,
D. D.
Blackwell
,
W. E.
Amatucci
,
G.
Ganguli
, and
L.
Rudakov
,
Phys. Plasmas
23
,
055707
(
2016
).
15.
X.
An
,
B. V.
Compernolle
,
J.
Bortnik
,
R. M.
Thorne
,
L.
Chen
, and
W.
Li
,
Geophys. Res. Lett.
43
,
2413
2421
, doi: (
2016
).
16.
T. M.
O'Neil
,
Phys. Fluids
14
,
1204
(
1971
).
17.
T. M.
O'Neil
and
J. H.
Winfrey
,
Phys. Fluids
15
,
1514
(
1972
).
18.
H. E.
Mynick
and
A. N.
Kaufman
,
Phys. Fluids
21
,
653
(
1978
).
19.
N.
Padhye
and
W.
Horton
,
Phys. Plasmas
6
,
970
(
1999
).
20.
Y.
Omura
and
H.
Matsumoto
,
J. Geophys. Res.
87
,
4435
, doi: (
1982
).
21.
A. R.
Soto-Chavez
,
A.
Bhattacharjee
, and
C. S.
Ng
,
Phys. Plasmas
19
,
010701
(
2012
).
22.
G. J.
Morales
,
Phys. Fluids
23
,
2472
(
1980
).
23.
E. G.
Evstatiev
,
P. J.
Morrison
, and
W.
Horton
,
Phys. Plasmas
12
,
072108
(
2005
).
24.
R. S.
Mackay
, “
Stability of equilibria of Hamiltonian systems
,” in
Nonlinear Phenomena and Chaos
, edited by
S.
Sarkar
(
Hilger
,
Bristol
,
1986
), pp. 254–270.
25.
R. C.
Davidson
,
Methods in Nonlinear Plasma Theory
(
Academic Press
,
New York
,
1972
).
26.
R.
Zhang
,
H.
Qin
,
R. C.
Davidson
,
J.
Liu
, and
J.
Xiao
,
Phys. Plasmas
23
,
072111
(
2016
).
27.
T. H.
Kho
,
A. T.
Lin
, and
L.
Chen
,
Phys. Fluids
31
,
3120
(
1988
).
28.
E.
Ott
,
Phys. Fluids
14
,
959
(
1971
).
29.
O. V.
Agapitov
,
V.
Krasnoselskikh
,
F. S.
Mozer
,
A. V.
Artemyev
, and
A. S.
Volokitin
,
Geophys. Res. Lett.
42
,
3715
, doi: (
2015
).
30.
Y.
Kusama
,
G. J.
Kramer
,
H.
Kimura
,
M.
Saigusa
,
T.
Ozeki
,
K.
Tobita
,
T.
Oikawa
,
K.
Shinohara
,
T.
Kondoh
,
M.
Moriyama
,
F. V.
Tchernychev
,
M.
Nemoto
,
A.
Morioka
,
M.
Iwase
,
N.
Isei
,
T.
Fujita
,
S.
Takeji
,
M.
Kuriyama
,
R.
Nazikian
,
G. Y.
Fu
,
K. W.
Hill
, and
C. Z.
Cheng
,
Nucl. Fusion
39
,
1837
(
1999
).
31.
K.
Shinohara
,
Y.
Kusama
,
M.
Takechi
,
A.
Morioka
,
M.
Ishikawa
,
N.
Oyama
,
K.
Tobita
,
T.
Ozeki
,
S.
Takeji
,
S.
Moriyama
,
T.
Fujita
,
T.
Oikawa
,
T.
Suzuki
,
T.
Nishitani
,
T.
Kondoh
,
S.
Lee
,
M.
Kuriyama
,
JT-60 Team
,
G. J.
Kramer
,
N. N.
Gorelenkov
,
R.
Nazikian
,
C. Z.
Cheng
,
G. Y.
Fu
, and
A.
Fukuyama
,
Nucl. Fusion
41
,
603
(
2001
).