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.

## I. INTRODUCTION

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 $\delta B/B\u223c0.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 whistlers^{5} 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 cone^{6} 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\u22a5>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 Malmberg^{16,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 waves^{19} 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 modes^{22} 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 scattering^{25} 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.

## II. THEORETICAL FORMULATION

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.

### A. Nonlinear model

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.

### B. Wave dynamics

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 $\Phi $, the electrostatic potential, is zero so that, $B=\u2207\xd7A$ and $E=\u2212(1/c)\u2202A/\u2202t$ and Maxwell's equations may be combined to form

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

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 $\chi s$ for species *s* to describe the response of background current to the field

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

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

where

and where $\Omega s=qsB/(msc)$ is the signed particle gyrofrequency for species *s* and $\omega ps2=4\pi 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

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\u2020$. Then we note that Eq. (7) should be viewed as an operator equation^{18} where *ω* is an operator so that the left hand side of Eq. (7) becomes

and using the fact that $el\u2020\xb7D(\omega l0,kl)\xb7el(t)=0$, we arrive at the evolution equation for the mode amplitude

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=\u2212\beta Jl\u2009exp(\u2212i\theta 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=\rho \u0302l\u2009exp(i\phi l)$, where $\rho \u0302l$ and $\phi 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

### C. Electron beam dynamics

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=(\u2212By)ey$, so that the single particle Hamiltonian may be written as

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

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

where

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

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

where $\u2212enbV/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 *i*th beam particle. There are 9 delta functions and 6 integrals, which leave three delta functions. Here we write

### D. Self-consistent dynamics

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

where

putting in our representation for *A _{l}* and equating real and imaginary parts of Eq. (18), and after some straightforward algebra we can write the evolution equations for $Jl$ and

*θ*as

_{l}and

### E. Hamiltonian formulation

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, $\theta 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

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, $\theta \u0307l=\u2202H/\u2202Jl$ and $J\u0307l=\u2212\u2202H/\u2202\theta l$ will give the appropriate wave dynamics if we choose

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.

## III. SINGLE PARALLEL WHISTLER WAVE

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

which can be solved for *k*^{2} as

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

so that

where *S* + *H* = *R* and

With two species and in the limit $\omega \u226b\Omega i$ and with $me/mi\u226a1$, the dispersion relation can be written as

where $\omega \xaf=\omega /\Omega e,\u2009k\xaf=ck/\Omega e$, and $\delta 2=\omega pe2/\Omega e2$. The self-consistent Hamiltonian reduces significantly to

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, $\psi i=kzi\u2212\omega kt\u2212\varphi i$. To take advantage of this, we introduce the type-2 mixed-variable generating function that takes us from $(t;zi,\varphi i;P\varphi i,pzi)$ to $(t;\Psi i,\xi i;P\Psi i,Ii)$ via

so that the transformation equations are

note that the generating function is time-dependent, thus the new Hamiltonian becomes $H\u2032=H+\u2202tF2$. The new Hamiltonian becomes

This Hamiltonian has *M* + 1 degrees of freedom, *I _{i}* are constants of the motion. We can make this dimensionless

where time is measured in Ω_{e}, energy with *mc*^{2}, particle momentum $P\psi i$ and wave action with $mc2/\Omega e$, and the new constant $\alpha =\delta 2(nb/n0)/(\Omega e/c2\u2202\omega \xafD)$. 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

From momentum conservation, we see that as the beam particles lose $p\psi 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 $\omega =0.225\Omega e$, and the parallel wave vector was, $kc/\Omega e\u22432.769$. The background plasma parameters were chosen typical of chorus so that $\delta \u22435.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 ($10\u221213$%). We initiated the particles around the fixed point (as discussed in Section III A) with $\psi i0=2\pi 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.

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\u2009\Omega 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\u2009\Omega 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\psi =0.04$. Conservation of momentum providing the mechanism of saturating the instability. In the third panel, at $t=180.0\u2009\Omega 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 model^{18} may be a fair approximation for times just after the nonlinear saturation. However, in the final panel, at $t=2500.0\u2009\Omega 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 ($\theta \u0307$ and the instantaneous growth rate $J\u0307/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.

### A. Hamiltonian analysis of fixed points

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=2J\u2009cos(\theta ),\u2009Y=2J\u2009sin(\theta )$, and *X* is a momentum variable and *Y* a position variable so that the Hamiltonian becomes

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

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=(\delta \psi i,\delta p\psi i,\delta Y,\delta Y)T$ and calculate the linearized motion around the fixed point, which may be written in the matrix form as $z.=Lz\u0307$, where

**L**is a block matrix of the form

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 $\xb1\sigma j$, imaginary pairs $\xb1i\sigma j$, zero, or complex quadruplets $\xb1a\xb1ib$. The solution to the linearized motion can be written as a sum $\u2211jcj\xi j\u2009exp(\sigma jt)$, where

*ξ*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.

_{j}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 $\delta =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 *s _{j}* which was computed by

where $\omega $ 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.

### B. Linear theory

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

where $\omega \xaf,k\xaf$, and *δ* are defined before. $V||,V\u22a5$ are the parallel and perpendicular beam velocity.

In Figure 5, we plot the linear growth rate Im$(\omega \xaf)$ as computed from Eq. (40) as a function of $k\xaf$ for two cases. Both cases consider a 50 keV beam propagating at a pitch angle of 17°, through a plasma with $\delta =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\xaf=2.04$ is indicated by a vertical grey line. As dictated by the nonlinear model and the requirements of the fixed point, this $k\xaf$ and corresponding $\omega \xaf$ 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.

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 $\delta =5.115$ with a beam density of $nb/n0=0.01$. We choose a $\omega \xaf=0.31407694$ and $k\xaf=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 $\omega \xaf$ and $k\xaf$ 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($\omega \xaf$) = 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.

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 experiments^{14} 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.

In Figures 9 and 10, we show the evolution of the instantaneous frequency ($\theta \u0307$ and the instantaneous growth rate $J\u0307/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}

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=14\u2009500\Omega 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=15\u2009000\Omega e$, these bar like structures have detached from each other, and by $t=25\u2009000\Omega 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.

## IV. MULTIPLE PARALLEL WHISTLER WAVES

### A. Laboratory experimental results

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\xb0$, ratio of beam to background plasma density $nb/n0=0.05\u22120.13$, and ratio of plasma to cyclotron frequency $\delta =1.75\u22129.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.

### C. Self consistent Hamiltonian

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

Note that the presence of at least two wavelengths with a gyrophase dependence $\varphi 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=(\omega 2\u2212\omega 1)\tau +\theta 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

where $\Delta k=k2\u2212k1$. In addition to *h*, examination of the equations of motion reveals that there are two additional constants of motion

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 $\alpha 22\u226a\Delta \omega $

where $\Delta \omega =\omega 2\u2212\omega 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.

### D. Nonlinear harmonic generation

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.

To model the harmonic generation seen in space observations, we choose the background plasma to have $\delta =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\xd710\u22124$. The unstable wave that resonates with this beam has a frequency of $\omega \xaf1=0.123$ and a wavelength of $k\xaf1=1.92$. We consider a second wave with $k\xaf2=2k\xaf1$ and a frequency $\omega \xaf=0.358$ that satisfies the dispersion relation. With these parameters, the dimensionless constants $\alpha 1=8.74\xd710\u22123$ and $\alpha 2=6.38\xd710\u22123$. 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 *J*_{1} (blue) and *J*_{2} (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 ($J\u223c10\u22123$), 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.

### E. Nonlinear evolution of two closely spaced waves

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 $\delta =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 $\omega \xaf=0.318$ and a wave-vector of $k\xaf=3.508$, while the non-resonant mode has a frequency of $\omega \xaf=0.402$ and a wave-vector of $k\xaf=4.21$. Both modes are linearly unstable. In Figure 15, we plot the evolution of the wave-action variables *J*_{1} and *J*_{2}. 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.

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.

## V. DISCUSSION

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 fields^{22}). 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-potential^{28} 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.

## ACKNOWLEDGMENTS

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.