The numerical simulation of the generation of lower-band VLF chorus using a quasi-broadband Vlasov Hybrid Simulation code

In this paper, we perform the numerical modelling of lower-band VLF chorus in the earth’s magnetosphere. Assuming parallel propagation the 1d3v code has one spatial dimension z along the ambient magnetic field, which has a parabolic z dependence about the equator. The method used is Vlasov Hybrid Simulation (VHS) also known in the literature as the method of Kinetic Phase Point Trajectories (Nunn in Computer Physics Comms 60:1–25, 1990, J Computational Phys 108(1):180–196, 1993; Kazeminezhad et al. in Phys Rev E67:026704, 2003). The method is straightforward and easy to program, and robust against distribution function filamentation. Importantly, VHS does not invoke unphysical smoothing of the distribution function. Previous versions of the VLF/VHS code had a narrow bandwidth ~ 100 Hz, which enabled simulation of a wide variety of discrete triggered emissions. The present quasi-broadband VHS code has a bandwidth of ~ 3000 Hz, which is far more realistic for the simulation of chorus in its entirety. Further, the quasi-broadband code does not require artificial saturation, and does not need to employ matched filtering to accommodate large spatial frequency gradients. The aim of this paper which has been achieved is to produce VLF chorus Vlasov simulations employing a systematic variety of triggering input signals, namely key down, single pulse, PLHR, and broadband hiss.


Introduction
VLF chorus is currently a topic of considerable scientific interest. It presents as a strong self-sustaining whistler mode emission, normally found at dusk/dawn outside the plasmapause. The generation region is near the equator of the earth's magnetosphere Li et al. 2013) and chorus is observed on satellites, for example Cluster , the Van Allen probes (Fu et al. 2014), the THEMIS mission (Bortnik et al, 2011;Li et al. 2013) or on the ground (Manninen et al. 2012;Macotela et al. 2019). Lower-band chorus waves have a main population that is quasi-parallel, but also a finite population of very oblique waves propagating close to the resonance cone angle (Agapitov et al. 2018;Li et al. 2013). Chorus normally has the form of a sequence of rising tone emissions, which may be narrow band and well-spaced or swishy and closely packed. Falling tone chorus is less frequently observed.
Numerous satellite observations suggest equatorial generation ) and parallel or quasiparallel propagation for lower-band chorus. Some satellite observations suggest though that upper band chorus has a wide angle between the k vector and local magnetic field. (Bortnik et al. 2007).
Chorus is usually divided into two bands, lower-band chorus below half the equatorial electron gyrofrequency, and upper band chorus above it Tsurutani and Smith 1974). In a recent detailed study (Teng et al. 2019), it was found that about 20% of events showed both lower and upper band chorus present with a distinct gap at half the gyrofrequency. Some 6% of events showed both bands present with no gap, and the remaining events were lower band only. Various theories have been advanced to explain this gap Ratcliffe & Watt 2017), but to date there is no totally convincing explanation. The gap is very precisely located and one would be looking for singular behaviour at half gyrofrequency in either propagation or wave-particle interaction dynamics.

The generation mechanism
Most historical theoretical and numerical work on chorus generation has assumed propagation parallel to the ambient magnetic field and would therefore be mainly applicable to lower-band chorus. Nearly all theoretical analyses assume that the VLF wavefield is narrow band, and it should be born in mind that the presence of sidebands or finite bandwidth will complicate the theory considerably. It is widely accepted that the generation process involves electron cyclotron resonance with the keV (up to MeV) hot electron population. The free energy required for generation of self-sustaining emissions derives from the hot electron distribution anisotropy, typically modelled as a bi-Maxwellian possibly with added loss cone. (Kennel & Petschek 1966).
A useful review of early research on nonlinear wave generation will be found in Omura et al. (1991). The complex spectral forms of chorus strongly imply that the wave-particle interaction process is fully nonlinear. The basic nonlinear process involved was expounded in Nunn (1974), and is nonlinear phase trapping of cyclotron resonant keV electrons in an inhomogeneous medium. A fully relativistic derivation of nonlinear resonant particle dynamics will be found in Omura et al. (2009).
Assuming parallel wave propagation, we define ψ as being the angle between the perpendicular velocity vector of a resonant electron and the perpendicular E field of the VLF wavefield, both vectors being perpendicular to the ambient field. It may then be shown that ψ obeys a pendulum equation with a constant force term S as follows: where z is the coordinate along the field line and S(z,t) is the so-called collective inhomogeneity factor, which has terms, in order of importance, in ambient magnetic field gradient, frequency sweep rate and gradient of cold plasma density: Here, B w is the local magnitude of the wave B field, and Bo(z) is the background magnetic field.
The quantity ω TR is the well-known trapping frequency or frequency of oscillation of electrons in the potential trap and is given by The solution to Eq. 1 is well known and has been often given in the literature (Nunn 1974). Full relativistic analysis will be found in Omura et al. (2009). Here, we will emphasise important aspects of the problem. If |S|> 1 inhomogeneity is too strong, and trapping will not be allowed. Electrons will sweep right through resonancethe so-called passing particles-and the problem will be far less nonlinear than the trapped particles and closely linear. This is because for passing particles interaction time with the wavefield is far less. For |S|< 1 electrons may be trapped-and they remain in resonance with the wave and relative phase ψ oscillates about the phase locking angle cos −1 S.
The conditions for such motion to be nonlinear are as follows. Not only must |S|< 1, but in the frame of the particle this must be satisfied for at least a trapping period. We also require that S itself is slowly varying on a time scale of the trapping period. If not then particle motion will be far more complex and not described by the above model.
Stable trapping in an inhomogeneous medium has profound implications for resonant particle distribution function, resonant particle current and thus nonlinear growth rates. For trapping times greater than a trapping period trapped electrons undergo large changes in invariants energy and magnetic moment, and they are dragged through invariant space. By application of Liouville's theorem we then expect to see a HOLE or a HILL in distribution function in phase space corresponding to the trap location as described in Omura et al. (2009). Now for rising frequency VLF emissions and rising chorus, the normal situation in the generation region is S < 0. Since the generation region is located at the equator where ambient field gradient is of order zero, then we expect S < 0 from the inhomogeneity due to the positive frequency sweep rate. We also require wave amplitude to be above the threshold for nonlinear growth ). With negative inhomogeneity, trapped particles undergo increases in energy and magnetic moment. Whether this gives a hole or a hill in distribution function depends of course on the shape of the unperturbed distribution function F 0 . For a positive linear cyclotron growth rate, it may be shown that S < 0 will give a hole in distribution function. The depth of the hole will depend on trapping time and the exact shape of F 0 . It must be realised that while some electrons remain trapped, others called passing particles stream right through resonance. Such particles undergo smaller changes in energy and magnetic moment, but there are far more passing particles than trapped and the total energy exchange with passing particles is of opposite sign to that of trapped particles and significant. This matter was examined in detail by Shklyar & Matsumoto (2009), who compared the contribution to nonlinear growth from trapped and passing particles.

The resonant particle currents J r (z,t) and J i (z,t)
The resonant particle current is evaluated from the familiar plasma theoretic expression which is the integral over velocity space of -ev times the distribution function F(v,r): Reference will be seen in the literature to 'particle bunching current' (Helliwell 1967). In phase space according to Liouville's theorem particles do not bunch, since particle density in phase space following a phase trajectory is conserved. This is true in the presence of a non-zero EM field. The concept of particle bunching as loosely defined in fact gives no indication of the result of the above integral and does not describe how the plasma current is produced.
The hill or hole in distribution function located at the particle trap in phase space will immediately give a sizeable nonlinear current, perpendicular to the ambient magnetic field, whose phase relative to the wave electric field will be controlled by the local inhomogeneity factor S (z, t) and hence will be either cos −1 S or cos −1 S + π. The component of current J r (z,t) parallel to the wave electric field will give wave growth for a linearly unstable zero order distribution function, and this may be several times the linear growth rate. The other component J i (z,t) parallel to the wave B field will modify wave phase and will be responsible for the sweeping frequency characteristic of individual chorus elements. These current fields will be controlled by the trapping dynamics in the inhomogeneous medium and will have a nonlinear dependence on the entire wave field and its history. The inhomogeneity and in particular the parabolic variation of B 0 (z) about the equator are absolutely key in understanding the chorus problem. Simulations in a homogeneous medium will give quite different results. In a homogeneous medium nonlinear trapping will give rise to growth rate saturation, as opposed to the enhanced growth rates of the inhomogeneous case.
Most theoretical analyses of particle trapping dynamics make the assumption that S is constant in the frame of a resonant particle. This will give a rather simplified view of the actual particle motion. In a parabolically varying zeroorder magnetic field S will have a roughly linear variation with z and thus time in the resonant particle frame. Particles first become trapped downstream from the equator at a point where S = − 1. Moving towards the equator S increases to S = 0 and the trap steadily increases in size, drawing in new particles. If a wavefield exists upstream from the equator S increases towards S = 1 moving away from the equator, and trap size decreases forcing once trapped particles out of resonance.
Another complicating factor when considering resonant particle motion comes from spectral broadening and sideband formation. It may be shown theoretically (Nunn 1986) that in a nonlinear trapping situation, assuming an unstable plasma, the wavefield is upper sideband unstable and lower sideband stable if S < 0, and vice versa if S > 0. The separation of resonant sidebands is of the order of the trapping frequency. Sidebands will resonate with the motion of trapped particles and drive them out of resonance. Where sidebands are of comparable amplitude to the main wave particle motion will become very complex and more chaotic in nature.

Numerical simulation of chorus generationprevious work
The paper by Golkowski et al. (2019) gives an excellent summary of active VLF experiments, and also compares various methodologies for the numerical simulation of VLF wave particle interactions. Previous simulations of VLF have used either PIC (particle in cell) methods or Vlasov methods. The PIC methodology is of course well known, widely used and a classical plasma simulation method. The essence of the method is as follows. The full set of Maxwell's equations are solved on a spatial grid. The linear cold electron population may be dealt with analytically or included as part of the hot population. The hot electron population is represented by a large number of weighted electrons distributed through phase space in accordance with the initial zero order distribution function. At each timestep the current appertaining to each simulation particle is distributed to the nearest spatial grid points using area weighting or higher order algorithm. This enables the spatial plasma current field to be computed and thus the forward push of the Maxwell set of equations.
The RISH group at Kyoto University and co-workers have developed over many years 1d3v PIC codes for the numerical simulation of VLF chorus generation (Katoh & Omura 2016;Hikishima et al. 2009Hikishima et al. , 2010aHikishima et al. , 2010b. These codes have one spatial dimension and assume the VLF wave field is parallel propagating, a reasonable approximation at least for lower-band chorus. The zeroorder magnetic field has a parabolic variation about the equator. The simulation is of nonlinear electron cyclotron resonance only as Landau and higher order resonances do not occur. The codes are fully relativistic, rigorously broadband, and cover the entire hot electron distribution function. With PIC codes it is not possible to limit the distribution function to a finite range in parallel velocity corresponding to a finite VLF frequency bandwidth. If this is done, at the point where the distribution is abruptly cut off, real plasma instabilities will be invoked thus invalidating the simulation. However, in the case of the ingenious and excellent PIC df algorithm (Sydora 2003;Tao et al. 2017), this will not apply and advanced df PIC codes should be able to consider finite ranges in parallel velocity.
PIC codes are in principle noisy as compared to Vlasov methods (Denavit 1972) and they require a very large number of simulation particles (~ 1 billion) in order to accurately represent the plasma dynamics. Successful runs were performed in which an initial input broadband noise developed into rising chorus waveforms consisting of a succession of rising elements. This represented the first chorus simulations using a fully broadband code. Although PIC methodology is computationally expensive PIC codes can be useful if one is investigating particle heating, diffusion and or particle precipitation, all very important issues (Hikishima et al. 2010b).
Thanks to use of greater computer power, the most recent 1D simulations with PIC code (Katoh and Omura 2016) are able to use magnetospherically correct parameters. The results are certainly the best so far with broadband noise input triggering a clean sequence of rising frequency chorus elements. As of now these PIC codes have been less successful in producing the less common falling chorus or hook shaped emissions. However, recently in realistic conditions falling tones have been simulated in a specific magnetic field gradient geometry. (Wu et al. 2020).
Successful simulations of chorus have been performed with the DAWN code (Tao et al. 2017). This code is PIC, but uses the dF method of Sydora 2003. The method is fully nonlinear, but particles are used to sample only the perturbed part of the distribution so it can significantly reduce the simulation noise. This method is very powerful and greatly reduces the noise level of PIC codes. The DAWN code has successfully simulated chorus with chirping signals.
Mention should be made of recent 2d3v PIC simulations of chorus generation in a dipole magnetic geometry, with two spatial dimensions (Ke et al. 2017). These simulations have been very successful at simulating rising chorus generated at the equator in quasi-parallel mode and which becomes increasingly oblique as it propagates away from the equator.
Simulations of discrete VLF emissions have been made using the VHS method, which is the subject of this paper (Nunn et al. 1997(Nunn et al. , 2009Nunn & Omura, 2012;Omura and Nunn 2011). The VHS codes used were with a limited bandwidth of the order of several trapping frequencies (~ 100 Hz bandwidth). A number of reasonable approximations are made, for example that group velocity remains constant over the bandwidth of the emission. The code can only simulate chorus if viewed as a succession of discrete emission events (Nunn et al. 2009). The central frequency of the simulation is continuously shifted as the emission frequency rises or falls. Risers are preferred but with higher growth rates fallers or hooks may be produced. Furthermore, because the solution for a chorus element with sweeping frequency involves substantial spatial dependence of frequency it was necessary to employ matched filtering to accommodate this. It was also necessary to employ artificial saturation to prevent wave amplitudes from achieving very large values.
The current code, which we may call quasi-broadband VHS has an order of magnitude increase of bandwidth to 3 kHz, and may dispense with matched filtering and artificial saturation. In this sense it is far more effective at simulating true chorus.

The Vlasov Hybrid Simulation method
Vlasov Hybrid Simulation (VHS) is a completely general method for the numerical simulation of collision-free plasma, where it is desired to resolve the distribution function in phase space. The method was first developed for the VLF-triggered emission problem (Nunn 1990(Nunn , 1993. The excellent paper by Kazeminezhad et al. (2003) proposed an identical algorithm, but named it the Method of phase point trajectories. The latter paper concerns electrostatic waves and covers various stability aspects not considered in the VLF papers. It is well worth reading.
The basic methodology of VHS is as follows: 1. Define the phase space simulation box which may be a function of time. Employ a phase space grid within the box. The grid may be adaptive and have varying resolution within the box. 2. Fill the box with simulation particles (SPs) with a density at least 1/grid cell. The density may be greater than this. Particles need not be located at grid points at the start or subsequently. 3. Follow particle trajectories which is the same as the phase space trajectories continuously forwards in time. Applying Liouville's theorem distribution function F l for the lth particle will be known and a constant value, being given by Fo at that particle's phase point at the start of the simulation. Simulation particles are embedded in the Vlasov fluid, and act as markers conveying information, namely the value of F at that point in phase space. 4. In order to compute accurately the plasma charge/ current field on the spatial grid, we need F defined on the regular phase space grid. It is therefore necessary to interpolate distribution function F from simulation particles onto the phase space grid. This process is fundamentally different from that of assignment of charge/current onto the nearest grid point, as occurs in PIC. The interpolation is easily achieved as follows.
Using the area weighting coefficients α k (as in PIC) the product α k F l is assigned to nearest grid points and the weighting coefficients themselves are also assigned. At each phase space grid point i the interpolant of F is then given by where for each grid point i the sum j is over all SPs that lie in the adjacent cells. The reader may be concerned that some grid points will have no adjacent SPs and thus no value of F. Applying Liouville's theorem to the assembly of particles which are the simulation particles it follows that density of simulation particles is conserved and particles do not bunch in phase space, in spite of references to 'phase bunching' in the literature (Helliwell, 1967). In this classic paper, a group of resonant particles were selected with the same energy and magnetic moment and equispaced gyrophase. It was found that as a result of a resonant wavefield the gyrophases of the particles bunched. It is not clear how this influences the plasma physical problem which is based upon Maxwell's equations and Liouville's theorem. The well-known hole/peak in distribution function does not arise from any bunching process but is understood by a direct application of Liouville's theorem. Particles are dragged through invariant space and retain the phase space density appropriate to their start point. To the author's knowledge no expression exists for a particle bunching current. By contrast the concept of phase trapping is well known and firmly based upon the particle equations of motion in an inhomogeneous medium. Trapped particles remain close to resonance for a long time, and the gyrophase relative to wave electric field oscillates about a phase locking angle dependent upon inhomogeneity. We may argue this point also as follows. Liouville's theorem concerns the interaction of a specified distribution of particles with a known wavefield. No consideration of self-consistency between particle charge/current and field is involved. For the total distribution function or subset thereof it states that density in phase space 'F' is conserved following a particle trajectory. Now there are two distribution functions involved here. One is the distribution function appropriate to the plasma physical problem to hand. The other is the distribution of simulation particles. This will normally be constant throughout the simulation box. Applying Liouville's theorem this distribution F will remain constant everywhere, except at the phase box boundary where Vlasov fluid is inflowing.
Particles move around in phase space maintaining a constant density. The VHS/VLF code in this paper continuously monitors occurrence of uncovered grid points and it is very low ~ 1/1000 with a particle density ~ 1.5 per cell. Increasing particle density will quickly reduce this figure to vanishingly small values. Where missed grid points occur a value of F can be found by interpolating from neighbouring grid points. One place where uncovered grid points occur is at the boundary of the phase box, where Vlasov fluid is flowing into the phase box. The VHS code must check boundary cells for absence of SP's at every timestep. Where this is the case a new SP must be inserted and assigned the zero-order value of F. The new particles must be inserted carefully to avoid excess particle density in the incoming Vlasov fluid. Placing them in the centre of the boundary cells was found to be satisfactory. Note that here there is a further restriction on simulation timestep. Where Vlasov fluid is incoming at the boundary a particle should not penetrate into the phase box by more than say 2 cells. The current code checks continuously for cells with no adjacent particles, and so if there is under density of SPs this will be flagged and the program halted. In many years of use the code has never halted in this way. Note that where Vlasov fluid has the value F = 0 it must still be represented by simulation particles, otherwise unresolved distribution function fine structure will give serious errors. In the VLF chorus problem in an inhomogeneous magnetic field there is a very strong flux of Vlasov fluid into the phase box as particles drift into resonance with the wave due to the inhomogeneity.
It is worthwhile to compare VHS with other algorithms. Although particles are followed continuously in time like PIC the algorithm is NOT the same as PIC as has been suggested. Application of Liouville's theorem to these phase space trajectories imparts far more information, namely the value of F, and allows the distribution function to be constructed. Indeed VHS resembles forward semi-Lagrangian more closely (Sonnendrucker et al. 1999;Besse & Sonnendrucker 2003), except the latter interpolates F back onto the phase space grid at each time-step and restarts trajectories at the grid points. This is not necessary, nor desirable as it results in nonphysical diffusion of the distribution function resulting from repeated interpolation onto the grid. VHS has the advantage that it is robust against distribution function filamentation, and can never produce negative distribution function. Distribution function fine structure is not resolved of course, but this has no effect on algorithm stability. VHS does not invoke unphysical diffusion of the distribution function, unlike conventional Vlasov methods (Cheng & Knorr 1976;Cheng 1977).
A detailed analysis of VHS algorithm stability will be found in Kazeminezhad et al, (2003). A separate 1D electrostatic VHS demonstrator code has been written at Southampton University, and it shows that energy and particles are conserved to a high degree of accuracy. Incidentally, in VHS codes particle density is not the number of simulation particles but the integral of F over the velocity space grid.
Excellent plasma simulations have been performed using the VHS method at North Western University in South Africa (Jenab & Brodin 2019). The important paper Killian et al. (2018) presents a highly optimised parallelised code for the VHS simulations of electrostatic waves and the reader is referred to this.

The quasi-broadband VHS VLF code
The quasi-broadband VHS/VLF code is a 1d3v simulation code with one spatial dimension for interaction between hot electrons (~ keV) and band-limited VLF waves, the assumption being that propagation is parallel to the ambient field direction. The ambient field is assumed to have a parabolic variation about the equator, as a function of coordinate z along the magnetic field. The field value at the equator and its spatial variation is appropriate to the L = 5 field line.
The nonlinear interaction takes place in the equatorial region. Well away from the equator, the large field gradient prohibits trapping and also linear cyclotron growth rates fall off quickly as cyclotron resonance velocity rises. The spatial simulation box extends about 7000 km either side of the equator and covers the region where nonlinear trapping is possible. Following each run it must be checked that the achieved wave amplitudes do not permit trapping outside the box.
The bandwidth of the simulation is about 3 kHz, centred on the base or central frequency, and the wave field is spatially band pass filtered to this band, by DFT/IDFT (discrete Fourier transform and its inverse). In response to changing wave spectrum, the band itself is moveable. The phase box dimension Vz is centred on the local resonance velocity at the current average frequency and has a width of at least 3 trapping widths, and must include the range of resonance velocities appropriate to the chosen simulation bandwidth. Of course if the average frequency changes, the range of Vz in the phase box will incrementally change. Only particles within the phase box, and with the prescribed range of Vz contribute to the nonlinear current.
The axes of the phase box are coordinates z, Vz, and gyro phase angle Ψ between the perpendicular velocity vector and wave electric field. The number of grid points along each axis are N z = 2048 in the z direction, N v = 500 in the Vz direction and N Ψ = 20 in gyro phase. The perpendicular velocity coordinate V ⊥ is relatively unimportant. The contribution to nonlinear growth comes from a relatively narrow range of perpendicular velocities centred on a pitch angle ~ 50 degrees. It is a reasonable approximation to consider each value of V ⊥ as a separate 1d2v VHS simulation, whence we may use quite a small number of V ⊥ values N V ⊥ = 10. Alternatively, one can have a larger value N V ⊥ = 30 say and run a 1d3v VHS simulation. With N V ⊥ = 10 the total simulation particle count comes to 300 M, assuming a density of 1.5 particles per phase space cell. Note that the relatively coarse resolution in V ⊥ is not permitted with PIC codes, which must use a full range of perpendicular velocities. However, where particle precipitation and heating are issues, PIC codes are very useful. This code is described as being quasi-broadband in that field push is achieved using the narrow band field equation appropriate to the triggering frequency rather than formally solving Maxwell's equations. This approximation is reasonable and avoids explicit computation of Maxwell's equations. Regarding the integration of the particle phase space trajectories, a simple leap frog scheme is used.
The simulation is initiated by introducing a triggering signal at the upstream boundary z = − 7000 km. The initiating signal may be a pulse or key down, and consist of either a CW single-frequency tone, or any number of frequencies. Broadband input is represented by a large number ~ 160 of single-frequency waves with small frequency separation and random phases. The broadband signal generated by the beginning and end of a pulse is removed by the filtering process in the code. In any case, the pulse ends are rounded off with an exponential shape function.
In this paper, we assume the zero-order distribution function Fo is multiple bi-Maxwellian with anisotropy factor A = 2.0 according to the usual definition as stated in Kennel and Petschek 1966. This is linearly unstable to VLF waves and provides the free energy to drive the nonlinear instability. The magnitude of resonant particle fluxes is fixed by specifying the linear electron cyclotron growth rate at the equator at the mean (base) frequency.

The broadband field equation
It is of some interest to spell out how the wave field is forward pushed in time through the broadband field equation. We assume parallel propagation coordinate z being the distance from the equator along the field line. There is no functional dependence upon coordinates x and y perpendicular to ambient field. We ignore small terms derived from duct curvature, duct inhomogeneity and ducting geometry. We first define wavefields perpendicular to the ambient magnetic field as and for the nonlinear resonant particle current and for the current due to the linear cold plasma From Maxwell's equations and the linear equation of motion of the cold plasma electrons, we may readily derive the differential equation governing the time development of the complex electric field: where Ω(z) is electron gyrofrequency and П(z) is electron plasma frequency. To develop a narrow band field equation, we now specify a base frequency ω normally the triggering or start frequency, but in the case of triggered risers could be selected at the riser half way point. We now divide out the fast base phase to give complex Nunn Earth, Planets and Space (2021) 73:222 electric field F and nonlinear current G that are slowly varying functions of z and t: It is now convenient to work in dimensionless units as follows: Define ω b = Ω e (0)/2 as half the equatorial electron gyrofrequency and k b = П e (0)/c then we de-dimensionalise with respect to t unit = 1/ ω b and z unit = 1/k b., v unit = ω b /k b . The dimensionless base frequency is given by ω 0 = ω/ω b and k 0 = k/k b . The dimensionless complex wave amplitude R and dimensionless complex current I are given by where ζ is the wave electric field phase relative to base phase. The dimensionless linear recursion relation now becomes The dimensionless group velocity V g is given by where in dimensionless units and From Eqs. 9-17, we may derive the narrow band dimensionless field equation at frequency ω 0 assuming amplitude R and current I are slowly varying with respect to z and t: g e = N e (z)/N e (0) = 1 + 0.5νz 2 . (18) The physical significance of the terms in the above equation are as follows. The complex wave amplitude is advected at the local group velocity and updated by the complex current I. Any errors in the constant term in curly brackets will have no great significance as this only corresponds to a change in resonant particle flux.
To develop the quasi-broadband field equation, we now allow complex wave amplitude R to have a bandwidth of ± 1500 Hz centred on the base frequency. The field is low-pass filtered to this band by DFT/IDFT, as Fourier components outside the band are not provided with resonant particles. We now ignore the z dependence of group velocity and use only the value of V g at the base frequency and at the equator. This means V g is a constant, which ignores dispersion. The quasi-broadband field equation now becomes Since V g is a universal constant, the complex wave amplitude may be defined on the field grid with grid spacing V g dt, where dt is the simulation timestep. At each step the complex field is advanced by one grid point, and is updated by complex I, cross interpolated from a separate grid on which particles and resonant particle current are defined. For the present runs timestep dt = 0.125 ms.
The broadband field equation is most accurate in the triggering phase when wave frequencies are close to the base frequency. Once chorus is triggered and wave frequencies are quite a bit different from the base frequency group velocity becomes inaccurate though the essential physics of the nonlinear wave particle interaction process is retained.
Since the wavefield is filtered by DFT/IDFT it is of course possible to perform the field push in the spatial Fourier domain as in Eq. 20, where subscript n denotes the nth Fourier component: In the actual code, the field is pushed using the time domain Eq. 19.

The resonant particle equations of motion
The basics of nonlinear wave-particle interaction between cyclotron resonant electrons and VLF waves has been extensively studied in the literature. In particular, the reader is referred to Katoh and Omura (2016) for a detailed study of the problem. Here, we summarise how the particle equations of motion are dealt with in this code. (20) Defining a complex dimensionless perpendicular velocity vector by we now define ψ as the gyrophase η relative to the base phase In the above dimensionless units, the equations of motion become where the dash symbol denotes d/dt. For perpendicular velocity we have The rate of change of gyrophase relative to base phase is given by where where V* is the difference between parallel velocity and resonance velocity at the base frequency and its time derivative given by where Q 0 (z) is factor due to the magnetic field inhomogeneity and cold plasma inhomogeneity: Since the zero-order distribution function F 0 (μ,W) is a function of magnetic moment μ and energy W, the change in these invariants is of special significance: The rate of change of particle energy W is easily shown to be (30) W ′ = −|R||V ⊥ | cos(ψ − ζ ) = 0.5ω 0 µ ′ , which gives the rate of change of magnetic moment to a very good approximation. In integrating simulation particle trajectories, which are phase space trajectories, the code follows the variables gyrophase ψ, V*, integrated energy change ΔW and hence integrated magnetic moment change Δμ. Also, the quantities |V ⊥ | and V z are required and need to be tracked. Applying Liouville's theorem the value of distribution function at the phase space point occupied by a particle is given by In the framework of the current formalism, we may recover the classical trapping equation. Define Θ as the phase between the perpendicular velocity vector and the electric field vector whence where the dash symbol represents d/dt in the frame of the resonant particle. The cyclotron resonance velocity V* res is given by the equation below and will be the centre point of the dimension V* in the phase space box: From Eqs. 27 and 33 we get the trapping equation where S is the collective inhomogeneity factor as discussed extensively in Omura et al 2009 andNunn 1990: where < ζ'' > represents a spatial average of the acceleration of additional phase viewed in the frame of the resonant particle. This quantity is related to frequency sweep rate through the approximate expression Here ω TR is the well-known trapping frequency of frequency of oscillation of electrons in the potential trap: Following the development of the chorus equations in Omura et al. 2009, we may use Eqs. 32-35 to derive an expression for frequency sweep rate. At the equator, which is the centre of the generating region, we have Q 0 (0) = 0 we may derive We assume strong nonlinear trapping in the generation region. Nonlinear growth maximises at S = − 0.4, so we assume this value at the equator. We have applied this equation to the riser in Fig. 1 assuming an rms wave amplitude at the equator of 50pT, which gives a sweep rate + 2709 Hz/s, quite close to the observed sweep rate. It would thus seem that the chorus equations are useful and valid and give a good estimator of sweep rate. It should be pointed out that Eq. 38 is not an exact science. Both trapping frequency and inhomogeneity are functions of pitch angle, and inhomogeneity S probably does not have to be exactly − 0.4 at the equator.

Results of the chorus simulations
The first run of the quasi-broadband VHS code will be of chorus triggered by weak hiss.
The data for this run are as follows. We assume L = 5 and take realistic values for ambient plasma values as follows. We take an equatorial cold plasma density of 5.4 electrons/cc., equatorial electron gyro frequency of 6700 Hz and electron plasma frequency of 20869 Hz. The variation of ambient field Bo(z) about the equator is parabolic and given in MKS units by Bo(z) = Bo(0) 1 + az 2 , where a = 9 10 -6 / (6370 L) 2 . The cold plasma density has a weaker parabolic z dependence The zero-order distribution function consists of two bi-Maxwellians, the first has T ⊥ = 44 keV and T = 15 keV, the second has T ⊥ = 192 keV and T = 60 keV. The distribution function is linearly unstable with anisotropy factor A ~ 2, and a linear equatorial growth rate at the start frequency of 130db/s. This growth rate corresponds to a ratio of hot electrons to cold electrons of Nh/Ne = 0.00521. The base frequency of the simulation is 2010 Hz.
The triggering signal is introduced at z = − 7000 km and is broadband covering the range 1500-2500 Hz. It is represented by 160 CW signals with equal separation and random phases. The signal is introduced continuously (key down) and is weak with an overall amplitude of 0.8pT.
The results of the simulation are presented in Fig. 1 as a frequency-time spectrogram of the complex amplitude sequence at the exit of the spatial simulation box at z = 7000 km. The spectrogram is derived from overlapping time DFTs with Hamming weighting. The spectral power plotted is in units of pT 2 Hz −0.5 . The same units are used throughout the paper. The start frequency of 2010 Hz is arbitrarily assigned and is below the frequency of maximum linear growth which is ~ 3 kHz. This gives better simulations than starting at the linear growth rate maximum. In reality chorus/ emission start frequencies may be determined as much by the availability of trigger signals as by the frequency of maximum growth. The frequency at which a rising emission stops is probably determined when the local linear growth rate falls below the threshold needed to produce a self-sustaining generation region (Nunn et al. 2009).
This simulation is of particular interest. The hiss band triggers a strong riser which then itself triggers a chorus sequence with realistic values of sweep rate of 2 kHz/s, element separation ~ 0.1 s. Note that the triggering hiss band does not figure in this plot as it is extremely weak. The achieved wave amplitudes are ~ 150pT. The quasi-broadband code does not require artificial saturation probably due to the increased fine structure in the wavefield which considerably reduces the extent of nonlinear trapping and nonlinear growth.
Ne(z) = Ne(0) 1 + 0.3az 2 . Fig. 1 Spectrogram of the exit wavefield at z = 7000 km. It shows a simulation of narrow band chorus, triggered by a strong riser, which in turn is triggered by a weak hiss band Figure 2 shows a short segment of the output wave field in pT observed at z = + 2000 km. It is seen that there is considerable fine structure with the output amplitude comprising a sequence of short wave packets of order 5-10 ms in length. Five recent papers (Nunn et al. 2021;Zhang et al. 2021Zhang et al. , 2020Zhang et al. , 2018Agapitov et al. 2013) comcern satellite observations of lower-band VLF chorus structure and compare with VHS simulation results and results from PIC codes. The probability distribution of packet lengths and frequency sweep rates of individual sub packets showed very good agreement.
Our second run has all the same plasma data as before, but the triggering signal is a key down narrow band signal at 2010 Hz and an amplitude of 7pT as in the well known Siple VLF experiment (Helliwell 1983). Figure 3 shows the frequency-time spectrogram of the data sequence at the exit point z = 7000 km. It is seen that three strong risers are produced with time spacings ~ 500 ms. This is the time taken to create a generation region starting from the trigger signal. These risers have peak amplitudes ~ 200pT and sweep rates of order 3 kHz/s. Interestingly, the first riser has a small downward hook at the top, and the second riser shows bifurcation in which secondary emissions are triggered by the main emission. All the emissions reveal sideband structure with rather variable sideband separation ~ 50-200 Hz. Figure 4 gives a time snapshot of the riser generation region, plotting wave amplitude, current Jr in phase with wave electric field and current Ji in phase with wave magnetic field, all as functions of z. The unit of current is the linear in phase current at the equator at the trigger frequency and trigger wave amplitude. The wave packet structure of the emission is clearly visible. The generation region is clearly located in the equatorial region, but the amplitude profile extends upstream from the equator by some 2000 km. The resonant particle currents are both negative, where Jr gives nonlinear growth and Ji winds up the phase giving sweeping frequency.
The third run has PLHR (Power line Harmonic radiation) as input trigger signal, modelled as five keydown CW signals each of 8pT, with 50 Hz separation extending from 1835 to 2135 Hz. All the remaining plasma data are the same. Figure 5 presents the f/t spectrogram at z = 7000 km. Clearly a substantial signal is triggered off the top of the original 5 lines. This has the appearance of rising chorus with separation ~ 10 ms.
The fourth run is a simulation of an emission with an N shaped frequency profile. Figure 6 shows a ground observation in Northern Finland of an unusual emission with an N spectrogram consisting of a riser, followed by a faller, then a riser. This entire emission then undergoes two 2 hop reflections which show successive dispersion and spectral broadening.
We have been able to simulate an N emission, but have not attempted to replicate the plasma data of the observed event as the L shell and cold plasma density are not known.
Our simulation of the N emission uses the same plasma data as previously except that the linear growth rate at the equator at the start frequency of 2 kHz is increased to 380db/s. The zero-order distribution function is a single bi-Maxwellian with T ⊥ = 62 keV and T = 20 keV or an anisotropy factor A ~ 2. The initiating signal is a key down doublet consisting of 2 waves of 80pT and a separation  Figure 7 shows the f/t spectrogram of the data sequence at z = 7000 km. The initial activity is complex with a self-triggered riser which triggers a slow faller off its mid-point. The faller then triggers a steep riser followed by a steep faller.
The main N event is triggered by the doublet at 2 kHz. This initial riser then becomes a faller. This is due to the very large growth rates at the turning frequency, causing the wave profile to creep upstream into the region occupied by a faller generation region. When the faller reaches the doublet frequency nonlinear trapping is compromised, and as a result of the reduced growth rate the wave profile slips downstream into the riser position, giving a riser. It might be noticed that the emissions go past the half-gyrofrequency, and there is no sign of a gap. All elements of the emission reveal substantial sideband structure, and appear to consist of steep falling elements stacked together.
The last run concerns development of spectral fine structure in hiss bands due to nonlinear wave particle interactions. It has been noted (Summers et al. 2014) that VLF hiss bands often exhibit fine structure. This run has plasma data the same as in run 1, except we have a lower linear growth rate of 90db/s. This is low enough Fig. 4 Time snapshot of the generation region of a strong riser from the previous simulation. The graphs show wave amplitude, resonant particle current Jr in phase with the wave electric field (blue curve) and current Ji in phase with the wave magnetic field (green curve). Note the strong sideband activity. The overall shape of the current curves is as expected for trapping in a parabolic inhomogeneity. Although the generation region is in the equatorial zone, the wave profile extends upstream from the equator by some 2000 km that full-scale triggering does not take place. The initial wavefield consists of a hiss band from 1500-2500 Hz represented by 80 CW waves with random phases. The RMS amplitude of the hiss band is 27pT and is large enough for nonlinear wpi effects to occur. Figure 8 shows the spectrogram of the exit field at z = 7000 km and shows substantial spectral structuring. This appears as stacked rising elements. This might resemble chorus but this is not a self-sustaining emission and should be regarded as spectral structure. Rudimentary risers appear at the band top, but there is insufficient growth rate for these to take off.

Conclusions
In this paper, we have simulated lower-band VLF chorus using the method of Vlasov Hybrid Simulation (VHS) also called the method of Kinetic Phase Point Trajectories (Kazeminezhad et al. 2003). Simulation particle trajectories are followed forwards in time continuously, which is synonymous with a phase space trajectory. From Liouville's theorem, distribution function F is then known at the phase points occupied by the particles. Interpolation of distribution function to the phase space grid is achieved with a simple low-order interpolator which does not permit negative distribution function. Higher order interpolators are not necessary as the interpolation process is only required to calculate the charge/current fields which are the integrals of distribution function over velocity space. Continuous forward integration in time is only allowable because in phase space particles do not bunch and leave grid points devoid of simulation particles. The algorithm is very simple and easily encoded. The trickiest part to program is where Vlasov fluid flows into the simulation box and particles need to be inserted into the Vlasov fluid with appropriate density. The VHS code has a finite bandwidth. The previous version had a fairly narrow bandwidth ~ 100 Hz or of the order of several trapping frequencies. The current version is quasi-broadband with a simulation bandwidth ~ 3 kHz, which is achieved by increasing the range of parallel velocity and thus N v . The wavefield is continuously spatially band pass filtered to the simulation band by DFT/IDFT, in order to remove out of band field which is not provided with resonant particles. This is not done to enforce algorithm stability. This is in contrast to PIC codes which model the entire distribution function and do not filter the wavefield. The code is termed quasibroadband because it uses the narrow band field equation to advance the wavefield.
The narrow band version of the code can only model discrete VLF emissions, and also required artificial saturation and use of matched filters to accommodate substantial spatial gradients of frequency and wave number.
By contrast the current quasi-broadband version is able to properly simulate chorus and complex triggering events involving considerable bandwidth as well as the situation where the trigger signal is broadband hiss. This code does not use matched filters and saturates naturally.
In contrast, the excellent 1d3v PIC codes of Omura and co-workers and Tao and the 2D pic codes of Ke et al. (2017) are all truly broadband and advance the wave field by solving Maxwell's equations directly. The characteristic time of field variables is now wave period and these codes have a much smaller timestep than VHS and are thus much more expensive. However, successful simulations of rising chorus have been achieved due to parallelisation of the codes.
Arguably, we have achieved our aim of simulating lower-band chorus with a Vlasov code. The quasibroadband version is far superior to the narrowband version and well suited to the chorus problem. Interestingly, the nature of the triggered event does depend on the trigger signal and to a large extent on linear growth rate/energetic particle flux. We have demonstrated chorus triggered by weak hiss, spectral structuring of strong hiss due to nonlinear wave particle interaction and triggering of a sequence of risers by a keydown CW signal.
Much work remains to be done in this field. A fully broadband VHS code needs to be developed that solves Maxwell's equations explicitly. This will be far more computationally intensive and require massively parallel code. Extension to 2D and even 3D codes is the next step and here Ke et al. (2017) have led the way. Vlasov codes are notoriously expensive at higher dimensionality but a 2d3v VHS code should at least be possible, particularly if the perpendicular velocity dimension is simplified. Eventually one hopes to see codes with three spatial dimensions. Modern PIC codes are very sophisticated and may be able to achieve this.

Funding
Partially supported by JSPS KAKENHI Grant JP17H06140. Supercomputing services on Iridis4 Computer provided by the University of Southampton.

Availability of data and materials
The single spectrogram of an N-shaped emission in Fig. 6 is obtainable from the author.