Annales Geophysicae On the numerical modelling of VLF chorus dynamical spectra

This paper presents a study of the use of a one-dimensional Vlasov Hybrid Simulation (VHS) computer code to simulate the dynamical spectra (i.e. frequency versus time spectrograms) of ELF/VLF chorus signals (from ∼a fraction to∼10 kHz). Recently excellent measurements of chorus have been made in the source region close to the geomagnetic equator aboard the four spacecraft Cluster mission. Using Cluster data for wave amplitude, which is up to 300 pT, local gyrofrequency, cold plasma density, and Lshell, observed chorus signals are reproduced with remarkable fidelity and, in particular, sweep rates in the range 1– 10 kHz result as observed. Further, we find that the sweep rate is a falling function of increasing cold plasma density, again in accord with observations. Finally, we have satisfactorily simulated the rather rare falling frequency elements of chorus which are sometimes observed aboard Cluster in the generation region. For both rising and falling chorus we have presented detailed structural analyses of the generation regions. The main contributor to the frequency sweep rate is primarily the establishment of wave number/frequency gradients across the generation region by the out of phase component of the resonant particle current. The secondary contributor is the shortening of the wavelength of resonant particle current relative to that of the wave field. In view of the close agreement between observation and simulation, we conclude that nonlinear electron cyclotron resonance is indeed the mechanism underlying the generation of chorus signals just outside the plasmasphere.


Introduction
Chorus, a fascinating phenomenon in radio physics, consists of naturally occurring, strong ELF/VLF discrete radio emissions propagating in the whistler mode through the Earth's magnetosphere (Helliwell, 1965).Their frequencies generally lie between a few hundred Hz and several kHz, and they are usually found outside (beyond) the plasmapause (Lemaire and Gringauz, 1998, p. 107) at L-shells of L=3-10 between local midnight and local noon.Chorus signals are commonly observed on the ground at such L-values, and excellent examples are seen at Sodankyla Geophysical Observatory (Manninen, 2005) and Halley Bay, Antarctica (Smith et al., 1991), following geomagnetic disturbances.Chorus emissions are commonly observed on space missions, particularly on board the Cluster spacecraft (Santolik et al., 2002), Geotail (Nunn et al., 1997) and the MAGION satellites (Kozelov et al., 2001), and going back to the early space missions such as OGO3 (Burtis and Helliwell, 1976) and Injun 3 (Oliven and Gurnett, 1968).Chorus has also been observed in the Jovian magnetosphere (Coroniti et al., 1984), as well as the Kronian magnetosphere (Hospodarsky et al., 2007).The reader will find a good overview of the phenomenon of chorus at Extremely Low Frequency (ELF, <3 kHz) and Very Low Frequency (VLF, 3-30 kHz) in the paper by Sazhin and Hayakawa (1992).
ELF/VLF chorus signals, "chirps" sounding like birds at the time of the dawn chorus, typically take the form of a sequence of rising frequency tones, with absolute values of the sweep rates being in the range of 1-10 kHz/s (Helliwell, 1965, p. 240-248).The sequence of elements may be tightly packed or they may often be well spaced, especially when observed at higher geomagnetic latitudes.The time separation may be fairly constant, giving a regular sequence of emissions, or rather random, giving what looks like a sequence of discrete and well separated emissions.Chorus is generally believed to be generated by a cyclotron instability in the neighbourhood of the geomagnetic equator by gyroresonant energetic electrons and whistler mode waves moving in opposite directions along a geomagnetic flux tube, typically with L∼4 to 6 (Burton and Holzer, 1974;LeDocq et al., 1998).Although there have been some controversies amongst theoreticians in the past about the generation mechanism, there is no doubt that the phenomenon is fully nonlinear (Trakhtengerts, 1999).Nonlinear cyclotron resonance is exemplified by the process of "phase trapping" (see, e.g., Chen, 1974, p. 209), with some electrons being trapped in the potential well of the electromagnetic wave (Nunn, 1974).It was shown in that paper that the magnetosphere is an inhomogeneous medium, due to the parabolic variation of the geomagnetic field strength with distance from the geomagnetic equator, and also due to the sweeping frequency of each element.The consequence of these facts is that the cyclotron resonance velocity is a function of space and time along a geomagnetic flux tube.This profoundly affects the resonant particle (electron) dynamics such that the transfer of power between particles and wave fields becomes very much more effective than for a linear interaction, the nonlinear growth rates being several times (∼2-6 times) the linear growth rate.The underlying reason for this is that in an inhomogeneous medium a much broader slice of electron parallel velocities may interact with the wave than in a homogeneous medium (Nunn, 1974).
In the frequency-time plane, Fig. 1 shows the dynamic spectrum df/dt of one chorus element triggered near the equatorial plane of the magnetosphere at a particular frequency.This could be the frequency corresponding to a particular magnetospheric line or line of Power Line Harmonic Radiation from the electricity grid system (see, e.g., Manninen, 2005) or the frequency where the growth rate of the cyclotron instability is a local maximum.Immediately following the rising frequency tone, there is a quiet band; that is discussed later in this paper.Shortly after that, the second element commences.
Theoretical models of chorus generation typically assume, for simplicity, that the ELF/VLF wave of interest is at all times propagating parallel to the geomagnetic field, at least in the equatorial generation region.There is a good deal of experimental evidence to support this idea (Hayakawa et al., 1984;Goldstein and Tsurutani, 1984).In the closely related problem of VLF emissions triggered by the signal from a ground based radio transmitter, such as the one at Siple Station, Antarctica, parallel propagation is assumed to arise due to ducting of the wave in cylindrical enhancements of plasma density (known as "ducts") within the plasmasphere.However, outside the plasmasphere, ducting is less common and the quasi parallel propagation of chorus in the region of generation is presumed to be due to the fact that parallel propagating waves are preferentially nonlinearly amplified.
The four spacecraft of the successful Cluster mission have provided excellent observations of chorus close to the source region (Parrot et al., 2003), at L-shells ∼L=4.4,with both the temporal and spatial resolution of the measured quantities being the best available to date.Simultaneous measurements have been made of the cold plasma density, the ambient geomagnetic field strength and the wave amplitude, as well as presentations of the frequency-time spectrograms (the dynamical spectra of the signals) showing clearly the sweep rates of the emissions.
The purpose of this paper is, under the assumption of parallel propagation and using a 1-D Vlasov Hybrid Simulation (VHS) code (Nunn, 1993), to confirm the nonlinear waveparticle interaction theory as applied to VLF chorus generation by modelling as closely as possible the frequency-time characteristics (i.e. the dynamic spectra) of specific chorus events as observed by Cluster.

The CLUSTER observations
Since its inception in 2001 the Cluster mission has provided impressive data allowing theoreticians to gain a much better understanding of the generation mechanism of chorus.Close to perigee, the Cluster spacecraft pass through the equatorial zone at L∼4-5 on the night or morning side of the Earth where they are ideally placed to observe chorus in its generation region.Here the individual Cluster satellites have separations of hundreds to thousands km; they are therefore in a position to investigate the three-dimensional structure of the chorus generation region.Particularly at magnetically disturbed times Cluster observes intense VLF chorus signals.High resolution electric and magnetic field waveforms are recorded by the wide band data (WBD) wave instrument, using the 88-m long electric dipole antenna (Gurnett et al., 2001).The wave magnetic field is measured by the STAFF (Spatio-Temporal Analysis of Field Fluctuations) instrument (Cornilleau-Wehrlin et al., 2001).The use of both electric and magnetic field measurements enables the wave Poynting vectors to be determined.The ambient cold plasma density is estimated from measurements made by the WHISPER sounder (Decreau et al., 2001) and the ambient geomagnetic field is observed using a fluxgate magnetometer (Balogh et al., 2001).The reader is referred to Santolik et al. (2002) and Santolik and Gurnett (2003) for a thorough explanation of the detailed analysis of these excellent sets of data.At times Cluster observes "upper band chorus", just above half the local equatorial gyrofrequency, which often disappears as the spacecraft move away from the equator.More frequently observed is the "lower band chorus" somewhat below half the electron gyrofrequency, w be /2.We shall concentrate on this type of emission in our numerical simulations.
Calculations of the Poynting flux on the four spacecraft allow the spatial extent of the chorus generation region (GR) to be mapped out (Santolik et al., 2003;Santolik and Gurnett, 2003).It appears that in the GR wave propagation is often close to being parallel to the geomagnetic field (the angle between the wave normal k and the geomagnetic field B, θ<10 • ).The generation region extends downstream (in the direction of group velocity) from the equator some 2000 km, approximately (Santolik et al., 2004).By examining corre-lated spectrogram features between the different spacecraft, some indication of the lateral extent of the GR may be obtained (Santolik et al., 2004).The width of the generation region across the geomagnetic field direction turns out to be ∼50-100 km, a result confirmed by Cluster data analysis in Inan et al. (2004) and Platino et al. (2006).The wave field envelope often consists of a sequence of "wave packets" of duration from a few ms to tens of ms.Such an envelope indicates the presence of strong sidebands with a separation ∼100 Hz.It is noted, though, that the sideband structure has a low correlation perpendicular to the ambient magnetic field (Santolik et al., 2004).Perhaps one of the most significant Cluster results found to date is the overall level of the wave amplitude, with maxima ∼200-300 pT being common (Santolik et al., 2004).At such large amplitudes strong trapping of cyclotron resonant electrons by the wave is expected, providing powerful support for the hypothesis that nonlinear electron cyclotron resonance trapping is the underlying plasma physical generation mechanism.Santolik (2008) has usefully reviewed these results, putting them in the current research context.
For simulation purposes we now concentrate on two specific data examples.Figure 2 shows the spectrograms on the four Cluster satellites at L=4.4, at 08:49 UT, on 18 April 2002 (Santolik et al., 2003).It was found that the wave vectors were closely parallel to the ambient geomagnetic field, with θ <10  the upper chorus band is largely non existent.The Cluster 1 and Cluster 2 spectrograms are very alike.Those from satellite 3 and 4 are more different because these spacecraft are separated more across (perpendicular to) the geomagnetic field line.Individual chorus elements, whose frequencies lie somewhere between 2.7 and 3.8 kHz (i.e. between 0.32 and 0.46 times the local electron gyrofrequency), are randomly spaced, some being quite isolated.The sweep rates are remarkably consistent between different chorus elements, and are in the region of +10 kHz/s.The maximum observed wave amplitudes are ∼30 mV/m, corresponding to 300-500 pT.
The second example shown in Fig. 3 is at 08:52 UT, some three minutes later on the same day.The lower band chorus is at lower frequencies here, and distinct upper band chorus is present, although embedded in a hiss band.The lower band elements are very clean and again well separated.From the point of view of their theoretical explanation some elements are indistinguishable in form from discrete, or triggered, emissions.The reader will find exact amplitude and frequency profiles of the elements indicated by the black arrows in Santolik et al. (2003).Here we consider that the maximum envelope amplitude is ∼25 mV/m and the sweep rates are again ∼10 kHz/s.In both examples the ambient plasma density derived from the WHISPER sounder is in the region of 2-4 electrons/cc, a rather low value.Clearly the spacecraft were outside the plasmapause.
More recently the Cluster analysis team (Macusova et al., 2007) has undertaken a detailed analysis of sweep rate probability distributions for the observed chorus, for both upper and lower band chorus.They have found that, significantly, for both rising and falling frequency chorus signals, the sweep rate increases with decreasing ambient plasma density.

Previous theoretical work on chorus
In discussing theories of the generation of ELF/VLF chorus, we first address the problem of triggered and discrete VLF emission generation.This is a very closely related problem since, to a first approximation, chorus may be viewed as a succession of discrete emissions.We are then left with only the issue of how the triggering of one element is controlled by the previous element, and with discovering the nature of the interaction between one element and the next, via the nonlinear wave-particle interaction process.
Theoretical work on triggered VLF emissions is usually directed at processes occurring inside the plasmapause, where such events are believed to take place in ducts of enhanced cold plasma density.Ducting keeps the wave vector essentially parallel to the ambient geomagnetic field vector, and all theories of triggered VLF emissions assume parallel propagation (Omura et al., 1991).Some fascinating work has been done on nonlinear particle dynamics for oblique whistler signals (Bell et al., 1982), but this work has not yet progressed to such a stage as to be able to consider the emission triggering process.
There seems to be general agreement that the chorus generation mechanism is due to the nonlinear electron cyclotron resonance between the narrow band whistler mode wave and energetic (some tens or hundreds of keV) electrons (Nunn, 1974;Omura andMatsumoto, 1982, 1985).Nonlinearity in the resonant charged particle trajectories is equivalent to the familiar concept of "phase trapping" in which the phase angle between the perpendicular velocity vector and the electric wave field vector oscillates about a phase locking angle (Nunn, 1990).The electron parallel velocity similarly oscillates about the locally defined cyclotron resonance velocity v res =(ωe )/k.
The medium is inhomogeneous (i.e. the gyroresonance velocity v res is a function of z and time t) mainly through the (approximately) parabolic variation of ambient geomagnetic field strength with distance z away from the equator, but also -and significantly -through the frequency sweeping of the emission itself (Nunn, 1974;Omura et al., 2008).An early phenomenological theory due to Helliwell (1967) identified the importance of the inhomogeneity, and identified the region of generation as being the point of net zero inhomogeneity; he introduced the term "second order resonance".However, the generation region is of finite size; also the explicit behaviour of the gyroresonant electrons was not considered in that paper.Indeed, wave-particle interactions are most effective when the inhomogeneity ratio |S|∼0.3-0.8, and not zero.(S is defined as the sum of all inhomogeneities divided by wave amplitude (Nunn, 1990;Omura et al., 2008)).For a useful review of competing theories, the reader may consult Omura et al. (1991).
The situation of nonlinear trapping dynamics in an inhomogeneous medium was investigated by Nunn (1972Nunn ( , 1974)).It was found that in an inhomogeneous medium for "trapped" electrons the phase of the perpendicular velocity relative to the wave electric field, ψ, oscillates about the phase locking angle P 0 =cos −1 S, with |S|<1, and that such particles undergo a steady change of energy and magnetic moment.Using Liouville's theorem to calculate the electron distribution function F (see Krall and Trivelpiece, 1973, p. 447), this results in the trapping region in velocity space being characterised by a large value of |dF |; it can be either a "hole" or a "peak", depending on the sign of the inhomogeneity ratio S.
The resulting resonant particle current J res will have its phase controlled by the phase locking angle, and thus by the inhomogeneity ratio S which varies in space and time.The local nonlinear growth rate is proportional to −J r , where J r is the component of resonant particle current parallel to E. The nonlinear growth rate is proportional to the linear growth rate, increasing to multiples of the linear growth rate with increasing trapping time.This real part of the resonant current provides the power to produce the chorus emission.The component J i parallel to the wave B field changes the phase of the wave field, and thus directly determines the observed frequency shifts.It should be noted that both linear and nonlinear systems require a source of free energy, namely an anisotropic zero order distribution function.
Trapped particles undergo relatively large energy and magnetic moment changes.In a chorus generation region these changes will normally be positive; there is therefore the possibility of very effective acceleration and/or heating of electrons by chorus.In a positive inhomogeneity region, namely on the upstream side of the equator, trapped particles will be de-energised.At this point we remark that untrapped, or so called "passing particles", that "pass through" the resonance condition, make a significant contribution to the overall energy transfer to the wave.It is quite wrong to assume that trapped particles are the sole contributors to wave growth as far as energy changes are concerned.
In the paper by Nunn (1974) the field equation was integrated using a model for J res assuming the dominance of trapped particles in the distribution function dF.Rising frequency emissions triggered by the ground-based VLF transmitter NAA were simulated.More recently a full Vlasov Hybrid Simulation code (VHS) has been developed that has simulated triggered risers, fallers and hooks observed at Halley Bay in Antarctica (Nunn and Smith, 1996;Smith and Nunn, 1998) and Power Line Harmonic Radiation (PLHR) triggered emissions recorded at Sodankyla Geophysical Observatory in Finland (Nunn et al., 1999;Manninen, 2005).These simulations pointed to the distinctly different spatial structures for riser and faller generation regions (GRs), these structures resembling quasi-static, yet fully nonlinear, self consistent entities that behaved in a very repeatable and quite stable manner.Hooks were interpreted as resulting from transitions between faller and riser type GRs.
Most recently the VLF emission code has been run, systematically, to investigate how the triggered emission form and sweep rate depend upon key parameters, such as the linear growth rate, cold plasma density, saturation wave amplitude, input wave amplitude and pulse length (Nunn et al., 2005).One result emerging from that study of considerable relevance to the chorus generation problem was that it was found that the frequency sweep rate was largely determined by the cold plasma density, becoming larger with decreasing density, for both risers and fallers.Risers were the commonest emissions, fallers being associated with either weak growth rates ("termination" fallers), or very strong growth rates, where oscillating tones sometimes appeared.Some recent very significant studies (Katoh and Omura, 2007b;Omura et al., 2007) have discovered that, in the case of cyclotron resonance with highly relativistic electrons in a parabolic inhomogeneous magnetic field, the resonance velocity can become positive (downstream, in the direction of group velocity), so allowing very long trapping times and very significant energisation of the electrons.This result has great implications for electron acceleration and heating by chorus and also for the wave-particle interaction process itself.These are yet to be fully explored.
We now turn to the chorus problem, which is more challenging than the triggered emission problem as we are dealing with packed successions of emissions.A global view of these phenomena was expounded by Trakhtengerts (1995), in which the source mechanism was viewed as being due to a Backward Wave Oscillator (BWO) instability.Recent work has shown good agreement between BWO theory and D. Nunn et al.: On the numerical modelling of VLF chorus dynamical spectra many details of the Cluster observations (Trakhtengerts et al., 2007).A feature of the BWO theory is that it postulates the existence of a "step" in the equatorial parallel velocity in the zero order electron distribution function, which might result from quasilinear diffusion at the upper frequency edge of a hiss band.Such a step gives a greatly enhanced peak in the linear growth rate for VLF waves whose resonant velocity coincides with the step velocity.This provides not only a natural starting frequency for each chorus element but also supplies the large anisotropy and initial growth rates required by the emission process.The recent monograph by Trakhtengerts and Rycroft (2008) goes into this and related mechanisms in considerable detail.Both Kozolev et al. (2008) and Santolik et al. (2008) have compared Cluster data with results of the BWO model, and found good quantitative agreement.
In the paper by Nunn et al. (1997), data from Geotail satellite observations at L=10 were input into the VHS code, and rising and falling chorus signals were accurately simulated.However, the chorus elements simulated were not closely packed; rather, they were more like individual emissions.One significant aspect of the Geotail data was that the observed amplitudes were as large as 160 pT, which implies strong nonlinear trapping.
Recently, important simulations of chorus producing sequences of narrow band rising tones have been achieved with a numerically intensive broadband Particle in Cell (PIC) simulation code developed at Kyoto University (Katoh andOmura, 2006, 2007a;Omura et al., 2008).This code has a bandwidth encompassing that of chorus; it can thus address the issue of the interaction between successive chorus elements via nonlinear wave-particle interactions and how one element controls the timing of the next element (see Fig. 1).It also addresses the question of why highly structured emissions are produced as a result of the plasma instability in preference to broadband hiss.It may be of some interest to compare this work with that of the present paper.
The broadband property of the electron hybrid code (Katoh and Omura, 2006Omura, , 2007a;;Hishikima et al., 2009) avoids the need for matched filtering and the simulations show saturation.However, the PIC methods are noisier than Vlasov methods and derive the distribution function from the density of simulation particles, in contrast to the Vlasov method which resolves the distribution function explicitly.The Vlasov code here uses plasma physical data direct from Cluster and aims to verify fundamental theory by performing a close simulation of actual Cluster events.The electron hybrid code (Katoh andOmura, 2006, 2007a) , however, applies a spatial rescaling of the wave magnetic field by a factor ∼12 in order to speed up the computation.Analysis of computer simulation results show that a narrow region near the equator acts as a source of the VLF wavelets, and it is here that the wave number gradient across the generation region is set up and the initial frequencies of the seed wavelets determined (Hishikima et al., 2009).Omura et al. (2008) introduce, for the first time, full relativistic trapped particle dynamics.They also derive an expression for the frequency sweep rate which is proportional to wave amplitude at the equator.They argue that maximum nonlinear growth rates are achieved when the net inhomogeneity factor S at the equator is −0.4 for rising chorus, assuming the average perpendicular velocity for the energetic electrons.This is an interesting estimator of df/dt, and future research might investigate the extent to which both actual chorus and simulations satisfy the Omura et al. (2008) expression.In a recent paper by Hikishima et al. (2009) PIC simulations of chorus generation showed very close agreement with the Omura expression for df/dt.
In recent work by Inan and Platino (Inan et al., 2004;Platino et al., 2006) analysis of Cluster data shows clearly identical elements on different spacecraft but with significant frequency shifts up to 1 kHz.They have put forward a model in which the generation region has dimensions transverse to the field line of order wavelengths, and highly localised along the field line with a scale length ∼100s kms.The model identifies observed frequency shifts with Doppler shifts of fast moving generation regions.

The computational model
In this research we shall use a Vlasov Hybrid Simulation (VHS) computer code that has been used to simulate triggered VLF emissions (Nunn, 1990(Nunn, , 1993;;Nunn et al., 1997Nunn et al., , 2005)).The code has been modified to enable it to model a sequence of chorus elements and to conform to the plasma parameters as observed by the Cluster mission.
The chorus generation region, as observed by Cluster (Santolik et al., 2003) in the equatorial plane of the magnetosphere, is a three-dimensional region that is "cigar shaped" and aligned with the geomagnetic field.This extends from near the equator to about 2000-3000 km downstream.Performing a full 3-D simulation is not feasible with the computer resources available anywhere in the world at present.Thus the simulation has to be performed in 1-D, considering parallel propagation, which is in good agreement with the predominantly quasi parallel propagation deduced from the Cluster observations (Santolik et al., 2003) made in the vicinity of the chorus generation region.Such an approach is a close approximation to the actual physical situation.In recent work Platino et al. (2006) and Inan et al. (2004) have analysed Cluster data and shown that the generation regions of chorus are very localised 3-D objects.Thus the 1-D simulation ignores the loss of wave energy due to radiation by the resonant particle current fields into modes with k vectors inclined at an angle to Bo.Here, such a nonlinear "spreading loss" has to be modelled phenomenologically, as an extra loss term, when the wave amplitude exceeds the predetermined saturation value B sat .Similarly, Landau damping arising from non parallel propagation is, to date, ignored.
The ambient geomagnetic field Bo(z) appropriate to L=4.4 has a parabolic variation with z, the distance from the equator (at z=0), where, in dimensionless units (Nunn, 1990), Numerical values are all given in SI units.The cold plasma density N e (z) is also taken to vary parabolically with z: The bare essentials of the code may be stated as follows.The phase space box and particle grid are defined in z, V * z , ψ space for each value of energetic electron pitch angle, where V * z is the electron velocity parallel to the geomagnetic field and ψ is the gyrophase angle between the perpendicular velocity vector and the wave E field.
The phase box is fixed in variable z and more than encompasses the region where electron trapping is possible by a whistler mode wave at the saturation amplitude.The range of V * values is centred on the local resonance velocity V res (z, t) and has a width ∼3 trapping widths (Nunn, 1990).In this application we use a spatial grid of 2048 points, a phase grid of 16 points, a V * grid of 40 points.Significant nonlinear wave-particle interaction occurs for quite a narrow range of V ⊥ , with electron pitch angles of from 40-60 degrees.This means that in a Vlasov formulation only a few grid points in V ⊥ are needed; we shall employ three or indeed only one.Simulation particles fill phase space with a density which must be greater than ∼1.2 per phase space cell.
The motions of these electrons are followed forwards continuously in time, integrating the relativistic equations of motion to second order accuracy.Particles leaving the phase box are deleted from the simulation and new particles are continuously introduced at the V * phase box boundary where the phase fluid is flowing in/out.The integrated energy change dW is calculated for each particle; applying Liouville's theorem, the distribution function F is then known at all the phase space points occupied by the particles.In order to calculate the resonant particle current, it is necessary to interpolate F from the particles to the phase space grid, using a simple low order interpolator where the sum l is over all simulation particles in the cells touching the grid point in question, and the β ij kl are the familiar "area weighting coefficients" used in PIC codes.The current is then readily obtained by integrating over the velocity space grid.The VHS method has many advantages, such as: 1.It is inherently low noise, and far superior to PIC methods in this respect (Cheng and Knorr, 1976).
2. The population of particles is dynamic, and is always confined to particles which are in resonance with the local wave field.
3. The method is very stable and robust against filamentation of the distribution function.
4. No artificial filtering or smoothing of the distribution function needs to be introduced.

5.
The method provides good diagnostics of the resonant particle distribution function.
6.It is permissible to use very few, or only one, value of perpendicular velocity, which is valid if the contribution to the linear growth rate as a function of perpendicular velocity is sharply peaked; running PIC codes with single values of perpendicular velocity implies a ring distribution, which is a quite different assumption.

Development of the field equation
In order to understand fully the computational results and the frequency shift mechanism, it is valuable here to derive the one-dimensional (1-D) field equations integrated by the code, under the assumption of parallel propagation (Nunn, 1990).With z being the spatial coordinate along the geomagnetic field line measured from the equator, the electron plasma frequency is given by where both the cold plasma density N e (z) and the electron gyrofrequency vary parabolically with z.The electron gyrofrequency is given by Helliwell (1965) as To derive the wave field equations we employ dimensionless units as follows: A dimensionless electric field amplitude E is defined in terms of the amplitude E (in V/m) by the expression and the dimensional resonant particle current by Assuming strictly parallel propagation, and ignoring terms ∼1/R E , where R E is the radius of the Earth, Maxwell's equations and the linear equations of motion of the ambient cold electrons give the following field equation in dimensionless units (dropping the primes henceforth): The dimensionless dispersion relation now becomes where k 0 (z) and ω 0 are the dimensionless frequency and dimensionless wave number, and the dimensionless group velocity is of course given by For a band limited ("narrow band") wave-particle interaction, with a bandwidth ∼150 Hz, we now define a base frequency ω 0 and the corresponding wave number k 0 (z) through the above dispersion relation.It is evident that, where the frequency of the simulated emission is rising or falling, the base frequency will have to be constantly redefined.This is achieved by performing a Fast Fourier Transform on the complex amplitude field, determining its mean frequency, and then transferring the accumulated phase from the complex amplitude to the base phase.Next, we define the dimensionless amplitude and resonant particle phasors R, J in which the fast phase variation is divided out where the phase appropriate to the base frequency is Under the narrow band assumption, the field and current phasors are slowly varying in space and time, and the general field equation above reduces to the form derived by Trakhtengerts and Rycroft (2000) (their Eq. 7) The last term above represents the nonlinear spreading loss which is zero below the saturation level and increases linearly above it.Defining again J r (z, t), the component of resonant particle current parallel with the local wave electric field, and J i (z.t) the component parallel with the local wave magnetic field, are If we express the wave phasor as where R(z, t) is the wave amplitude and N (z, t) is the wave field additional phase, then the field equations may be expressed in the following useful form (Nunn, 1990;Trakhtengerts and Rycroft, 2000) Clearly the in phase component of the resonant particle current J r amplifies the wave amplitude and provides the growth/power input into the local wave field, and the out of phase current J i modifies the local additional phase.Resonant particle trapping in an inhomogeneous medium produces a very prominent J i field (Nunn, 1990); it is this that causes the frequency sweep of triggered VLF emissions and of chorus.However, this may only be demonstrated convincingly by a fully self consistent simulation, as is done here.Straightforward manipulation of Eq. ( 19) gives an expression for the localised rate of change of wave frequency (Nunn, 1990) The first term represents purely the advection of higher frequencies/shorter wavelengths from upstream.The sum of these three terms must be roughly constant through the generation region and outside it.The second term represents a wavenumber/frequency shift due to the differential rotation of the wave field vector from the gradient with respect to z of the quantity J i /R.From another viewpoint it can be seen that the wavelength difference between the resonant particle current and the wave field causes the steady frequency change.Numerical experimentation (Nunn, 2005;Omura et al., 2008) revealed that the advective term appears to be dominant.Such wave number gradients have to be set up naturally as a result of the self consistent interaction involving J i (z, t).Omura et al. (2008) show that such gradients are established in the equatorial region, probably due to detrapping of trapped particle bunches.The last term, which corresponds to a non linear modification of the dispersion relation, on its own can never be very significant since, on integration, it can only give frequency shifts ∼J i/Rω 0 V g /k 0 , which is ∼20 Hz.However, this term will not be negligible at any given point due to the variability of the other two.
In this code the wave field is band limited with ∼150 Hz bandwidth at each time step.The reason for this is as follows.Without filtering the wave spectrum outside the desired band increases steadily.This unwanted wavefield is not provided with resonant particles and is not handled properly by the grid provided, and it would cause numerical diffusion in the integration of resonant particle trajectories.Now the wavefield must be allowed to develop wavenumber gradients since these are significant in providing the sweeping frequency.At each time step the developed wavenumber in the whole complex wave amplitude field is determined as a function of z.A least means squared linear fit to this function is determined in which the penalty function is weighted by the magnitude of the wave amplitude.This prevents developed wavenumbers in the upstream part of the phase box, where wave amplitudes are small, from influencing the fitting process.The upstream boundary condition cannot influence the developed gradient as once the emission has started the incoming wavefield is zero.The reader should note that any developed wavenumber/frequency gradients and the resulting frequency sweep rates are determined by the main generating region where the wave field amplitudes are significant.The phase due to this measured wavenumber gradient θ(z)=exp(i(dk/dz)z 2 ) is then subtracted from the wavefield phase, and then the resultant is filtered by a spatial Discrete Fourier Transform (DFT) and IDFT, where the filter frequency response is boxcar fitted to the desired passband width.The phase θ (z) is then added back to this complex field.This matched filtering process corresponds to a passband filter in which the passband centre is a linear function of position.Of course, if the simulation bandwidth were wide enough to accommodate the range of frequencies present in the spatial box at any one time (∼1000 Hz), then matched filtering would not be necessary.However, the computer run times would then be ∼1000 times longer.
Note that the code does not integrate the wave amplitude and phase update equations (Eqs.18 and 19) separately; it is the complex field Eq. ( 15) which is integrated as follows.A separate field grid is defined with separations V g t, where t is the simulation time step and V g the initial wave group velocity.At each step the field values are advanced by one grid position and updated using a weighted sum of current values calculated on the counterstreaming particle grid.

Production of sequences of chorus elements
The basic code simulates each element in turn, having a bandwidth ∼150 Hz.We postulate the following model for producing sequences of elements, as illustrated in Fig. 1.Following the work of Trakhtengerts (1995) and co-workers, we assume the existence of a diffuse step in the equatorial parallel velocity at a value corresponding to the gyroresonance velocity for the start frequency of each emission.This gives a very enhanced linear growth rate, relative to the energetic particle flux available, and thus defines the natural starting frequency as being that of maximum instability.
In this work we select a zero order distribution function Fo (µ, W ) consisting of multiple bi-Maxwellians with anisotropy A∼2, with a superimposed smoothed step in the equatorial parallel velocity V z:  ..
where d=0.7, σ =0.04 are suitable numerical values.Here µ is the dimensionless magnetic moment, W is the dimensionless energy, V res is the resonance velocity at the equator at the base frequency, V z is the parallel velocity, and the Tj are the corresponding parallel and perpendicular temperatures.The coefficients Aj are the amplitudes of the bi-Maxwellians, selected so that the linear growth rate at the equator has the specified value.Figure 4 plots the zero order distribution function as a function of equatorial parallel and perpendicular velocities.The diffuse step is only just visible as an inflection in the contour lines.For Case 1 of the simulations to follow Fig. 5 plots the linear equatorial growth rate as a function of normalised frequency, relative to the equatorial gyrofrequency.The growth rate itself is normalised to its initial value, which in this case is 1200 dB/s at a normalised frequency of 0.35.The growth rate is peaked just above the starting frequency as expected.As the emission element frequency moves away from the starting frequency, the linear growth rate, and hence also the nonlinear growth rate, decrease until they fall below the level at which a self sustaining generation region is possible.This element then "dies", all the wave energy having propagated out of the equatorial zone.The base frequency is then reset to the start frequency, so allowing a new element to grow from the background plasma turbulence.Such a precise functional form of Fo is by no means necessary; any distribution function giving a reasonably broad peak in the linear growth rate could be used in this model.With the above Fo the contribution to J res as a function of perpendicular velocity is usually peaked in the pitch angle range 40-60 degrees.As this code is band limited the new element cannot be produced until the previous one is terminated.This code thus best describes well separated, or isolated, chorus elements.In the case of tightly packed chorus the situation is as illustrated in Fig. 1; the next element starts to grow before the last one has ended.It is, however, well known (Nunn, 1986) that linear wave-particle interactions in parabolic inhomogeneous media exhibit an upper sideband instability and lower sideband damping, with the peak growth/damping rates at separations of the order of the trapping frequency (∼100 Hz here).We thus expect that the initial growth of the following element will be suppressed by this nonlinear quiet band effect until the frequency separation at a specific time and place is more than 2-3 trapping frequencies.This model will have an element separation which is the duration time of an element plus the time to grow the next element out of noise.These times turn out to be not unreasonable, but probably they are an overestimate.A broadband code would be needed to model fully the process whereby one element triggers the next.It should be noted again that, although chorus elements are often "tightly packed" with nearly constant element separations, the examples presented in Figs. 2 and 3 exhibit elements which are spaced rather widely and randomly.These are fairly consistent with this model.

Case 1; Ne=2/cc: low frequency chorus
We first present a numerical simulation of chorus elements which are similar to those observed by Cluster and shown in Fig. 2. The cold plasma density appropriate to this spectrogram is in the region of 2-4 electrons/cc.In this simulation we take Ne=2/cc, L=4.4 and an equatorial electron gyrofrequency of 10.3 kHz.The bandwidth of the simulation is 186 Hz.This simulation uses the relativistic version of the VHS code.In order to achieve a sufficient number of trapping oscillations the saturation amplitude needs to be quite large, in this case 200 pT, which is in line with the Cluster observations.The code uses a linear equatorial growth rate of 1200 dB/s, and large growth rates of this order are certainly needed with these parameters to maintain a self sustaining generation region.Indeed Santolik et al. (2003) have observed maximum temporal growth rates of 1560 dB/s.In fact the true linear growth rate will always be greater than the temporal growth rate, assuming that the wave amplitude increases downstream.This run employs a single value of energetic electron pitch angle, 45 degrees.
The triggering signal is broadband noise filling the bandwidth of the simulation, with an rms amplitude of 5 pT.This is modelled by 10 continuous wave (CW) signals with random phases and frequencies spaced equally across the simulation band.Figure 6 shows the spectrogram of two successive chorus elements, produced by overlapping block shaded Discrete Fourier Transforms (DFT) of the exit wave field data stream, emerging from the downstream end of the simulation box.The spectrogram data are in arbitrary units.We note that the successive triggering of rising frequency elements is completely repeatable and stable, and that the code will produce any number of elements.The sweep rate is not an exact quantity but was measured at 17.5 kHz/s, in very good agreement with the Cluster data spectrograms.The element separation time of 0.2 s is typical of this data set.If the second element were to start its growth at 0.1 s, where it would  Fig. 8. Spatial structure, as functions of the distance z from the geomagnetic equator, of a rising frequency chorus element generation region (GR) for the case Ne=3 electrons/cc.It is evident that the quasistatic wave amplitude envelope exhibits modulation due to resonant sidebands.The structure is sustained by the J r field (blue curve), and the J i field (green curve) continuously shortens the wavelength through the term d/dz(J i /R), which leads to an emission of rising frequency.A substantial contributor to the sweeping frequency is the gradient of wavenumber/frequency established across the GR.
be below the quiet band, the repetition interval would drop to about 0.1 s.
It should be noted that, due to approximations in the code, the turning point accelerations noted by Katoh and Omura (2007b) are not modelled here.Further extensive modifications of this code for highly relativistic cases are an important next step.

Case 2; Ne=3/cc
We continue to model the events shown in Figs. 2 and 3, using the relativistic code, though now assuming that Ne=3/cc.Again the electron gyrofrequency is 8.2 kHz, but the bandwidth is 110 Hz, the saturation amplitude 270 pT and the linear growth rate 1400 dB/s.A single perpendicular beam is employed with a pitch angle of 42 degrees; the resonant energy is now 110 keV.Rycroft (1976) illustrated graphically how the gyroresonance condition determines that the gyroresonant electron energy and the cold plasma density are related in an inverse manner as is evident here when cases 1, 2 and 3 are compared.
The spectrogram presented in Fig. 7 shows regularly repeatable riser elements produced with sweep rates ∼7 kHz/s.We note that amplitude modulation is present with a period ∼30 ms.This is also seen in Fig. 8 which gives a time "snapshot" of the situation at t=1342 ms, in the generation region (GR) of the riser.The top panel shows the wave amplitude profile.We see straight away that the GR is a quasistatic structure of finite size, some 5000 km long here, and located almost entirely on the downstream side of the equator.This explains why Cluster observes Poynting flux vectors directed away from the equator (Santolik et al., 2005).Amplitude modulation (i.e.sideband activity) is observed, with a spatial length of about 1000 km.The bottom panel shows the in phase current J r (blue line) and the out of phase current J i (green line) as functions of z.The currents are expressed as multiples of J unit which is the linear in phase resonant current appertaining to the saturation amplitude at the equator.All subsequent plots of current are in these units.The J r field maintains the amplitude profile in a roughly constant position Fig. 10.Wave amplitude history (in pT) of the simulation as a contour plot in the z, t plane, for the case of Ne=5 electrons/cc.The stable character of the GR is again evident.As power input falls with rising frequency, the structure slips downstream, away from the geomagnetic equator at z=0, at the group velocity which is in the upward direction.Fig. 11.History of the in phase current J r of the simulation as a contour plot in the z,t plane, for the case of Ne=5 electrons/cc.Since J r is normalised to unity for linear growth at the equator at the saturation amplitude and at the starting frequency, nonlinear growth rates peak at about 1.7 times the linear value.Note that the detailed structure of each element is by no means the same.while the J i field differentially rotates the wave phase in some average sense to give a shortening of the wavelength and thus a rising frequency signal.

Case 3; Ne=5/cc
We now consider the case of a cold plasma density Ne=5/cc and use the relativistic code.The parameters used for this run are a linear growth rate of 900 dB/s, a saturation amplitude of 250 pT, and an equatorial electron gyrofrequency of 8211 Hz.A single pitch angle beam, with a pitch angle of 31 degrees, is used here.The simulation bandwidth is 100 Hz and the resonant electron energy 81 keV at the start.
Figure 9 shows the spectrogram of simulated chorus with four clean and repeatable rising elements with a separa-48 Fig. 12. Fig. 12.For Ne=5 electrons/cc, the history of the out-of-phase current J i as a contour plot in the z, t plane for the whole simulation run.Note the large peak in J i of maximum value 2, giving a positive value to d/dz(J i /E) and hence positive frequency sweep rate.49 Fig. 13.Fig. 13.For Ne=5 electrons/cc, the plot of smoothed localised frequency in the z, t plane for the first element only.The development of a marked frequency/ wave number gradient across the GR is noted.Hence the advective term in df/dt is very significant.tion ∼400 ms.The frequency gradients were found to be 4450 Hz/s, again in accord with the measurements.Figures 10 to 13 give a concise overview of the entire history of the simulation.Figure 10 plots the wave amplitude as a function of position z and time.The closely quasistatic nature of the generation region is apparent here.When the frequency reaches >3200 Hz, the power input falls and the wave profile slips downstream at the group velocity.Figure 11 is the corresponding plot for the in phase current J r .This current is nearly always negative (indicating growth) and reaches maximum values ∼1.6 units.This implies nonlinear growth rates which are ∼1.6 times the linear values.Figure 12 is a plot of the J i field, where we note that J i is also generally negative, but somewhat larger than J r , with a maximum value ∼2.1.Figure 13 shows the local instantaneous frequency Fig. 14.Snapshot of the term in d/dz(J i /R) (Eq.20) in units of Hz/s during a rising element in for the case Ne=5 electrons/cc.The considerable variability is due to the sideband activity permitted within the simulation bandwidth.However, it is seen that the mean value is positive and hence this term is a significant contributor to df/dt.
(smoothed) plotted in the z, t plane; for reasons of simplicity, only the first element is shown.Subsequent elements reveal a very similar structure.The establishment of a significant frequency/wave number gradient across the GR will be noted.The advective term is thus a significant contributor to the sweep rate df/dt.To further clarify this, Fig. 14 gives a snapshot of the d/dz(J i /R) term in Eq. ( 20), in Hz/s as a function of z.There is some variability due to sideband waves but a positive mean exists ∼400 Hz/s, showing that this term is also a contributor to the sweep rate, but in this case one order of magnitude smaller than the advective term.

Case 4; Ne=20/cc
Increasing the cold plasma density further, Fig. 15 shows the spectrogram of simulated chorus for a cold plasma density of Ne=20/cc.This run and subsequent runs employ the non relativistic version of the code.This run has a considerably smaller linear growth rate of 300 dB/s and a smaller saturation of amplitude of 100 pT, which nonetheless gives a comparable degree of trapped particle nonlinearity, i.e. the same number of trapping oscillations in the equatorial zone.The chorus elements are stable, although varying somewhat, with sweep rates ∼4100 Hz/s.  on 13 September 2001 when the cold plasma density was ∼17/cc, on 21 October 2001 (at 20:10 UT) when the cold plasma density was ∼191/cc and on 22 December 2001 when the cold plasma density was ∼12/cc.We shall here simulate the lower band chorus for 21 October.
For this run we have Ne=191/cc.At higher plasma densities the linear growth rate and wave saturation amplitude may be set to lower values, at 14 pT and 220 dB/s, respectively, which gives a greater level of particle nonlinearity than for the risers and thus strong nonlinear growth rates.The   magnitude of the step in the distribution function is reduced by putting d=0.9 in Eq. ( 21). Figure 16 shows the computed spectrogram of the exit field.The falling elements are stable and very repeatable, with sweep rates ∼−1193 kHz/s.Each element shows an initial rise of some 100 Hz before the main falling frequency segment.Many numerical simulations all show this unobserved feature.The observed sweep rates given in Macusova et al. (2007) for lower band chorus at Ne=191/cc show a wide spread of negative sweep rates, from −1 to −6 kHz/s, with an average ∼−3 kHz/s.Our sweep rate is thus rather on the low side, but lies within the range observed.Interestingly the average sweep rate for upper band fallers is only ∼−1.2 kHz/s.
Because the wave-particle interaction time scales at Ne=191 /cc are much longer than in the earlier simulations, the element separations are inevitably longer, being about 3 s in this simulation.Macusova et al. (2007) report element separations up to 2 s, but averaging at about 0.8 s.Spectrograms of falling frequency chorus recorded on Cluster (Macusova, private communication) reveal tight packing of the elements, with each new element starting a quarter of the way through the present element, well before the termination of the previous element.A broadband code should be able to simulate this tight packing density.
Figure 17 shows the amplitude history of the first element of the simulation.It will be noted that the generation region extends some 4000 km upstream from the equator.Figure 18 plots the history, for the first element only, of in phase current J r , which becomes positive near the equator due to the rotation of phase trapping angle through the B-direction, causing a region of negative nonlinear growth and a marked dip in amplitude in this region.Figure 19 shows the corre-  sponding history of J i .It will be seen that J i becomes positive upstream of the equator in accordance with basic trapping theory, and it is this that causes the wave phase to be wound round the opposite way, and so giving a faller.The J i field sets up a positive gradient of frequency across the GR as shown in Fig. 20 which plots the smoothed localised frequency f (z, t).It is this gradient that is the dominant mechanism behind falling tones.It should be noted that the magnitudes of the resonant particle currents are somewhat greater by a factor ∼10 than in the case of the rising chorus.This is due to the higher level of nonlinearity and nonlinear growth rate, and also the weaker step discontinuity.
Since the falling frequency element generation region is postulated to extend further upstream, we now make a prediction which should be measurable in the data.Where chorus is bidirectional in the equatorial zone, compared to solely rising frequency chorus, we would expect a much wider region, ∼5000 km, where the net Poynting flux is indeterminate.The observations will be searched to see if this prediction is borne out in practice.

Case 6; Ne=191/cc: rising frequency chorus
In the high plasma density Ne=191/cc case of 21 October studied in Macusova et al. (2007), lower band rising frequency chorus is often observed.We now simulate this case.Since a weaker power input is required to generate risers we set the linear the growth rate at only 161 dB/s and the saturation amplitude at only 9 pT.The simulation is driven by a single pitch angle electron beam at a pitch angle 63 degrees, and the resonant energy is only about 4 keV.The bandwidth of the simulation is 64 Hz, compared with the trapping frequency of 40 Hz, thus allowing considerable   sideband activity.Figure 21 shows the simulation spectrogram.Repeatable and stable rising frequency elements are formed with sweep rates of up to ∼1.4 kHz/s.This compares with observed average sweep rates at this time of 1.87 kHz/s.The chorus element separation is somewhat variable, being up to 2.2 s, which is rather larger than the observed range of up to 1.3 s.As noted already, rising frequency chorus elements can strongly overlap so that this code will always take longer to generate the next element.Variability in spacing is also often noted in the Cluster data, and, from a theoretical viewpoint, this is not surprising at all.As the bandwidth is now greater than the trapping frequency, strong unstable upper sidebands appear as may be seen by the vertical structure of the spectrum along the frequency axis, which gives the amplitude modulation, or wavelets, as observed aboard Cluster.
Finally, for all the rising chorus simulations undertaken in this study, in Fig. 22 we plot sweep rates against values of the cold plasma density to investigate whether this functional dependence is in accord with that found by Macusova et al. (2007).Indeed we find that the sweep rate is a decreasing function of Ne, the shape and position of the curve being very similar to that illustrated by Macusova et al. (2007), within the limits of experimental error and the variability of the theoretical and simulation parameters.Note that all the rising chorus simulations are performed at roughly a constant level of particle nonlinearity, namely the saturation wave amplitude is decreased with increasing cold plasma density so that the approximate number of trapping oscillations experienced in the equatorial zone is the same.We could equally well state that the sweep rate is an increasing function of saturation amplitude at constant nonlinearity, a statement in accord with the sweep rate expression given in Omura et al. (2008),       distribution function, etc., do have an effect, but they are not so strong neither are they so readily quantified.Particularly at low densities, it is found that the predicted sweep rates do vary over quite a range, as indeed they are observed to do in practice.

Discussions and conclusions
We have considered ELF/VLF chorus signals observed aboard the four Cluster spacecraft at L∼4.4 in the equatorial generation region, and have simulated specific events occurring at frequencies less than half the electron gyrofrequency using a self consistent VHS code, inputting the observed plasma parameters, particularly the saturation amplitude, Lshell and the cold plasma density.The code has one spatial dimension (1-D) and assumes parallel propagation, which is in good accord with observations.In view of the large wave amplitudes observed, up to 300 pT, the code is fully nonlinear and allows trapping of electrons by the whistler mode wave in the equatorial zone.However, it considers all particles in the resonance region, including linear untrapped electrons.
First we considered the modelling of the generation region of a single rising frequency element of chorus, a "riser", resembling the situation for specimens of chorus where the elements are well scattered in time, with random spacings.We have successfully modelled rising frequency elements for a range of cold plasma densities from Ne=2 to 200/cc.The first comment is that, in every case, the wave amplitudes observed are sufficient to cause strong trapping of energetic electrons in the equatorial region.Secondly, the observed riser sweep rates are in excellent agreement with the observations, ∼+10 kHz/s at Ne=2/cc.However, the sweep rates of simulated falling frequency elements are a little lower than observed, being ∼−1 kHz/s instead of ∼−3 kHz/s at Ne=191/cc.In particular, the frequency sweep rates were found to decrease as the cold plasma density increases, just as observed (Macusova et al., 2007).The fact that the cold plasma density was the dominant determinant of sweep rate, but not the only one, was noted by Nunn et al. (2005) in their study of Siple triggered emissions.
For the simulations at L∼4.4 discussed here, the zero order electron distribution function Fo is not well known, so a multiple bi-Maxwellian model with an anisotropy of A∼2 is taken.The simulations are not sensitive to the choice of Fo, the main factors of interest being the pitch angle range providing the most power input (∼40 to 60 degrees here) and the linear equatorial growth rate.For each simulation there is a minimum value of the linear growth rate below which a self sustaining generation region is not possible.In each case presented here, we choose a value which is a factor ∼2 above this minimum.A summary of the parameters chosen to simulate different chorus elements is given in Table 1.It is evident that this chorus is simulated at frequencies somewhere between 2.4 and 3.8 kHz (i.e.0.29 and 0.46 of the electron gyrofrequency at the equator).As the plasma density increases, the time intervals between successive elements increase, as does the "dead" time, defined as the time interval when the program does not produce any emission.And, as Ne increases, the electron beam energy decreases, the linear equatorial growth rate decreases, and the saturation amplitude also decreases.
The fact that individual elements of chorus, especially for the predominant rising frequency tones, may be so closely modelled provides powerful confirmation that nonlinear cyclotron resonance (with trapping) is the underlying mechanism.Indeed, there is no other plausible plasma physical mechanism on offer.The simulations provide detailed information on the structure of the generation region (GR) along the geomagnetic field line.For a riser, the GR is of finite length, ∼several thousand km, embracing the trapping region and extending about 3000 km downstream from the equator.This is in good agreement with the Cluster observations that the chorus Poynting flux is directed away from the equator, except in a narrow (∼500 km) region which is very close to the equator.The simulations indicate that for a faller the GR extends some 2000 km across the equator.The theory therefore predicts a broader region about the equator for a faller than for a riser in which there is no net Poynting flux; this could be useful test of the theory.The 1-D theory has nothing to say about the lateral structure perpendicular to the ambient geomagnetic field, except to note that 1-D simulations have excess power input whose dissipation can only be due to a nonlinear spreading loss mechanism where the GR is a thin, field-aligned structure.
Table 1.Summary of input and output data for the six simulation runs.The dead time, the time between the end of one chorus element and the start of the next, is the duration of the quiet band plus the delay before the next elements starts.GR-left and GR-right are the approximate boundaries of the generation region relative to the equator.A highly significant result of these simulations is the repeatability and the stability of the production of the rising frequency elements.The sweep rate holds remarkably constant during each element and is repeated in the next element, the sweep rate being a function only of the ambient plasma parameters, particularly the cold plasma density.This opens up the possibility of producing successions of very similar elements.That only leaves the question of how the triggering of the next element is controlled.In this model the "start frequency" is assumed to be that of maximum instability, or rather less plausibly that of a strong triggering signal such as an MLR line.The "termination frequency" is that at which the linear growth rate falls to below the minimum required for a GR.
We have produced successions of elements here by adopting a zero order distribution function consisting of a multiple bi-Maxwellian with A∼2 superimposed on a smoothed out discontinuity in parallel velocity like a step in the distribution function.This produces a linear growth rate which is sharply peaked in frequency, although any distribution producing a peak in linear growth rate will suffice.In our simulations, each element grows from small amplitude random turbulence, the initial growth being linear.The exact mechanism whereby the next element is initiated is not properly addressed by this narrow band code.In practice we may expect a complex interaction between one element and the next, via the nonlinear wave-particle interaction process.In particular each element will exhibit lower sideband stability (a "quiet band") which is several trapping frequencies wide.This will delay the growth of the next element until a sufficient frequency separation is achieved.Furthermore, the wave-particle interaction process in one element may generate broadband "seed fields" for the growth of the next element.Where this is the case, we might expect an even spacing of elements.Where the seed field is ambient noise, a more random element separation might be expected.
Regarding further work, we note the very interesting simulations of Omura et al. (2008) which use a broadband PIC code to produce chorus emissions.Currently, due to the ex-treme computational demands of the code, the plasma parameters used are not representative of those actually observed.In principle, a broadband VHS code could be developed, but since Vlasov codes demand the proper resolution of distribution function throughout the phase space simulation box such a code would be extremely expensive to run.
Further in the future the strictly 3-D nature of the generation region needs to be addressed, and 3-D plasma simulations will be needed to gain a full understanding of the plasma physics of chorus.

Fig. 1 .
Fig. 1.Postulated mechanism for the retriggering of successive elements of chorus, shown in the frequency-time domain diagrammatically.Each element is developed from noise at the start frequency, e.g., corresponding to the frequency of maximum growth rate/instability.The element dies when the growth rate becomes too low.Spontaneous growth of the next element is delayed by the nonlinear "quiet band" effect.
Fig. 2 .Fig. 2. Spectrograms (i.e.dynamical spectra) of lower band chorus at frequencies below half the local electron gyrofrequency as indicated, as observed on all four Cluster satellites, measured by the Wide Band Detector (WBD) electric field instrument at an L-value=4.4,for six seconds at a Universal Time of 08:49:56 on 18 April 2002.MLat shows the magnetic latitude; these observations are made very close to the geomagnetic equator.
Fig. 3. Fig. 3. Spectrograms of chorus, in the same format as shown in Fig. 2, at 08:52:28 UT on 18 April 2002.The two black arrows are explained in the text.

Fig. 6 .Fig. 6 .
Fig. 6.Fig.6. Spectrogram of the exit wave field data stream for a cold plasma density of Ne=2 electrons/cc at L=4.4.The saturation amplitude is 200 pT, the initial linear equatorial growth rate 1200 dB/s, and the local electron gyrofrequency 10.3 kHz.The sweep rate of the chorus elements is 17.5 kHz/s.

Fig. 7 .Fig
Fig. 7. Spectrogram of the exit wave field data stream for a cold plasma density of Ne=3 electrons/cc.The saturation amplitude is 270 pT, the linear equatorial growth rate 1400 dB/s, and the local electron gyrofrequency 8.2 kHz.The sweep rate is ∼7 kHz/s. 45

Fig. 9 .
Fig. 9. Spectrogram of rising frequency chorus element for the case of a cold plasma density Ne=5 electrons/cc.The consistency, reproducibility and stability of the different elements are noteworthy, as is their constant time separation.

Fig. 15 .FigFig. 16 .
Fig. 15.Simulated chorus sequence spectrogram for a cold plasma density of N e =20 electrons/cc, with a sweep rate of 4.1 kHz/s.

Fig. 17 .
Fig. 17.Wave amplitude history (in pT) of the simulation as a contour plot in the z, t plane, for the case of falling chorus.The structure has high nonlinear growth rates and the wave amplitude extends 4000 km upstream from the equator.

Fig. 18 .
Fig. 18.History of the in phase current J r of the simulation as a contour plot in the z, t plane, for the case of falling chorus with Ne=191 electrons/cc.Note the low values near the equator due to the change of sign of inhomogeneity factor S. 55

Fig. 19 .
Fig.19.For falling chorus with Ne=191 electrons/cc, the history of the current J i for the whole simulation run.Note that J i becomes positive in the upstream sector where the inhomogeneity S>0.

Fig. 20 .
Fig. 20.For falling chorus with Ne=191 electrons/cc, the plot of smoothed localised frequency in the z, t plane for the first element only.The development of a marked negative frequency/wave number gradient across the GR is noted, which is the dominant cause of the frequency sweep. 57

Fig. 21 .
Fig. 21.Simulated lower band rising chorus sequence for a cold plasma density of Ne=191 electrons/cc.The linear growth rate is lower than the previous case at a value 160dB/s, which favours riser generation.There is considerable variability in element structure and duration.

Fig. 22 .
Fig. 22.The computer simulation sweep rates as a function of the cold plasma density.This figure presents the results from all the simulations carried out.The sweep rate of the chorus elements is a decreasing function of cold plasma density.Note that all these runs are at a constant degree of nonlinearity, and thus the saturation amplitude is itself a decreasing function of the cold plasma density.