The nonlinear evolution of whistler-mode chorus: modulation instability as the source of tones

We review the modulation stability of parallel-propagating/field-aligned whistler-mode chorus (WMC) waves propagating in a warm plasma from a formal perspective with a focus on wave–particle interactions via ponderomotive forces. The modulation instability criteria are characterised by the group velocity dispersion, $d c_g/dk$, for whistler-mode waves and a condition on the ratio between the group velocity $c_g$ and the electron sound speed $c_{s,e}$. We also demonstrate that in order to investigate the spatiotemporal evolution of the envelope and the formation of packets (according to this mechanism), one necessarily needs to account for the motion of ions within the system, leading to an ionic influence on the modulation instability threshold determined by the ion fraction of the plasma. Finally, we demonstrate that chirping may be captured when higher-order effects are included within the spatiotemporal evolution of the amplitude. This yields not only an explicit expression for the sweep rate but also identifies a possible origin for the power band gap that occurs at half the electron gyrofrequency. Numerical validation demonstrates that the interaction between wave packets is a source for the emergence of tones observed within mission data, and such interactions may be a major source of the electron energisation which WMC are responsible for.


Introduction
Whistler-mode chorus (WMC) waves play a significant role in determining energetic electron dynamics within terrestrial and magnetospheric plasmas [Horne, 2005, Thorne et al., 2010, Artemyev et al., 2016, Woodfield et al., 2019].WMC are one particular manifestation of the so-called 'whistler-mode' electromagnetic/plasma wave [Stix, 1992], and are particularly noteworthy for their role in rapid electron energisation and pitchangle scattering [Bortnik et al., 2008, Omura et al., 2008a, Albert, 2010, Artemyev et al., 2018, Zhang, 2022].Discussions of the role of WMC as one driver among many (within the general context of energetic charged particle dynamics in the inner magnetosphere) can be found in e.g. Green & Kivelson [2004], Thorne [2010], Bortnik et al. [2016], Li & Hudson [2019], Lejosne et al. [2022].In perhaps overly simplistic terms, one can consider two main and contrasting challenges to achieving a full understanding of the role of WMC in magnetospheric plasma dynamics.While they are somewhat contrasting, both of these challenges are fundamentally united by the critical role that is played by wave-particle interactions [Brice, 1964, Kennel & Petschek, 1966, Tsurutani & Lakhina, 1997, Summers et al., 1998].
One challenge is to determine the direct impact of WMC on electrons that would otherwise evolve adiabatically as a geomagnetically trapped particle [Shklyar & Matsumoto, 2009, Albert et al., 2022].One of the most impressive manifestations of this approach ('wave effects on particles only') is the application of the resonant diffusion limit of the quasilinear theory (e.g.Kennel & Engelmann [1966], Summers [2005], Allanson et al. [2022]) to global-scale numerical modelling of the terrestrial and planetary radiation belt populations using Fokker-Planck radiation-belt models (e.g.Li [2014], Glauert et al. [2018], Wang et al. [2020], Allison et al. [2021].One of the numerous outstanding problems in this area is to understand and incorporate the role of so-called 'nonlinear wave-particle interactions' [Artemyev et al., 2021[Artemyev et al., , 2022]], with WMC playing a very important role. Another challenge is to instead try to solve for one or more of the generation, interaction and subsequent evolution of the WMC wave modes, as a function of e.g. a given initial plasma condition, and perhaps with some external driving or particle sources/injections.Studies of this nature try to understand the evolution of both the wave amplitude (i.e.amplitude amplification and modulation) and the structure in frequency space (i.e.either rising-or falling-tones, or even more exotic forms such as 'hooks').We should point out that the most general definition of WMC includes a variety of spectral forms, including comparatively structureless emissions (sometimes known as 'hiss-like chorus', e.g.see Tsurutani & Smith [1974], Li et al. [2012], Tsurutani et al. [2013], Gao et al. [2014], Shumko et al. [2018]), as well as the more well-known and coherent/structured/'chirping' rising-and falling-tones [Burtis & Helliwell, 1976, Koons & Roeder, 1990, Li et al., 2011, Li et al., 2012, Taubenschuss et al., 2014, Santolík et al., 2014, Teng et al., 2019].Approaches of this kind usually necessitate some form of 'self-consistent approach', in which one ultimately solves some variation of the Vlasov-Maxwell system [Schindler, 2007] given a number of constraints.Therefore one is likely solving first for the influence of unstable particle distributions on waves [Gary, 1993], and possibly also for subsequent resulting turbulence induced by wave-particle interactions and/or wave-wave interactions [Kadomtsev, 1965, Sagdeev & Galeev, 1969].
A number of important open questions remain regarding both of these challenges, and of course the separation into these two contrasting approaches (which for the purposes of this discussion have been crudely polarised as 'wave effect on particle only' and 'wave evolution only') is an approximation to the complex, dynamic and networked energy pathways of the inner magnetospheric plasma [Jaynes, 2015, Li & Hudson, 2019, Ripoll et al., 2020, Koskinen & Kilpua, 2022].
There have been a number of thorough recent reviews and discussions of WMC generation and evolution [Gołkowski et al., 2019, Tao et al., 2020, 2021, Zonca et al., 2021, Omura, 2021] and so we do not do a complete literature review, instead directing the reader to those references and therein.It suffices to say that the inhomogeneity of the background magnetic field is frequently invoked to play a key role in the chirping mechanism for WMC [Helliwell, 1967, Sudan & Ott, 1971, Nunn, 1974, Vomvoridis et al., 1982, Trakhtengerts, 1995, Omura et al., 2008a, Tao et al., 2021], facilitating resonant particle trapping, bunching, and the formation of 'resonant currents'.However there are some proposed mechanisms that do not rely upon the inhomogeneous background to drive the chirping behaviour [Zonca et al., 2021, Zonca et al., 2022], and in particular, we note recent particle-in-cell numerical experiments that demonstrate chirping behaviour within the context of a uniform background magnetic field [Wu et al., 2020].
The most significant contribution of this work is to demonstrate the key role that the ponderomotive force can play in driving chirping behaviour in WMC within the context of a homogeneous background field.This phenomena could now be considered in addition to other aforementioned mechanisms.
In this paper we simultaneously consider the role of wave-particle interactions both on the evolution of WMC and on the particle populations themselves due to ponderomotive forces.We will do so via the theory of modulations and weakly nonlinear theories, through which the evolution of the wave amplitude is coupled to the variations of number density, to observe how the interplay between the two manifests at the onset of nonlinear effects.Ultimately this will build upon the ideas introduced in previous theoretical treatments, and crucially those introduced by Omura et al. [2008b], Omura [2021] which investigated the evolution of WMC given a background particle population, now accounting for the simultaneous evolution of these species alongside the wave motion.However, we reiterate that the derivations presented in this work are limited to the case in which the background magnetic field is infinite and uniform.Therefore, in the context of e.g. the Earth's radiation belts, this implies that the analysis applies close to the geomagnetic equator.
The approach that we take is facilitated by several observations.The first is that coherent WMC waves are known to be narrow banded, with the bandwidth being approximately 10% of the local gyrofrequency [Santolík et al., 2003, Santolik et al., 2008], meaning that one may restrict the study of the dynamics to a single wave mode.Further, the asymptotic picture is simplified further by the observation that the majority of nonlinear generation processes and amplification events occur near-equatorially in the magnetosphere and confined to very limited latitudes [Lauben et al., 2002, LeDocq et al., 1998, Meredith et al., 2020], meaning it is not unreasonable to restrict ourselves to a fluid system in cartesian co-ordinates without latitude considerations and that we are within the remit of weak curvilinear magnetic field effects.Finally, by considering the field-aligned case the reduction procedure is much more straightforward owing to the fact that the Lorentz force v × B no longer contributes to the nonlinear effects.This is due to the fact that when the wave and the velocity field are oriented in the same direction, the Lorentz force does not generate anything beyond linear terms in an asymptotic theory.A consequence of this is that higher harmonics, that is non-zero integer frequency multiples of the carrier wave (e.g.±2ω, ±3ω, . ..), do not contribute to the wave motion.Instead, the wave evolution consists of simply the carrier wave interacting with the particle populations, and it is solely the wave-particle interaction via the ponderomotive forces that drives the nonlinearity observed in the wave's evolution.
With these simplifying factors considered, we are able to explore the emergence of wave packets in WMC using two main nonlinear approaches.The first of these is via classical modulation theory [Whitham, 2011], which postulates that the parameters of the wave such as the amplitude, frequency and wavenumber all evolve slowly over the course of many wave periods in a similar fashion to classical WKB theory.Such approaches have been successfully used within space plasmas [Mjølhus, 1976, Gribben & Parkes, 1977, Mjølhus & Wyller, 1986, Eliasson & Shukla, 2005, Tracy et al., 2014, Omura et al., 2008b, Omura, 2021], with the closest related work to this paper considering WMC to deduce the modulation stability of these waves [Tam, 1969], but with wave-particle interactions not fully accounted for and by use of a simplified version of of the dispersion relation.When a similar approach is also utilised for the electron beam velocity and number density with ponderomotive effects accounted for, we obtain our first insight into the ponderomotive-driven wave-particle interactions influencing WMC modulation.It is these interactions influencing the emergence of wave packets through electron-acoustic effects.Ultimately this builds upon the ideas introduced in theoretical treatments of the second approach, crucially those introduced by Omura et al. [2008b], Omura [2021] which investigated the evolution of WMC given a background particle population, and now account for the simultaneous evolution of these species alongside the wave motion.
The second approach is to undertake a formal multiple scales analysis to derive a direct evolution equation for the spatiotemporal evolution of WMC amplitude, taking the form of the Nonlinear Schrödinger (NLS) equation.This equation has emerged from heuristic arguments in prior works [Karpman & Washimi, 1977, Stenflo et al., 1986], with the key work of Krafft & Volokitin [2018] highlighting its emergence within a waveparticle interaction framework.However, one of the results of this paper is to demonstrate that a formal asymptotic procedure reveals that one should be cautious using such approaches due to an inconsistency that emerges in the induction equation, which is overlooked in previous approaches.As a consequence it highlights that the motion of ions, and not just a population that ensures neutrality, must necessarily be considered to resolve this inconsistency, and in doing so we find that these wave-particle interactions are augmented.Our paper demonstrates that although some of the qualitative conclusions of Krafft & Volokitin [2018] are the same, namely that the curvature of dispersion plays a role in the formation of solitons, the ion number fraction plays a critical role in whether WMC elements are to be observed.Therefore one of the main conclusions of this work is the statement that ion motion cannot be neglected in such problems.
With the understanding that arises from the above theoretical approaches, we are then able to augment these ideas within this paper by capturing WMC chirp, the mechanism behind rising and falling tones and one of the most intriguing features of WMC.This phenomenon presents itself as significant repetitive sweep of the dominant frequency peak of the wave.It is known that the NLS equation does not admit such behaviours and requires higher order effects to be introduced to account for this behaviour.This is precisely what we obtain as part of this work, extending the multiple scales analysis we demonstrate from a formal perspective that this process arises from the wave-particle interactions.A surprising consequence of this is that the extended model may provide an explanation for the observed band gap in WMC waves at half the gyrofrequency [Fu et al., 2014, Gao et al., 2019, Li et al., 2011, Chen et al., 2022] as the terms responsible for frequency variations vanish at precisely this frequency.As such, the theory suggests that sweep rates decrease for waves which approach the band gap before arresting completely.Overall, this analysis provides expressions for the sweep rate of a single WMC element, which demonstrates that within isolation a WMC element cannot produce monotonic tones with a net change in the frequency of the wave.Instead, it suggests that the chirping behaviour observed originates from the interaction of several WMC elements, which we corroborate with numerical experiments.These produce repetitive WMC tones, and the space-time series demonstrates that their interaction also generates pulsations in the wave envelope which manifest in the particle dynamics as amplifications in the energy density.This, we speculate, may shed light on which stage of their propagation WMC might be energising the electron populations they trap during transit.
The outline of this paper is as follows.We begin with a review of the modulation instability of parallel propagating WMC in §2, a process entirely driven by the ponderomotive wave-particle interactions, outlining where such waves become unstable and form subpacket structures.Subsequently, we derive evolution equations for the wave envelope of WMC in §3, leading to a Nonlinear Schrödinger (NLS) equation and revealing that ions play a crucial role in the envelope dynamics and alter the expected modulational stability transition.Owing to a lack of chirping behaviour, we add correction terms to the NLS that capture such effects in §4, leading to explicit expressions for chirping within a single element.Finally, we use numerical simulations of this model in §5 to demonstrate how the interactions between wave packets is the main driver of the chirping seen within the mission data.Concluding remarks are given in §6.

Review: Modulation Instability of Whistler Mode Waves in Electrononly Plasmas
To understand the formation of WMC wave packets, we must first analyse the necessary conditions that permit their formation.This is done from the viewpoint of modulation instability, the process under which uniform wavetrains destabilise and undergo amplitude modulations, ultimately to form several wave elements as wave energy clusters within packets [Ablowitz & Segur, 1981, Billingham & King, 2000, Whitham, 2011, Treumann & Baumjohann, 2001, Chen et al., 1984].The approach to identify this instability is to derive quasilinear modulation equations governing the slow evolution of wave, typically amplitude and wavenumber, and determine when this system possesses complex eigenvalues.Quasilinear modulation equations have been derived previously for the wave without particle interaction effects, with notable works relevant to our approach including [Tam, 1969, Omura et al., 2008a, Omura, 2021], but a key extension of this work will be to introduce modulation equations governing the electron number density and beam velocity.Such effects make a significant difference to the stability transition of the wave and thus it is pertinent to include such evolution simultaneously with the electromagnetic wave.Throughout this paper, we will be considering collisonless, warm, isothermal plasmas from a fluid description.In the first instance, we will be considering a non-relativistic plasma purely comprising of electrons, neglecting any ion influences for the moment but we note these will be accounted for in later sections.Thus, we will be concerned with the following equations of motion: In the above, B, E, V and n represents the magnetic field, electric field, electron velocity field and electron number density respectively.The parameters q = −e and m represent the charge and mass of an electron, µ 0 is the magnetic permittivity constant, c is the speed of light and c 2 s = k B T /m is the speed of sound for the purely electron plasma.The ponderomotive force F P acting on each electron is given by (see, for example, [Chen et al., 1984, Nicholson, 1983, Lamb & Morales, 1983]) where ω pe is the electron plasma frequency and ω is the frequency of the electromagnetic wave.In essence, this force determines the mean drift of electrons over rapid gyrofrequency oscillations due to amplitude modulations emerging within the monochromatic wavetrain.There are a number of choices one can make for this force depending on the plasma environment [Treumann & Baumjohann, 1997, Krafft & Volokitin, 2018], but in this body of work we consider the simplest such force corresponding to the lowest order ponderomotive effect.This is to illustrate that the presence of such a force, even in its most rudimentary form, is crucial for wave-particle interactions and as a consequence the formation of chorus wave packets and elements.
The framework in which this will be achieved is through the use of WKB theory, or equivalently Whitham modulation theory.The starting point for this will be to consider the following Stokes wave ansatz for a wave-particle solution, representing a parallel propagating right polarised wave in the presence of a uniform magnetic field with strength B 0 : where z is in the direction of ẑ.The parameters v || and n 0 represent the constant parallel velocity and reference electron number density respectively.The wave amplitude B W , mean velocity perturbation V and number density perturbation N are initially assumed to be constant but small, so that |B W | ≪ 1. Typically one is able to characterise this smallness by comparing linear and leading order nonlinear terms, which for WMC is achieved by comparing the convective term with the Lorentz force in the momentum equation: where ϕ is the wave normal angle, that is the angle between the wave and background magnetic field.It follows that this ordering of magnitude is trivially satisfied for any choice of the system parameters for parallel propagating WMC.Thus there is considerable freedom regarding the magnitude of waves this theory can consider, but must still be small enough to separate linear and nonlinear scales.The expansion procedure requires that N, V = O(|B W | 2 ), owing to the fact that these oscillation-free terms must balance the oscillation-free terms generated by the ponderomotive force.
The approach is to substitute this ansatz into the governing equations (1) and consider terms up to O(|B W | 3 ).The details of this calculation can be found within appendix A, but we summarise the key elements of the approach here.The carrier wave terms, when substituted into the governing equations, generate the dispersion relation which vanishes whenever the frequency ω = ω 0 (k, v || ) satisfies the typical Whistler-mode dispersion curve By continuing the analysis to higher powers of the amplitude to include amplitude-dependent weakly nonlinear effects, one finds the result ( which can be used to extract the nonlinear dispersion relations (6) These are denoted as nonlinear dispersion relations due to the presence of the wave amplitude, mean beam velocity and number density variation as corrections to the linear dispersion.The system (5) can be thought of as that obtained by Omura et al. [2008a], Omura [2021] but accounts for higher order nonlinear effects and inherent wave-particle interactions due to ponderomotive effects.In theory, the two approaches could be combined to augment the existing theory of the previous references, but this is not the focus of this work.Instead, to complete our analysis here, we demonstrate the above system can be cast in variational form, to make the subsequent analysis closer to classical wave modulation theory [Whitham, 2011].This is done by introducing with the factor nonvanishing for Whistler waves, allowing the system can be written as The above system of equations is generated by the B and N variations of Lagrangian density This Lagrangian can be thought of as that which is averaged over one period of the Whistler mode wave.
It is from this Lagrangian density that we will derive the conditions for the Whistler wave to undergo a modulational instability, associated with wave packet generation.Ultimately this will signpost this criteria in the electron case, which we will develop within the more applicable but involved ion-electron plasma case.

Modulation Instability
We may now study the Lagrangian (7) to determine when the parallel Whistler wave is expected to be stable and remain close to monochromatic.If it is not, it is expected to form packets (also known as elements) and as a result generate larger amplitude events.We do so by appealing to Whitham modulation theory [Whitham, 1970[Whitham, , 2011] ] which has in the past been utilised in plasma contexts for the same purpose (see [Mjølhus, 1976, Gribben & Parkes, 1977, Mjølhus & Wyller, 1986] for some key examples).Notably we will be augmenting the existing literature on plasma modulations to account for wave-particle interaction effects, requiring the consideration of an additional phase for the background velocity.The crux of the approach is to introduce the 'rapid' wave phase and velocity potential, This allows the wave parameters k, ω, v || and γ to vary slowly in space and time: where subscripts denote partial derivatives with respect to the subscripted variable.This admits the phase consistency conditions We can replace ω and γ using ( 6) to obtain the first two modulation equations The first of these equations is referred to as the conservation of waves, whereas the second is the classical conservation of momentum for the plasma in the presence of the wave.The subsequent two equations, which close the system, come from taking the Θ and Φ variations of L , which when simplified give These equations represent the conservation of wave action and the conservation of mass respectively.These 4 equations can then be written in quasilinear form [Whitham, 2011]: This quasilinear system of equations encapsulate the wave-particle interactions for a single wave -the second equation determines how the wave amplitude is modified by the presence of variations to the number density, whilst simultaneously the fourth equation describes how the local number density is altered due to the wave.The remaining equations, the first and third, dictate how the local frequencies and velocity field responds to the wave-particle interaction respectively and thus alter the properties of the wave propagation.These modulation equations for the wave parameters can now be analysed for their stability, which is achieved by investigating small perturbations to some fixed state.We take this constant state to be U 0 = (k 0 , B 0 , 0, N 0 ), noting that the choice of velocity does not meaningfully impact the results that follow.By considering perturbations of the form U = U 0 + δ Û e i(Z−CT ) , then the leading order perturbation is governed by the eigenvalue problem The resulting eigenvalues C ultimately determine the stability of this system -if they are all real the constant state U 0 is stable and the monochromatic wave perseveres, however when any of these eigenvalues is complex (occurring in complex conjugate pairs) there is an exponentially growing mode that causes the perturbation to rapidly diverge from the monochromatic wave state [Whitham, 2011].The emergence of these complex eigenvalues is more commonly referred to as a modulational instability.The characteristic polynomial for this problem admits 4 roots in general, which can be categorised by their values as B 0 N 0 → 0: We focus on the latter set of roots, as it transpires that the eigenvalues associated with the sound speed can be shown to be real but those associated with the group velocity can become complex.This is typical of problems involving waves coupled to a mean field, where it is the wave mode driving the instability [Bridges & Ratliff, 2022, Tam, 1970], and so is not unexpected here either.For this problem, these latter roots can be expanded in powers of B 0 and N 0 , again assumed small, to give where the effective nonlinear frequency correction ω 2 is given by .
It is clear that these eigenvalues are complex whenever ω ′′ 0 ω 2 < 0. To determine when this occurs more readily, we introduce the nondimensional scalings giving The only sign changes that happen within the Whistler wave interval 0 < W < 1 are the root of ω ′′ 0 , which for large α typical in the earth's radiation belts approaches W = 1 4 , a quarter of the gyrofrequency.The other sign change is due to the factor V 2 − ν 2 passing through zero, which occurs when the Whistler wave's group velocity goes from subsonic (c g < c s ) to supersonic (c g > c s ).Thus we have the following criterion for the modulation stability of Whistler waves: Modulational instability when This information is summarised in figure 1.It should be noted there are cases where |V | ∼ ν where the story is partially more complex and introduces a further stability boundary, but this is not generic and will not be discussed in detail here.As the speed of sound scales linearly with temperature within the setting considered, we can infer that the supersonic case is more prevalent in cold plasmas, with the subsonic case being expected in warmer plasmas.However, some caution should be noted here as the isothermal approximation is known to be poorer for warmer plasmas [Chen et al., 2012, Gao et al., 2014, Li, 2010], so a more technical theory to confidently conclude anything regarding this limit.
In summary, we have used classical modulation theory to explore the stability of monochromatic WMC waves with the additional consideration of the particle effects.Criterion for modulation instability, associated with the formation of wave packets, is deduced by exploring the nature of the eigenvalues of the 4 × 4 quasilinear system (8), revealing that both the group velocity dispersion ω ′′ 0 and the difference between the squares of the group and sound speed c 2 g − c 2 s .This identifies where one should look for the subelement structures typical of WMC waves, and the analysis of their evolution forms the remainder of this paper.To do so we will rely on the classical perturbative approach for the evolution of the wave envelope which we derive in the next section.

Ion Effects and Packet Generation
The classical Whitham modulation approach of the previous section grants us insight into the criteria necessary for WMC waves to develop into packets.It is the case, however, that to formally derive an equation for the spatiotemporal evolution of these packets one must also include the effects of ions within the analysis.It is not a priori obvious that this is necessary, especially as their role is greatly overshadowed by effects due to electron motion, but this analysis will highlight that ionic features, particularly the number fraction of ions present in the plasma, has a crucial role in the modulation stability of WMC.
With ion motions included, we will consider the two-fluid plasma description [Baumjohann & Treumann, 2012] with the assumption of isothermality for each particle species.This gives the system of equations ∇ × B = µ 0 (q e n e V e + q i n i s,e n e ∇n e = q e m e (E + V e × B) ) where the subscripts i, e denote fields and quantities which describe the ion and electron populations respectively.The ponderomotive forces for each species will be taken as The ponderomotive force for the ion equation is of the same form as the electron, which can be obtain by following the derivation in Chen et al. [1984].As in the purely electron case, alternate forms of this force may be supplied instead but this primary form of the ponderomotive force captures the essence of the wave-particle interaction.For the spatiotemporal analysis, we will undertake a formal weakly nonlinear analysis in the small parameter ε ≪ 1.This parameter is a characterisation of the wave amplitude, which in turn determines the strength of the nonlinear effects.Following typical amplitude equation approaches, we postulate the following expansions: where the slow variables Z and T , encoding two time scales, are defined as whose scalings are chosen so that the terms of the derived evolution equation are all of the same order in ε.
The smallness of the parameter ε is once again determined by the separation between linear and nonlinear scales, and following a similar calculation as that in §2, one finds that ε must satisfy the ordering which for field aligned WMC with a beam velocity parallel to the magnetic field is automatic, as sin ϕ = 0.Such scalings are typical within weakly nonlinear theories describing spatiotemporal amplitude evolution and ultimately indicate that the evolution equation for the amplitude is of Nonlinear Schrödinger type.Indeed, substitution of the above expansions into the governing equations and solving the resulting problems at each order of ε confirm this, as we will outline below.At leading order, O(ε), we find that ω, k satisfy the relation which has a root corresponding to the Whistler dispersion relation Waves along this dispersion branch require that The next order of the analysis generates both first harmonic and zero harmonic terms.The former of these may be solved to show that β 1 = 0 and The first of the zero harmonic (i.e.oscillation-free) problems arises from (12) and is simply This is to say that the parallel propagating wave does not induce an additional current at this order of the analysis.The second arises from (17) and gives that which is equivalent to quasineutrality holding in the presence of slow deviations.This allows one to write the electron number density variations according to where Z c is the ion charge number.Both of the conditions ( 18) and ( 19) highlight the importance of ionic effects, as without the ionic contribution this equation would necessarily yield that N e , V e must vanish and lead to no nonlinear effects emerging from the weakly nonlinear analysis.Thus this order makes it clear that ionic effects must be included in the study of the multiple scale WMC evolution and should not be neglected.The problem at O(ε 3 ) is where the analysis terminates, and only the first harmonic and zero harmonic terms need to be considered to develop the evolution equation that results here.The zero harmonic terms are more involved at this order, and these read where primes denote derivatives with respect to Z.The equation ( 18) to show that is we add q e times the first equation to q i times the second, we necessarily have that c g (q e N e + q i N i ) ′ = 0 , =⇒ q e N e + q i N i = constant.
From ( 19) we can see that it follows that this constant is zero.We will recast the system (20) by taking the third equation of ( 20) and subtracting qi qe = −Z c times the fourth.Overall, this gives that This system of equations can be inverted to show that the modulation of the mean particle quantities are related to the carrier wave by The arbitrary functions resulting from the integration have been ignored here, as these simply correspond to a shift in the frequency of the carrier wave which evolves much more slowly than the envelope.At O(ε 3 ), we require that all terms proportional to the first harmonic vanish, else the analysis will generate secular terms.It can be shown that in order to do so, the amplitude B W must satisfy the following Nonlinear Schrödinger (NLS) equation where the nonlinear frequency correction Γ is given by We conclude our formal derivation with a few marks about the validity and limitations of the above model.Primarily we expect the NLS equation above to be a good representation whenever the variations in the wave amplitude occur over scales much larger than the wave period, which is to say that the model is representative whenever the typical packet contains many waves.Observation from the Van Allen probes and THEMIS mission support this being typical of WMC [Zhang et al., 2019, Artemyev et al., 2022].It is also known that such envelope models are valid for evolution times up to O(ε −2 ), so one can expect much longer predictions for lower amplitude packets.However, nonlinear events in WMC typically happen on scales of less than a second [Santolik et al., 2008], with the most repetitive emissions taking place within windows of several seconds [Gao et al., 2022], suggesting that these events lie within the timespan of model validity.Finally a consequence of the scalings of the moving co-ordinate we have a narrow spectrum requirement of |k − δk|, |ω − δω| ∼ O(ε), where δk, δω represent the wavenumber and frequency associated with the amplitude B W , a constraint that mission data suggests WMC satisfies [Santolík et al., 2003, Santolik et al., 2008].

Influence of Ions on Modulation Stability
In order to analyse the modulation instability of Whistler mode waves, we will again nondimensionalise the wave quantities according to (9), and introduce the further nondimensionalisations Additionally, to simplify the analysis in the first instance, we will chose that v || = u || = 0. We note that the parameter r is typically small, with its largest value occurring for hydrogen ions where it takes the value r = 1 1836 .Thus we treat r ≪ 1.This simplifies the coefficients of the NLS (22) to It becomes clear that although the vast majority of effects due to ions characterised by r vanish, more or less recovering the modulation instability criterion of the electron-only plasma presented in §2.1.There is, however, a non-negligible ion effect that is crucial in the expression for ∆ associated with the subsonic to supersonic wave transition.The remaining expression involving the ion number density now controls one of the stability boundaries for the Whistler waves, as was demonstrated in §2.This expression indicates that the proportion of ions present in the plasma directly alter this boundary, and so the super-and subsonic regimes are altered to transition at ratios that can be much lower than unity.This is since Further, in the limit as η i → 0 where the electron fluid (i.e.static electrons) description is expected to be employed, one finds that the sound speed term is eliminated and only the supersonic case is operational, suggesting that parallel WMC would only rise in tone and parallel falling tones are expected to be nonexistent in this regime.This is in strong agreement with the mission data on parallel WMC waves, where rising tones are closely aligned with the magnetic field whereas falling tones are to be expected nearly orthogonal to it [Taubenschuss et al., 2014, Teng et al., 2019].We make this inference primarily based on the location of the region within which the wave is stable, however dynamically speaking we do not have a term which dictates the movement of the main spectral peak in frequency space.The nonlinear Schrodinger will need to be extended in order to accommodate such terms and test part of this hypothesis, which is precisely the aim of the next section of our analysis.

Spectral Asymmetry and the Emergence of Chirp
It is known that NLS models do not admit chirping behaviour for monochromatic waves or wave packets (see, for example, the standard textbooks Ablowitz & Segur [1981], Billingham & King [2000], Agrawal [2000] on the subject).This owes to the fact that the growth of sidebands due to modulation instability are symmetric and thus it has no preferred spectral shift.This limits the effect of the NLS solutions on the original carrier wave to simply a constant, time-independent frequency shift.To break this spectral symmetry, it is necessary to appeal to higher order effects within the weakly nonlinear theory that induce self-frequency shifting as has been explored in optics [Palacios et al., 1999, Goyal et al., 2011, Triki et al., 2016, 2022].We return our discussion to the approach outlined in §3.To obtain these higher order effects, it is necessary to analyse the contributions from the first harmonic terms at O(ε).The terms these generate are then included into the NLS equation ( 22) and treated as correction terms (as these are, strictly speaking, of a lower order than the original terms).The result of doing this for parallel propagating Whistlers is the modified NLS equation where the additional nonlinear correction has coefficient We now demonstrate how this additional term alters the growth of sidebands due to a modulation instability.This may be investigated by perturbing the Stokes wave solution, an exact monochromatic wavetrain solution to this amplitude equation, by for Stokes wave amplitude A, sideband wavenumber δk and P is a perturbation assumed small enough that quadratic terms in it are negligible.Upon substitution, it can be shown that the perturbation P admits the linear spectrum and thus there is a band of spectral wavenumbers κ which generates perturbations that grow in time whenever We can see that the modulation criterion ω ′′ 0 Γ < 0 may persist if the bracketed expression is positive.This is ensured even if δk > 0 as for parameter ranges operational in the magnetosphere one has εQ/Γ ≪ 1, and so δk must be significantly large in order to do so which would violate the assumptions used to derive (25).Thus the modulational instability has exactly the same thresholds as the NLS equation.The asymmetry in the wavenumber spectrum can then be determined by studying the maximum growth rate for the perturbed Stokes wave.This can be calculated as It is apparent now that the sideband wavenumber δk introduces bias to one side of the spectrum via this higher order term characterised by Q.This is determined by the sign of Q, which like Γ is determined by the factor ∆. Thus, terms with a lower sideband wavenumber will grow at a faster rate than their upper sideband counterpart.What happens subsequently cannot be inferred by this linear stability analysis, and requires further nonlinear approaches that we will not consider here.This observation is the first hint that this higher order term will be the driver of the chirping behaviour that we seek to understand.We will continue to explore this extended NLS equation from the perspective of nonlinear wave packets, the structures observed within WMC waves that ultimately generate the rising and falling tones we are interested in understanding.These structures will be derived and explored analytically and numerically in later parts of this section to explore how these may encourage a similar spectral asymmetry in WMC waves.

Emergence of the Band Gap
Of significant note within this weakly nonlinear theory, when accounting for the higher order terms, is that the classical band gap at half the gyrofrequency is a natural consequence.We can identify this by using the nondimensionalisations ( 9) and ( 24), where r ≪ 1, to simplify the coefficient of the term responsible for chirping to In this limit there is clearly a root of Q at W = 1 2 , namely ω = Ω e /2, which is at the typical band-gap for WMC waves.The only other sign change of Q is determined by ∆ (which ultimately leads to a breakdown of the weakly nonlinear theory as a whole) as the remaining factors result in a positive definite quantity in the lower band WMC range 0 < W ≤ 1 2 .This observation of this root of Q at W = 1 2 is an important one.It demonstrates that at half the gyrofrequency the amplitude evolution equations once again reduce to the NLS equation, so the wave mode neither biases a rise or fall in the frequency.Thus there is no emergent sweep rate for waves at this frequency, and it is clear that waves with frequencies close to this band gap will have very weak sweep rates too.Overall, this observation suggests that a WMC element that rise or fall in frequency will have its sweep rate continuously reduced as its spectral peak approaches this band gap until the wave becomes incoherent or its energy is exchanged with another wave of a different frequency.
Unfortunately, the narrow-band picture of the nonlinear theory is insufficient to fully test this hypothesis.For energy exchanges between different frequencies and how the wave responds, a multi-mode interaction system is required in the spirit of either a Coupled NLS/Manakov system [Manakov, 1974, Kourakis et al., 2005b,a, Baronio et al., 2012] or a Zakharov/Wave Turbulence model [Galtier et al., 2000, Newell & Rumpf, 2011, David & Galtier, 2022], both much more complex dynamical models.Therefore a full validation of this continual sweep rate reduction postulation is reserved for later study.We however note that the former of these extensions is essentially a coupled system of NLS-like equations and will retain many of the coefficients derived in this paper, and thus the argument regarding an arrested sweep rate in the proximity of the band gap we hypothesise to hold true.

Wave Packets with Chirp
Now that we have a modification to the NLS, as well as some insight into how chirp might come about for sidebands, let us formally derive a solution for a Chorus element that undergoes chirping.To do so, we follow Dias & Iooss [1993] and consider a solution of (25) of the form which is formed of three parts -the amplitude function R, the sideband wave component represented by the second term and a phase function ϕ, which will ultimately be the source of chirp.This is because the local frequency of this solution can be defined as the negative of time derivative of the total phase function, From this one may identify a sweep rate, defined as the rate of change of the frequency in time: Our aim to to relate this sweep rate with properties of the wave packet solution to determine the resulting frequency change, which is achieved though substitution of this ansatz into (25) and solving the resulting ODEs.
Starting with our guess at a solution, and splitting the resultant system into real and imaginary parts, we have the equations We start by taking the second equation, ( 28), and multiplying by R, which makes it an exact derivative.Integrating gives for some constant I.It will be conventient herein to introduce U = R 2 > 0, and note that for solitary wave packets one must have that U, U ′ and U ′′ tend to zero as ξ → ±∞.This yields that I = 0 and allows one to manipulate (29) to show that the sweep rate S is In modulationally unstable regions, this prefactor is negative, meaning that a positive sweep rate occurs when U ′ > 0, i.e. at points where the envelope of an element is increasing.We conclude the discussion here with an exact solution describing the envelope of a single Whistler packet/element.To do so, we note that we can set ν = 0 to simplify the analysis, and it can be reintroduced under suitable mappings.In this case, we get a quartic potential: A homoclinic connection ( that is a trajectory that begins and terminates at the same fixed point) within this system, corresponding to a solitary wave solution is only possible so long as there is an interval for which V (U ) > 0 .As B 4 > 0, we need B 2 < 0 and so Since σ is a free real-valued parameter, this is always possible to satisfy.Following [Kamchatnov et al., 2012] let us factorise this ODE as follows: where the roots are given by The negative subscript denotes the negative root and plus the positive root so that U − < 0 < U + .Thus, the expression in (31) under the square root is positive in the interval U − < U < U + , but as U = R 2 > 0 the only interval of interest to us will be 0 ≤ U < U + .The solution to this ODE is of the form of a Gardner/extended KdV soliton [Kamchatnov et al., 2012, Grimshaw et al., 2010]: and so the amplitude of this wave packet is A/B = U + .In summary, by considering higher order terms in the amplitude model we find that these terms are the source of chirping behaviour in both monochromatic wavetrains and within wave packets.In the case of the latter, however, we can observe that there is no net shift in the frequency due to the symmetry of the packet.Chirping due to WMC elements would appear to arise, therefore, upon packet interaction as this leads to asymmetry within the wave envelope and thus generates a net frequency shift.This is now explored using numerical means, where we will generate multiple packets, observe their interactions and extract their frequency spectrum to deduce the overall chirp produced by their interplay.

Numerical Simulation of Interacting Whistler Mode Elements
It is clear from the available spacecraft mission data that Whistler and Chorus wave packets are rarely isolated and propagate in groups of multiple packets.As a result, these packets interact and cause changes to their frequency and amplitude that the above insight of a solitary packet cannot provide.To investigate these interactions, we resort solving (25) using a timestepping procedure.For computational ease and to better identify the effects of system parameters on the sweep rate, we rescale (25), highlighting two cases dependant of the sign of ω ′′ 0 .We employ the scalings This reduces the problem to a single tunable parameter, which we can from the nondimensionalisations earlier that thus reducing the tunable parameters to just ε, W, Ω (or equivalently, B 0 ) and α which are determined by the plasma environment.For the majority of our simulations we will be using parameter choices representative of those found at L-shell L ∼ 6.This corresponds to α = 7.2 , B 0 = 1.4 × 10 −7 , n e = 1 × 10 7 , ε = O(10 −5 ).
In addition to this, we focus on wave frequencies lying in the lower frequency band range of 0.1Ω-0.5Ω.We advance (32) in time using an exponential time differencing method with fourth order Runge-Kutta timestepping (ETDRK4) [Cox & Matthews, 2002, Kassam & Trefethen, 2005], using periodic boundary conditions to take advantage of the speed and spectral accuracy of Fourier-based schemes.To initialise a multi-element solution and remain close to the analytic solution found, we initialise the simulations with the wave packet which has amplitude 1 A0−1 , width Λ and sideband wavenumber p.The width of these has to be suitably large (of the order 10) in order to generate multiple packets, with smaller Λ recovering a single wave packet with small amplitude wave radiation and larger widths generating many packets that partially fission before interaction.In our simulations, we allow this initial profile to evolve and we then observe the dynamics of this fissioned structure in time to determine how the emergent wave packets interact with each other and the resulting Fourier spectrum.
The result of simulating (32) yields wave packets that generate frequency sweeps, with examples appearing in figure 2. These emerge after the initial packet begins to split and interact with the new subelements, and develop further as the packets separate.This separation drives an envelope asymmetry, which from our observations in (30) would appear to the source for the frequency sweeps that emerge numerically.The sweep rate decreases over each simulation, as the subelements grow further apart, further reinforcing that it is the packet interactions that are the source of rising tones.We find over the course of our numerical investigations that higher magnitude of the ratio in (33) enhance the onset of frequency sweeping, which further points to the role of the modified term in the NLS plays in generating frequency shifts in the WMC waves.The structure of the envelope generated by the interaction of subelements is also noteworthy.In our simulations we find that the packets exhibit an almost time-periodic breathing behaviour, depicted in figure 3a) akin to the Kuznetsov-Ma soliton found in the NLS [Ma, 1979, Akhmediev et al., 1985], albeit with a zero background.The structure of maximum amplification for widths that generate 2-3 main packets also bears some semblance to higher order rogue wave solutions [Chabchoub et al., 2012, Slunyaev et al., 2013], a fact that is not unsurprising given that such solutions are the infinite-period limit of breathers.This seems to suggest, therefore, that isolated rising-tone WMC waves may fall into the category of rogue waves themselves, and the repetitive WMC emission events are breather events.This connection between hydrodynamic rogue waves and WMC has not yet been made in the literature and it may be useful to explore this connection further in the future, given also that the occurrence rate of large-amplitude WMC is significant in the magnetosphere.
We may also use ( 21) to explore the effect on both the number density and the electron beam velocity.From the latter, we may also extract the following energy density, a combination of kinetic and thermal energy of the electrons induced by the wave: (n e − n e,0 ), which allows us to determine which stage of the WMC evolution energises the electrons.We visualise these quantities in figure 3c), demonstrating that the maximum amplification of the wave is the stage which imparts the most energy to the particles, increasing the energy by several electronvolts.We may also correlate this moment with the frequency activity of the wave, as we do in figure 4. In it, we can observe that the amplification of the wave develops as one rising tone begins to terminate and another, lower in frequency, rising tone begins to initiates.The kinetic energy density peaks at exactly the point in time where the power of the higher frequency rising tone and the lower frequency rising tone are equal, and beyond this point the lower rising tone takes over as the most powerful part of the signal.This observation would therefore suggest that repetitive rising tone events, that are typical of chorus wave activity [Tsurutani & Smith, 1974, Li et al., 2012, Gao et al., 2014], are accelerating electrons most when the tones overlap and there is an exchange in wave power between these elements.It is difficult to conclusively identify this as the amplification mechanism however, as it could simply be that if the rising tone of higher frequency have been left to climb the amplification may have been much greater.Further study beyond that of this paper will be required to fully explore this interaction and its consequences.
We conclude this section with some final commentary on what aspects of WMC we cannot recreate with our current level of modelling, but could be captured with suitable alterations.Primarily, the MNLS equation appears to be unable to produce a single isolated rising tone event that has been successfully created via other theoretical modelling (for example, Nunn et al. [1997], Omura et al. [2008a], Tao [2014]).We attribute this to the fact that we have no source term corresponding to a hot electron population which ultimately drives the WMC waves in these existing theories.We anticipate that once such effects are accounted for properly, the simulations here should be able to reproduce these results.Further, another aspect of our simulations we are not able to reproduce is the non-sequential gap between events, as all WMC events in our simulations happen one after another.One can see from the mission data (cf.Agapitov et al. [2011], Gao et al. [2022]) that this need not be the case, and there can be considerable separation between WMC occurrences.We hypothesise that once hot electron growth effects are accounted for, these too may be emergent from simulations of our amplitude model.

Concluding Remarks
This paper has provided an overview of the formation and dynamics of parallel propagating WMC wave with rising tone.We have outlined the mechanism for packet formation is an instability of modulational type, whose transition is marked by a critical point in the group velocity as well as the sign of c 2 g − c 2 s .Further investigations demonstrated that the role of ions within the system is to reduce the latter instability threshold and suggests that their effects, although contributing little elsewhere, are non-negligible and must be accounted for in any analysis of WMC.
We have observed here that rising tones emerge from the NLS with higher order effects, and the degree at which the spectral peak changes can be attributed to either the strength of the background magnetic field and/or the frequency of the original carrier wave.We have also numerically investigated the evolution of multiple WMC wave packets, confirming that the interaction between WMC elements seems to be the mechanism for the emergence of their frequency sweeping.The envelope asymmetry appears to be the reason for this as this interplay overcomes the symmetry of the theoretical solitary WMC element solution which seemingly does not allow for.This suggested mechanism is supported by the data from the THEMIS and Van Allen missions, where the wave profile for rising tone elements contains both multiple elements and asymmetric wave envelopes.
The assumption of narrow-banded waves implicit in the derivation of the envelope equations captures the phenomenology of WMC.We do, however, note that formally the range of frequencies that the packets can sweep is limited by ε and should not cover a broad range of frequencies simultaneously.This means that the complete picture of a rising tone WMC will involve more complex models which couple multiple frequency bands together and transfer wave energy between one another.This has been seen in other systems, the simplest being coupled NLS models with the limit of such couplings being the Zakharov equation.For parallel propagating WMC, the theoretical analysis should remain tractable due to the lack of second harmonic terms in the perturbation analysis and yield further insight into the stability and evolution of nonlinear Chorus waves.
An important extension of the approach outlined here is to consider the more general case of WMC propagating obliquely to the magnetic field.This is more reflective of what is observed, where WMC is seen to most commonly propagate at angles of 10 • (close to parallel) or 70 • to the magnetic field.The work here remains reflective of the former of these cases by virtue of the obliqueness remaining small, but the latter case would require a revised version of the perturbative approach and result in a different modified NLS.
The mission data would also suggest that this version of the dynamics would instead lead to wave dynamics which admit falling tone WMC, a feature which may be linked to a negative version of the ratio (33).In either case, a quantification of the effect of oblique propagation should be determined and explored in a similar fashion to the field-aligned waves considered in this paper.
A Details of the WKB analysis of the electron plasma First order terms proportional to e iθ are given by the linear matrix problem The determinant of this matrix needs to be zero for nontrivial solutions and thus Useful for later are the following derivative results along the branch of solutions: For Whistler waves, one has that The left eigenvector of this matrix is also required in order to determine criterion for when can be solved for A 1 , namely that the Fredholme alternative is to hold and the right hand side vanishes when projected in the direction of the left eigenvector.This is given by The next order of the analysis only requires consideration of zero harmonic modes, associated to waveparticle interactions.This involves terms in which can be written in terms of gradients, and so by considering the expressions under the gradients we have the system Combining these gives This constant can be deduced by noticing that the wave-free flow is potential,as there is no indication that such waves have vorticity due to the constant underlying background magnetic field, so that This would mean that the overall version of this bulk equation would be At the final order we consider, we must now examine the terms proportional to e iθ , given by Now, we project the right hand side using l gives Now, we could impose that this is zero, but the analysis would be determined to be trivial.Instead, we assume that there is a correction to the frequency, δω, that has amplitude effects within it.Assuming ω = ω 0 + δω, the matrix system (34) has the amplitude dependent part δω Thus, the two equations that emerge from the WKB theory are

Figure 1 :
Figure 1: Signs of the right hand side term in (blue) and instability (red) for parallel propagating Whistler waves for subsonic (left) and supersonic (right) waves.Dashed line marks the asymptote ω = Ωe

Figure 3 :
Figure 3: Snapshots of the time series of a) the magnetic wave b) number density and c) kinetic energy density (right) at several spatial points, demonstrating the breather-like evolution of the wave envelope as the WMC wave travels.

Figure 4 :
Figure 4: Left: Snapshots of the magnetic field wave versus the power spectra for the time series snapshot for the parameter choice (ε, n e , α, Ω, W, k B T ) = (4 × 10 −5 , 1 × 10 7 , 7.2, 2.49 × 10 4 , 0.35, 7.373 KeV).Right: Comparison between the wave envelope (with maxima shifted to the same point in slow time t s ) over the time frame of one pulsation (the fast time t f ) and its short time Fourier transform.The white dashed line denotes the time at which the envelope achieves its maximum.