2. Review: modulation instability of whistler mode waves in electron-only 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 wave trains destabilise and undergo amplitude modulations, ultimately to form several wave elements as wave energy clusters within packets (Ablowitz & Segur Reference Ablowitz and Segur1981; Chen et al. Reference Chen1984; Billingham & King Reference Billingham and King2000; Treumann & Baumjohann Reference Treumann and Baumjohann2001; Whitham Reference Whitham2011). The approach to identify this instability is to derive quasilinear modulation equations governing the slow evolution of waves, 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 (Reference Tam1969), Omura et al. (Reference Omura, Katoh and Summers2008) and Omura (Reference Omura2021), but a key extension of this work will be to introduce modulation equations governing the electron number density and parallel 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. The non-relativistic assumption applies to the fluid velocities $v\ll c$ (with $v$ the magnitude of the velocity of a given fluid bulk motion), as is standard practice in such domains (Omura Reference Omura2021). Furthermore, the mathematics that follows in this paper that does consider individual particle/kinetic quantities (e.g. the gyrofrequency $\varOmega = |q|B_0/m$) does assume non-relativistic particle velocities. This is valid for the non-relativistic particle energies that we consider, for example, electron energies $\le$511 keV. Thus, we will be concerned with the following equations of motion:

(2.1)\begin{gather} \boldsymbol{\nabla} \times {\boldsymbol{E}} =-\frac{\partial {\boldsymbol{B}}}{\partial t}, \end{gather}

(2.2)\begin{gather}\boldsymbol{\nabla} \times {\boldsymbol{B}} = \mu_0 q n {\boldsymbol{V}}_e+\frac{1}{c^2}\frac{\partial {\boldsymbol{E}}}{\partial t}, \end{gather}

(2.3)\begin{gather}\frac{\partial {\boldsymbol{V}}_e}{\partial t}+({\boldsymbol{V}}_e\cdot \boldsymbol{\nabla}){\boldsymbol{V}}_e +\frac{c_s^2}{n}\boldsymbol{\nabla} n= \frac{q}{m}\left({\boldsymbol{E}}+{\boldsymbol{V}}_e \times {\boldsymbol{B}}\right)+{\boldsymbol F}_P, \end{gather}

(2.4)\begin{gather}\frac{\partial n}{\partial t}+\boldsymbol{\nabla} \cdot (n{\boldsymbol{V}}_e) = 0. \end{gather}

In the above, ${\boldsymbol B}$ and ${\boldsymbol E}$ represent the magnetic field and electric field, respectively. The electron velocity field and electron number density, ${\boldsymbol {V}}_e$ and $n$, respectively, are fluid quantities generated from the first two moments of a given distribution function (Chen et al. Reference Chen1984; Baumjohann & Treumann Reference Baumjohann and Treumann2012),

(2.5a,b)\begin{equation} n = \int f({\boldsymbol v}) \, {\rm d}{\boldsymbol v}, \quad n{\boldsymbol{V}}_e = \int {\boldsymbol v}f({\boldsymbol v}) \, {\rm d}{\boldsymbol v}. \end{equation}

The parameters $q = -e$ and $m$ represent the charge and mass of an electron, $\mu _0$ is the magnetic permittivity constant, $c$ is the speed of light and $c_s^2 = k_BT/m$ is the speed of sound for the purely electron plasma. The presence of this sound speed, and by virtue the $\boldsymbol {\nabla } n$ term in the momentum equation, are the result of the warm plasma and isothermal assumptions where $T$ is finite (Chen et al. Reference Chen1984). Such terms would be neglected within the cold plasma regime where $T\to 0$, whereas hot plasmas would be more accurately treated via a kinetic (Vlasov) approach rather than the fluid description of this paper (Chen et al. Reference Chen1984; Baumjohann & Treumann Reference Baumjohann and Treumann2012). The ponderomotive force ${\boldsymbol F}_P$ acting on each electron is given by (see, e.g. Lamb & Morales Reference Lamb and Morales1983; Nicholson Reference Nicholson1983; Chen et al. Reference Chen1984)

(2.6)\begin{equation} {\boldsymbol F}_P =-\frac{\omega_{pe}^2}{2\mu c^2 \omega^2}\boldsymbol{\nabla} \left( \langle |{\boldsymbol E}|^2\rangle\right), \end{equation}

where $\omega _{pe}$ is the electron plasma frequency and $\omega$ 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 wave train. There are a number of choices one can make for this force depending on the plasma environment (Treumann & Baumjohann Reference Treumann and Baumjohann1997; Krafft & Volokitin Reference Krafft and Volokitin2018), 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$:

(2.7)\begin{equation} \left. \begin{aligned} {\boldsymbol{B}} & = B_0\hat{\boldsymbol{z}}+(\hat{\boldsymbol{x}}-{\rm i}\hat{\boldsymbol{y}})B_W{\rm e}^{{\rm i}\theta},\\ {\boldsymbol{E}} & = \alpha_1(\hat{\boldsymbol{x}}-{\rm i}\hat{\boldsymbol{y}})B_W{\rm e}^{{\rm i}\theta},\\ {\boldsymbol{V}}_e & =v_{\|}\hat{\boldsymbol{z}} + \alpha_2(\hat{\boldsymbol{x}}-{\rm i}\hat{\boldsymbol{y}})B_W{\rm e}^{{\rm i}\theta}+{\rm c.c.}+V\hat{\boldsymbol{z}}\\ n & = n_0+N, \quad \theta = \omega t - k z, \end{aligned} \right\} \end{equation}

where $z$ is in the direction of $\hat {\boldsymbol {z}}$. The parameters $v_{\|}$ and $n_0$ represent the constant parallel velocity and reference electron number density, respectively.Footnote ¹ The constants $\alpha _1, \alpha _2$ characterising the wave modes of the electric and electron velocity, respectively, can be found as part of the linear theory and are explicitly given by

(2.8a,b)\begin{equation} \alpha_1 =-\frac{{\rm i} \omega}{k}, \quad \alpha_2 =-\frac{q}{mk}\frac{\omega-v_{\|}k}{\omega-v_{\|}k-\varOmega_e}, \end{equation}

as detailed in the Appendix (A), where $\varOmega _e = eB_0/m$ is the electron charge with $e$ being the elementary electron charge. 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|/B_0\ll 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,

(2.9)\begin{equation} \left\vert\frac{B_W}{B_0}\right\vert\sin \phi \ll \left \vert\frac{\omega -v_{\|}k}{\varOmega}\right \vert, \end{equation}

where $\phi$ 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 (2.1) and consider terms up to $O(|B_W|^3)$. The details of this calculation can be found within the 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

(2.10)\begin{equation} D(\omega, k,v_{\|}) = \frac{1}{c^2}\left[(\omega-v_{\|}k-\varOmega_e)(c^2k^2-\omega^2)+\omega_{pe}^2(\omega-v_{\|}k)\right], \end{equation}

which vanishes whenever the frequency $\omega = \omega _0(k,v_{\|})$ satisfies the typical whistler-mode dispersion curve

(2.11)\begin{equation} c^2k^2 = \omega_0^2-\frac{\omega_{pe}^2(\omega_0-v_{\|}k)}{\omega_0-v_{\|}k-\varOmega_e}. \end{equation}

The presence of the electron velocity within this dispersion relation is a natural consequence of the ansatz, rather than a particular assumption within our fluid approach, and acts like a Doppler shift on the plasma response portion of the dispersion curve. Typically this velocity is ignored within cold non-relativistic plasma fluid theory, which is permissible as $|v_{\|}|\ll c$, as it adds no new physics to the linear wave motion. It is known to have consequences on the nonlinear portion of the motion due to bulk effects (Krafft & Volokitin Reference Krafft and Volokitin2018), as we shall see within our analyses. We will make use of the fact we can neglect the particular value of $v_{\|}$ when evaluating later coefficients within the theory of this paper, as the choice of velocity does not meaningfully impact the results of this paper. The theory here may proceed so long as one avoids the choice that leads to the first cyclotron resonance,

(2.12)\begin{equation} v_{\|} = \frac{\omega_0-\varOmega_e}{k}, \end{equation}

where alternative treatments are required (Omura et al. Reference Omura, Katoh and Summers2008; Omura Reference Omura2021). By continuing the analysis to higher powers of the amplitude to include amplitude-dependent weakly nonlinear effects, one finds the result

(2.13)\begin{equation} \left. \begin{gathered} D B_W+\frac{\omega_{pe}^2}{c^2n_0}\frac{(\omega-v_{\|}k)^2-\varOmega_e \omega}{\omega-v_{\|}k-\varOmega_e}NB_W = 0,\\ \frac{c_s^2-v_{\|}^2}{n_0}N+\frac{\omega_{pe}^2}{\mu c^2 k^2} |B_W|^2 = \gamma -\frac{v_{\|}^2}{2}-c_s^2\ln n_0, \end{gathered} \right\} \end{equation}

which can be used to extract the nonlinear dispersion relations

(2.14)\begin{equation} \left. \begin{gathered} \omega = \omega_0(k,v_{\|}) +\frac{\omega_{pe}^2\left((\omega-v_{\|}k)^2-\varOmega_e \omega\right) }{\omega_{pe}^2\varOmega_e+2\omega(\omega-v_{\|}k-\varOmega_e)^2}\frac{N}{n_0},\\ \gamma = \frac{v_{\|}^2}{2}+c_s^2\ln n_0+\frac{c_s^2-v_{\|}^2}{n_0}N+\frac{\omega_{pe}^2}{\mu c^2 k^2} |B_W|^2. \end{gathered} \right\} \end{equation}

These are denoted as nonlinear dispersion relations due to the presence of the wave amplitude, mean velocity and number density variation as corrections to the linear dispersion. 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 Reference Whitham2011). This is done by introducing

(2.15)\begin{equation} \mathscr{B} = \frac{k^2}{n_0}\frac{(\omega-v_{\|}k)^2-\varOmega_e \omega}{\omega-v_{\|}k-\varOmega_e}|B_W|^2, \end{equation}

with the factor non-vanishing for whistler waves, allowing the system can be written as

(2.16)\begin{equation} \left. \begin{gathered} D +\frac{\omega_{pe}^2}{c^2n_0}\frac{(\omega-v_{\|}k)^2-\varOmega_e \omega}{\omega-v_{\|}k-\varOmega_e}N\equiv D+QN = 0,\\ \frac{c_s^2-v_{\|}^2}{n_0}N+Q\mathscr{B} = \gamma -\frac{v_{\|}^2}{2}-c_s^2\ln n_0. \end{gathered} \right\} \end{equation}

The above system of equations is generated by the $\mathscr {B}$ and $N$ variations of Lagrangian density

(2.17)\begin{equation} \left. \begin{gathered} \mathscr{L} = D\mathscr{B}+Q\mathscr{B}N+c_s^2N\ln n_0+\left(\frac{v_{\|}^2}{2}-\gamma\right)(N+n_0)+\frac{1}{2}\frac{c_s^2-v_{\|}^2}{n_0}N^2,\\ {\rm with} \quad Q = \frac{\omega_{pe}^2}{c^2n_0}\frac{(\omega-v_{\|}k)^2-\varOmega_e \omega}{\omega-v_{\|}k-\varOmega_e}. \end{gathered} \right\} \end{equation}

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 criterion in the electron case, which we will develop within the more applicable but involved ion–electron plasma case.

2.1. Modulation instability

We may now study the Lagrangian (2.17) 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 Reference Whitham1970, Reference Whitham2011) which has in the past been utilised in plasma contexts for the same purpose (see Mjølhus (Reference Mjølhus1976), Gribben & Parkes (Reference Gribben and Parkes1977) and Mjølhus & Wyller (Reference Mjølhus and Wyller1986) 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,

(2.18)\begin{equation} \theta = \varepsilon^{-1}\varTheta(Z,T), \quad \phi = \varepsilon^{-1} \varPhi(Z,T),\quad Z = \varepsilon z,\quad T = \varepsilon t,\quad \varepsilon \ll 1. \end{equation}

This allows the wave parameters $k, \omega, v_{\|}$ and $\gamma$ to vary slowly in space and time,

(2.19)\begin{equation} k =-\varTheta_Z, \quad \omega = \varTheta_T, \quad v_{\|} = \phi_Z, \quad \gamma =-\phi_T, \end{equation}

where subscripts denote partial derivatives with respect to the subscripted variable. This admits the phase consistency conditions

(2.20)\begin{equation} k_T+\omega_Z = 0, \quad (v_{\|})_T+\gamma_Z = 0. \end{equation}

We can replace $\omega$ and $\gamma$ using (2.14) to obtain the first two modulation equations

(2.21)\begin{equation} \left. \begin{gathered} k_T+\left(\omega_0(k,v_{\|}) +\frac{\omega_{pe}^2\left((\omega-v_{\|}k)^2-\varOmega_e \omega\right) }{\omega_{pe}^2\varOmega_e+2\omega(\omega-v_{\|}k-\varOmega_e)^2}\frac{N}{n_0}\right)_Z = 0,\\ (v_{\|})_T+\left(\frac{v_{\|}^2}{2}+c_s^2\ln n_0+\frac{c_s^2-v_{\|}^2}{n_0}N+\frac{\omega_{pe}^2}{\mu c^2 k^2} |B_W|^2\right)_Z = 0. \end{gathered} \right\} \end{equation}

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 $\varTheta$ and $\varPhi$ variations of $\mathscr {L}$, which when simplified give

(2.22)\begin{equation} \left. \begin{gathered} \mathscr{B}_T+\left(c_g \mathscr{B}-\frac{Q_k}{D_\omega}\mathscr{B}N\right)_Z = 0,\\ N_T+\left(v_{\|}(N+n_0)-\frac{v_{\|}}{n_0}N^2+Q_{v{\|}}N\mathscr{B}\right)_Z = 0. \end{gathered} \right\} \end{equation}

These equations represent the conservation of wave action and the conservation of mass, respectively. Thus, the complete modulation system is given by

(2.23a)\begin{gather} k_T+\left(\omega_0(k,v_{\|}) +\frac{\omega_{pe}^2\left((\omega-v_{\|}k)^2-\varOmega_e \omega\right) }{\omega_{pe}^2\varOmega_e+2\omega(\omega-v_{\|}k-\varOmega_e)^2}\frac{N}{n_0}\right)_Z = 0, \end{gather}

(2.23b)\begin{gather}\mathscr{B}_T+\left(c_g \mathscr{B}-\frac{Q_k}{D_\omega}\mathscr{B}N\right)_Z = 0, \end{gather}

(2.23c)\begin{gather}(v_{\|})_T+\left(\frac{v_{\|}^2}{2}+c_s^2\ln n_0+\frac{c_s^2-v_{\|}^2}{n_0}N+\frac{\omega_{pe}^2}{\mu c^2 k^2} |B_W|^2\right)_Z = 0, \end{gather}

(2.23d)\begin{gather}N_T+\left(v_{\|}(N+n_0)-\frac{v_{\|}}{n_0}N^2+Q_{v{\|}}N\mathscr{B}\right)_Z = 0. \end{gather}

Before proceeding, it is worth noting the connection between the system derived here and that of the system considered (Omura et al. Reference Omura, Katoh and Summers2008; Omura Reference Omura2021). When the electron fluid responses are neglected (so that $N, \gamma, v_{\|} = 0$) the system (2.23) reduces to a form that is identical to that of Omura when hot electron and magnetic curvature effects are neglected, consistent with the warm plasma and uniform field assumptions of this paper. Thus, the theory here could be used to extend parts of Omura's quasilinear wave analysis to account for higher-order nonlinear effects and the inherent wave–particle interactions due to ponderomotive effects, essentially taking the form of perturbations to the cold plasma distribution function. However, the aforementioned assumptions of this paper do not permit the formation of the phase space holes necessary to support the formation of rising tones within Omura's work. Instead, the formulation of this paper is a more local, fluid driven analysis focussing on a constant average electron velocity instead of a hot electron distribution. As a consequence, the mechanism which generates rising tones within this paper is of a different nature than that of Omura, but still rooted within wave–particle interactions, and is the subject of later sections within this paper.

The four equations in (2.23) can then be written in quasilinear form (Whitham Reference Whitham2011),

(2.24)\begin{align} \left. \begin{gathered} {\boldsymbol U}_T+{\boldsymbol A}({\boldsymbol U}){\boldsymbol U}_Z = 0, \\ {\rm with} \quad {\boldsymbol U} = \begin{pmatrix}k\\ \mathscr{B} \\ v_{\|}\\ N\end{pmatrix}, \quad {\boldsymbol A} = \begin{pmatrix} c_g-\left(\dfrac{Q}{D_\omega}\right)_k N & 0 & \dfrac{\partial \omega_0}{\partial v_{\|}}-\left(\frac{Q}{D_\omega}\right)_{v_{\|}}N & -\dfrac{Q}{D_\omega}\\ \omega_0''\mathscr{B} & c_g-\dfrac{Q_k}{D_\omega}N & (c_g)_{v_{\|}} & -\dfrac{Q_k}{D_\omega}\mathscr{B}\\ Q_k\mathscr{B} & Q & 0 & \dfrac{c_s^2}{n_0}\\ 0 & Q_{v_{\|}}N & n_0+N-\dfrac{N^2}{n_0} & Q_{v_{\|}}\mathscr{B} \end{pmatrix} \end{gathered} \right\}. \end{align}

This quasilinear system of equations encapsulates 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 ${\boldsymbol U}_0 = (k_0,\mathscr {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 ${\boldsymbol U}={\boldsymbol U}_0+\delta \hat {U} \textrm {e}^{\textrm {i}(Z-CT)}$, then the leading-order perturbation is governed by the eigenvalue problem

(2.25)\begin{equation} ({\boldsymbol A}-C{\boldsymbol I})\widehat{\boldsymbol U} = {\boldsymbol 0}. \end{equation}

The resulting eigenvalues $C$ ultimately determine the stability of this system – if they are all real the constant state ${\boldsymbol 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 Reference Whitham2011). The emergence of these complex eigenvalues is more commonly referred to as a modulational instability. The characteristic polynomial for this problem admits four roots in general, which can be categorised by their values as $\mathscr {B}_0\,N_0 \to 0$:

(2.26)\begin{equation} C ={\pm} c_s, \quad C = \frac{\partial\omega_0}{\partial k} \equiv c_g \ \mbox{(multiplicity 2)}. \end{equation}

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 Reference Bridges and Ratliff2022; Tam Reference Tam1970), and so is not unexpected here either. For this problem, these latter roots can be expanded in powers of $\mathscr {B}_0$ and $N_0$, again assumed small, to give

(2.27)\begin{equation} C = c_g\pm \sqrt{\omega_0''\omega_2\mathscr{B}_0}+O(\mathscr{B}_0,N), \end{equation}

where the effective nonlinear frequency correction $\omega _2$ is given by

(2.28)\begin{equation} \omega_2 = \frac{\omega_{pe}^4\omega^2( \varOmega_e-\omega)(\varOmega_e\omega_{pe}^2 + 2\varOmega_e\omega^2 - 2\omega^3)}{c^2(c_g^2-c_s^2)(2\varOmega_e^2\omega + \varOmega_e\omega_{pe}^2 - 4\varOmega_e\omega^2 + 2\omega^3)^2)}. \end{equation}

It is clear that these eigenvalues are complex whenever $\omega _0''\omega _2<0$. To determine when this occurs more readily, we introduce the non-dimensionalscalings

(2.29a–e)\begin{equation} \omega = \varOmega_e W, \quad k = \frac{\varOmega_e}{c}K, \quad c_g = cV, \quad \omega_{pe} = \alpha \varOmega_e, \quad c_s = c \nu, \end{equation}

giving

(2.30)\begin{equation} \left. \begin{gathered} \omega_0'' = \frac{2c^2\alpha^2(W - 1)^2(4W^4 - 4W^3 - 4W\alpha^2 + \alpha^2)}{\varOmega_e(2W^3 - 4W^2 + \alpha^2 + 2W)^3},\\ \omega_2 = \frac{\alpha^4\varOmega_e^4W^2(W - 1)(2W^3 - 2W^2 - \alpha^2)}{c^4(V^2-\nu^2)(2W^3 - 4W^2 + \alpha^2 + 2W)^2}. \end{gathered} \right\} \end{equation}

The only sign changes that happen within the whistler wave interval $0< W<1$ are the roots of $\omega _0''$,

(2.31)\begin{equation} \frac{W^3(W-1)}{W-\frac{1}{4}} = \alpha^2 , \end{equation}

which for large $\alpha$ typical in the earth's radiation belts approaches $W = \tfrac {1}{4}$, a quarter of the gyrofrequency. The other sign change is due to the factor $V^2-\nu ^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:

(2.32)\begin{align} \text{modulational instability when } \begin{cases} W^3(W-1)<\alpha^2\left(W-\frac{1}{4}\right), & |V|>\nu \quad ({\rm supersonic}),\\ W^3(W-1)>\alpha^2\left(W-\frac{1}{4}\right), & |V|<\nu \quad ({\rm subsonic}). \end{cases} \end{align}

This information is summarised in figure 1. It should be noted there are cases where $|V|\sim \nu$ 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 (Li et al. Reference Li, Thorne, Nishimura, Bortnik, Angelopoulos, McFadden, Larson, Bonnell, Le Contel and Roux2010; Chen et al. Reference Chen, Li, Bortnik and Thorne2012; Gao et al. Reference Gao, Li, Thorne, Bortnik, Angelopoulos, Lu, Tao and Wang2014), so a more technical theory to confidently conclude anything regarding this limit.

Figure 1. Signs of the right-hand side term in blue and instability in red for parallel propagating whistler waves for subsonic (a) and supersonic (b) waves. The dashed line marks the asymptote $\omega = {\varOmega _e}/{4}$.

In summary, we have used classical modulation theory to explore the stability of monochromatic WMC waves with the additional consideration of the particle effects. A criterion for modulation instability, associated with the formation of wave packets, is deduced by exploring the nature of the eigenvalues of the $4\times 4$ quasilinear system (2.24), revealing both the group velocity dispersion $\omega _0''$ and the difference between the squares of the group and sound speed $c_g^2-c_s^2$. 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.

3. 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 those ionic features, particularly the number fraction of ions present in the plasma, which have a crucial role in the modulation stability of WMC.

With ion motions included, we will consider the two-fluid plasma description (Baumjohann & Treumann Reference Baumjohann and Treumann2012) with the assumption of isothermality for each particle species. This gives the system of equations

(3.1)\begin{gather} \boldsymbol{\nabla} \times {\boldsymbol{E}} =-\frac{\partial {\boldsymbol{B}}}{\partial t}, \end{gather}

(3.2)\begin{gather}\boldsymbol{\nabla} \times {\boldsymbol{B}} = \mu_0 (q_e n_e{\boldsymbol{V}}_e+q_in_i {\boldsymbol{V}}_i) +\frac{1}{c^2}\frac{\partial {\boldsymbol{E}}}{\partial t}, \end{gather}

(3.3)\begin{gather}\frac{\partial {\boldsymbol{V}}_e}{\partial t}+({\boldsymbol{V}}_e \cdot \boldsymbol{\nabla}){\boldsymbol{V}}_e +\frac{c_{s,e}^2}{n_e}\boldsymbol{\nabla} n_e= \frac{q_e}{m_e}\left({\boldsymbol{E}}+{\boldsymbol{V}}_e \times {\boldsymbol{B}}\right)+{\boldsymbol F}_{P,e}, \end{gather}

(3.4)\begin{gather}\frac{\partial {\boldsymbol{V}}_i}{\partial t}+({\boldsymbol{V}}_i \cdot \boldsymbol{\nabla}){\boldsymbol{V}}_i +\frac{c_{s,i}^2}{n_i}\boldsymbol{\nabla} n_i= \frac{q_i}{m_i}\left({\boldsymbol{E}}+{\boldsymbol{V}}_i \times {\boldsymbol{B}}\right)+{\boldsymbol F}_{P,i}, \end{gather}

(3.5)\begin{gather}\frac{\partial n_e}{\partial t}+\boldsymbol{\nabla} \cdot (n_e{\boldsymbol{V}}_e) = 0, \end{gather}

(3.6)\begin{gather}\frac{\partial n_i}{\partial t}+\boldsymbol{\nabla} \cdot (n_i{\boldsymbol{V}}_i) = 0, \end{gather}

(3.7)\begin{gather}\boldsymbol{\nabla} \cdot {\boldsymbol E} = q_en_e+q_in_i, \end{gather}

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

(3.8)\begin{equation} {\boldsymbol F}_{P,j} =-\frac{\omega_{pj}^2}{2 \mu m_j n_j c^2 \omega^2}\boldsymbol{\nabla} \langle |{\boldsymbol E}|^2\rangle), \quad j = e, i. \end{equation}

The ponderomotive force for the ion equation is of the same form as the electron, which can be obtained by following the derivation in Chen et al. (Reference Chen1984). 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 $\varepsilon \ll 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:

(3.9)\begin{equation} \left. \begin{gathered} {\boldsymbol{B}} = B_0\hat{\boldsymbol{z}}+\varepsilon(\hat{\boldsymbol{x}}+{\rm i}\hat{\boldsymbol{y}})\left(B_W(Z,T)+\varepsilon\beta_1(B_W)_Z\right) {\rm e}^{{\rm i}\theta},\\ {\boldsymbol{E}} = \varepsilon(\hat{\boldsymbol{x}}+{\rm i}\hat{\boldsymbol{y}})\left(\alpha_1B_W(Z,T)+\varepsilon \beta_2(B_W)_Z\right){\rm e}^{{\rm i}\theta},\\ {\boldsymbol{V}}_e =v_{\|}\hat{\boldsymbol{z}} + \varepsilon(\hat{\boldsymbol{x}}+{\rm i}\hat{\boldsymbol{y}})\left(\alpha_2B_W(Z,T)+\varepsilon \beta_3(B_W)_Z\right){\rm e}^{{\rm i}\theta}+\varepsilon^2V_e(Z,T)\hat{\boldsymbol{z}},\\ {\boldsymbol{V}}_i =u_{\|}\hat{\boldsymbol{z}} + \varepsilon(\hat{\boldsymbol{x}}+{\rm i}\hat{\boldsymbol{y}})\left(\alpha_3B_W(Z,T)+\varepsilon \beta_4(B_W)_Z\right){\rm e}^{{\rm i}\theta}+\varepsilon^2V_i(Z,T)\hat{\boldsymbol{z}},\\ n_e = n_{e,0}+\varepsilon^2N_e(Z,T), \quad n_i = n_{i,0}+\varepsilon^2N_i(Z,T), \quad \theta = kz-\omega t, \end{gathered} \right\} \end{equation}

where the slow variables $Z$ and $T$, encoding two time scales, are defined as

(3.10)\begin{equation} Z = \varepsilon(z-c_gt), \quad T = \varepsilon^2 t, \end{equation}

whose scalings are chosen so that the terms of the derived evolution equation are all of the same order in $\varepsilon$. The smallness of the parameter $\varepsilon$ is once again determined by the separation between linear and nonlinear scales, and following a similar calculation as that in § 2, one finds that $\varepsilon$ must satisfy the ordering

(3.11)\begin{equation} \varepsilon \left\vert\frac{B_W}{B_0}\right\vert\sin \phi \ll \left \vert\frac{\omega -v_{\|}k}{\varOmega}\right \vert \end{equation}

which for field-aligned WMC with a velocity parallel to the magnetic field is automatic, as $\sin \phi =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 NLS type. Indeed, substitution of the above expansions into the governing equations and solving the resulting problems at each order of $\varepsilon$ confirm this, as we will outline below.

At leading order, $O(\varepsilon )$, we find that $\omega,k$ satisfy the relation

(3.12)\begin{align} D & = \frac{1}{c^2}\left[ (c^2k^2-\omega^2)(\omega-v_{\|}k-\varOmega_e)(\omega-u_{\|}k+\varOmega_i)\right.\nonumber\\ & \left.\quad +\omega_{pe}^2(\omega-v_{\|}k)(\omega-u_{\|}k+\varOmega_i)+\omega_{pi}^2(\omega-u_{\|}k)(\omega-v_{\|}k-\varOmega_e)\right] = 0 \end{align}

which has a root corresponding to the whistler dispersion relation

(3.13)\begin{equation} c^2k^2-\omega^2+\frac{\omega_{pe}^2(\omega-v_{\|}k)}{\omega-v_{\|} k-\varOmega_e} +\frac{\omega_{pi}^2(\omega-u_{\|}k)}{\omega-u_{\|} k+\varOmega_i} = 0. \end{equation}

Waves along this dispersion branch require that

(3.14)\begin{equation} (\alpha_1,\alpha_2, \alpha_3) = \left(\frac{ {\rm i} \omega}{k}, -\frac{q_e}{m_e k}\frac{\omega-v_{\|}k}{\omega-v_{\|}k-\varOmega_e}, -\frac{q_i}{m_i k}\frac{\omega-u_{\|}k}{\omega-u_{\|}k+\varOmega_i}\right). \end{equation}

The next order of the analysis generates both first-harmonic and zero-harmonic terms. The former of these may be solved to show that $\beta _1 = 0$ and

(3.15)\begin{align} & (\beta_2,\beta_3,\beta_4)\notag\\ & \quad = \left( \frac{c_gk-\omega}{k}, -\frac{{\rm i} q_e ((\omega-v_{\|}k)^2-\varOmega_e(\omega-c_g k) )}{m_e k^2(\omega-v_{\|}k - \varOmega_e)^2}, -\frac{{\rm i} q_i ((\omega-u_{\|}k)^2+\varOmega_i(\omega-c_g k) )}{m_i k^2(\omega-u_{\|}k + \varOmega_i)^2} \right). \end{align}

The first of the zero-harmonic (i.e. oscillation-free) problems arises from (3.2) and is simply

(3.16)\begin{equation} q_e(n_{e,0}V_e+N_ev_{\|})+q_i(n_{i,0}V_i+N_iu_{\|}) = 0. \end{equation}

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 (3.7) and gives that

(3.17)\begin{equation} q_eN_e+q_iN_i = 0, \end{equation}

which is equivalent to quasineutrality holding in the presence of slow deviations. This allows one to write the electron number density variations according to

(3.18)\begin{equation} N_e =-\frac{q_i}{q_e}N_i = \mathcal{Z}_cN_i, \end{equation}

where $\mathcal {Z}_c$ is the ion charge number. Both of the conditions (3.16) and (3.17) 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(\varepsilon ^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

(3.19)\begin{equation} \left. \begin{gathered} (v_{\|}-c_g)N_e'+n_{e,0}V_e' = 0,\\ (u_{\|}-c_g)N_i'+n_{i,0}V_i' = 0,\\ (v_{\|}-c_g)V_e'+\frac{c^2_{s,e}}{n_{e,0}}N_e' =- \frac{\omega_{pe}^2}{\mu m_en_{e,0} c^2 k^2}(|B_W|^2)',\\ (u_{\|}-c_g)V_i'+\frac{c^2_{s,i}}{n_{i,0}}N_i' =-\frac{\omega_{pi}^2}{\mu m_i n_{i,0} c^2 k^2}(|B_W|^2)', \end{gathered} \right\} \end{equation}

where primes denote derivatives with respect to $Z$. We may use equation (3.16) to show that, if we add $q_e$ times the first equation of (3.19) to $q_i$ times the second of (3.19), we necessarily have that

(3.20)\begin{equation} c_g(q_eN_e+q_iN_i)' = 0, \quad \implies \quad q_eN_e+q_iN_i = {\rm constant}. \end{equation}

From (3.17) we can see that it follows that this constant is zero. We will recast the system (3.19) by taking the third equation of (3.19) and subtracting ${q_i}/{q_e} = -\mathcal {Z}_c$ times the fourth. Overall, this gives that

(3.21)\begin{equation} \left. \begin{gathered} \mathcal{Z}_c(v_{\|}-c_g)N_i'+n_{e,0}V_e' = 0,\\ (u_{\|}-c_g)N_i'+n_{i,0}V_i' = 0,\\ (v_{\|}-c_g)V_e'+Z(u_{\|}-c_g)V_i'+\mathcal{Z}_c\left(\frac{c_{s,e}^2}{n_{e,0}}+\frac{c_{s,i}^2}{n_{i,0}}\right)N_i' =-\frac{Z\omega_{pi}^2m_en_{e,0}+\omega_{pe}^2m_in_{i,0}}{\mu m_em_i n_{i,0} n_{e,0} c^2 k^2}(|B_W|^2)'. \end{gathered} \right\} \end{equation}

This system of equations can be inverted to show that the modulation of the mean particle quantities are related to the carrier wave by

(3.22)\begin{equation} \left. \begin{gathered} \begin{pmatrix} N_i\\ V_e\\ V_i \end{pmatrix} = \frac{\mathcal{Z}_c\omega_{pi}^2m_en_{e,0}+\omega_{pe}^2m_in_{i,0}}{\mathcal{Z}_c\mu m_em_in_{e,0}n_{i,0} c^2 k^2\varDelta}(|B_W|^2) \begin{pmatrix} n_{e,0}n_{i,0}\\ \mathcal{Z}_c(c_g-v_{\|})n_{i,0}\\ (c_g-u_{\|})n_{e,0}\\ \end{pmatrix}, \\ {\rm with} \quad \varDelta = \left((c_g-u_{\|})^2-c_{s,i}^2\right)n_{e,0}+\left((c_g-v_{\|})^2-c_{s,e}^2\right)n_{i,0}. \end{gathered} \right\} \end{equation}

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(\varepsilon ^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 NLS equation:

(3.23)\begin{equation} {\rm i}B_T+\frac{\omega_0''}{2}B_{ZZ}-\varGamma |B|^2B = 0, \end{equation}

where the nonlinear frequency correction $\varGamma$ is given by

(3.24)\begin{align} \varGamma & = \frac{\mathcal{Z}_c\omega_{pi}^2m_en_{e,0}+\omega_{pe}^2m_in_{i,0}}{\mathcal{Z}_c\mu c^4m_em_in_{e,0}n_{i,0} k^2\Delta D_\omega}(\omega-v_{\|}k-\varOmega_e)(\omega-u_{\|}k+\varOmega_i)\nonumber\\ & \quad \times \left(\frac{\mathcal{Z}_c\omega_{pe}^2 n_{i,0}}{(\omega-v_{\|}k-\varOmega_e)^2}\left((\omega-v_{\|}k)^2-\varOmega_e (\omega-c_gk)\right)\right.\nonumber\\ & \left.\quad +\frac{\omega_{pi}^2 n_{e,0}}{(\omega-u_{\|}k+\varOmega_i)^2}\left((\omega-u_{\|}k)^2+\varOmega_i(\omega-c_gk)\right)\right). \end{align}

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. Reference Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky2019; Artemyev et al. Reference Artemyev, Mourenas, Zhang and Vainchtein2022). It is also known that such envelope models are valid for evolution times up to $O(\varepsilon ^{-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. Reference Santolik, Macusova, Titova, Kozelov, Gurnett, Pickett, Trakhtengerts and Demekhov2008), with the most repetitive emissions taking place within windows of several seconds (Gao et al. Reference Gao, Chen, Lu, Chen, Chen and Wang2022), suggesting that these events lie within the time span of model validity. Finally a consequence of the scalings of the moving coordinate we have a narrow spectrum requirement of $|k-\delta k|, |\omega - \delta \omega | \sim O(\varepsilon )$, where $\delta k, \delta \omega$ represent the wavenumber and frequency associated with the amplitude $B_W$, a constraint that mission data suggests WMC satisfies (Santolík et al. Reference Santolík, Gurnett, Pickett, Parrot and Cornilleau-Wehrlin2003; Santolik et al. Reference Santolik, Macusova, Titova, Kozelov, Gurnett, Pickett, Trakhtengerts and Demekhov2008).

3.1. Influence of ions on modulation stability

In order to analyse the modulation instability of whistler mode waves, we will again non-dimensionalise the wave quantities according to (2.29a–e), and introduce the further non-dimensionalisations

(3.25)\begin{equation} r = \frac{m_e}{m_i}, \quad c_{s,e} = c \nu. \end{equation}

Additionally, to simplify the analysis in the first instance, we will choose 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 = \tfrac {1}{1836}$. Thus, we treat $r \ll 1$. This simplifies the coefficients of the NLS (3.23) to

(3.26)\begin{equation} \left. \begin{gathered} \omega_0'' = \frac{2\alpha^2c^2}{\varOmega}\frac{(W - 1)^2(4W^4 - 4W^3 - 4W\alpha^2 + \alpha^2)}{(2W^3 - 4W^2 + \alpha^2 + 2W)^3}+O(r),\\ \varGamma = \frac{\alpha^4 n_{i,0} \varOmega}{\mu m_en_{e,0} \varDelta}\frac{(W - 1)^2(2W^3 - 2W^2 - \alpha^2)}{ (2W^3 - 4W^2 + \alpha^2 + 2W)^2(W^2 - \alpha^2 - W)}+O(r),\\ {\rm with} \quad \varDelta = c^2(n_{e,0}+n_{i,0})\left[V^2-\frac{n_{i,0}}{n_{e,0}+n_{i,0}}\nu^2\right]+O(r). \end{gathered} \right\} \end{equation}

In addition, our perturbations to the electron number $N_i$ and their velocity characterised by $\alpha _4, V_i$ become, in this limit,

(3.27)\begin{equation} \left. \begin{gathered} N_i = \frac{\alpha^2 n_{i,0}}{\mathcal{Z}_c\mu m_e c^2 k^2\varDelta}(|B_W|^2)+O(r),\\ V_i = \frac{\alpha^2 V }{\mathcal{Z}_c\mu m_e n_{e,0} c k^2\varDelta}(|B_W|^2)+O(r),\\ \alpha_4 = O(r). \end{gathered} \right\} \end{equation}

It becomes clear that although the vast majority of effects due to ions characterised by $r$ vanish, save for the zero harmonic/bulk variations, 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 $\varDelta$ 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 supersonic and subsonic regimes are altered to transition at ratios that can be much lower than unity; this is since

(3.28)\begin{equation} 0\leq \frac{n_{i,0}}{n_{e,0}+n_{i,0}} \equiv \eta_i \leq 1. \end{equation}

Further, in the limit as $\eta _i\to 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 non-existent 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. Reference Taubenschuss, Khotyaintsev, Santolík, Vaivads, Cully, Contel and Angelopoulos2014; Teng et al. Reference Teng, Tao and Li2019). 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 NLS 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.

4. 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 (Reference Ablowitz and Segur1981), Agrawal (Reference Agrawal2000) and Billingham & King (Reference Billingham and King2000) 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. Reference Palacios, Guinea, Fernandez-Diaz and Crespo1999; Goyal et al. Reference Goyal, Gupta, Kumar and Raju2011; Triki et al. Reference Triki, Porsezian, Choudhuri and Dinda2016, Reference Triki, Sun, Biswas, Zhou, Yıldırım, Zhong and Alshehri2022).

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(\varepsilon )$. The terms these generate are then included into the NLS equation (3.23) 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

(4.1)\begin{equation} {\rm i}B_T+\frac{\omega_0''}{2}B_{ZZ}-\varGamma |B|^2B-{\rm i}\varepsilon Q |B|^2B_Z = 0, \end{equation}

where the additional nonlinear correction has coefficient

(4.2)\begin{align} Q & = \frac{\mathcal{Z}_c\omega_{pi}^2m_en_{e,0}+\omega_{pe}^2m_in_{i,0}}{\mathcal{Z}_c\mu m_em_in_{e,0}n_{i,0} c^4 k^3\Delta D_\omega}(\omega-v_{\|}k-\varOmega_e)(\omega-u_{\|}k+\varOmega_i)\nonumber\\ & \quad \times\left[\frac{\mathcal{Z}_c\omega_{pe}^2n_{i,0}}{(\omega-v_{\|}k-\varOmega_e)^2}\left((\omega-v_{\|}k)^2-\varOmega_e(\omega-c_gk)\right)\left(1+\frac{(c_g-v_{\|})k}{\omega-v_{\|}k-\varOmega_e}\right)\right.\nonumber\\ & \quad \left.+\frac{\omega_{pi}^2n_{e,0}}{(\omega-u_{\|}k+\varOmega_i)^2}\left((\omega-u_{\|}k)^2+\varOmega_i(\omega-c_gk)\right)\left(1+\frac{(c_g-u_{\|})k}{\omega+\varOmega_i}\right)\right]. \end{align}

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 wave train solution to this amplitude equation, by

(4.3)\begin{equation} B = A(1+P)\exp({\rm i} [\delta k Z- (\omega_0'' \delta k^2/2)+ A^2(\varGamma-\varepsilon Q \delta k) T]), \end{equation}

for Stokes wave amplitude $A$, sideband wavenumber $\delta 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

(4.4)\begin{equation} \left(\sigma -\varepsilon \kappa QA^2\right)^2 = \kappa^2\left[\omega_0''(\varGamma-\varepsilon Q \delta k) +\frac{\omega_0''^2}{4}\kappa^2\right] \end{equation}

and thus there is a band of spectral wavenumbers $\kappa$ which generates perturbations that grow in time whenever

(4.5)\begin{equation} \omega_0''\varGamma \left(1-\frac{\varepsilon Q \delta k}{\varGamma}\right) <0. \end{equation}

We can see that the modulation criterion $\omega _0''\varGamma <0$ may persist if the bracketed expression is positive. This is ensured even if $\delta k>0$ as for parameter ranges operational in the magnetosphere one has $\varepsilon Q/\varGamma \ll 1$, and so $\delta k$ must be significantly large in order to do so which would violate the assumptions used to derive (4.1). 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

(4.6)\begin{equation} G = \max({\rm Im}(\sigma)) = |\varGamma+\varepsilon Q \delta k|, \quad {\rm at} \quad \kappa ={\pm} \sqrt{-\frac{2(\varGamma-\varepsilon Q \delta k)}{\omega_0''}}. \end{equation}

It is apparent now that the sideband wavenumber $\delta 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 $\varGamma$ is determined by the factor $\varDelta$. 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.

4.1. 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 non-dimensionalisations (2.29a–e) and (3.25), where $r \ll 1$, to simplify the coefficient of the term responsible for chirping to

(4.7)\begin{equation} Q = \frac{\alpha^4cn_{i,0}}{\mu m_e n_{e,0} \Delta K}\frac{(W - 1)^2(2W - 1)(2W^2 - \alpha^2 - 2W)(2W^3 - 2W^2 - \alpha^2)}{ (W^2 - \alpha^2 - W)(2W^3 - 4W^2 + \alpha^2 + 2W)^3}. \end{equation}

In this limit there is clearly a root of $Q$ at $W = \tfrac {1}{2}$, namely $\omega = \varOmega _e/2$, which is at the typical band-gap for WMC waves. The only other sign change of $Q$ is determined by $\varDelta$ (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\leq \tfrac {1}{2}$.

This observation of this root of $Q$ at $W=\tfrac {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 rises or falls 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 multimode interaction system is required in the spirit of either a coupled NLS/Manakov system (Manakov Reference Manakov1974; Kourakis, Shukla & Morfill Reference Kourakis, Shukla and Morfill2005b,Reference Kourakis, Shukla and Morfilla; Baronio et al. Reference Baronio, Degasperis, Conforti and Wabnitz2012) or a Zakharov/wave turbulence model (Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2000; Newell & Rumpf Reference Newell and Rumpf2011; David & Galtier Reference David and Galtier2022), 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.

4.2. 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 (Reference Dias and Iooss1993) and consider a solution of (4.1) of the form

(4.8)\begin{equation} B = R(\xi){\rm e}^{{\rm i} (\nu Z-\sigma T)}{\rm e}^{{\rm i} \phi(\xi)}, \quad \xi = Z-VT, \end{equation}

which is formed of three parts – the amplitude function $R$, the sideband wave component represented by the second term and a phase function $\phi$, 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,

(4.9)\begin{equation} \tilde{\omega} = \frac{\partial}{\partial t}(\sigma T-\nu Z-\phi) = \varepsilon(c_g\nu+\varepsilon\sigma) +\varepsilon(c_g+\varepsilon V)\phi'. \end{equation}

From this one may identify a sweep rate, defined as the rate of change of the frequency in time,

(4.10)\begin{equation} \mathscr{S} = \frac{\partial \tilde{\omega}}{\partial t} =-\varepsilon^2(c_g+\varepsilon V)^2\phi''. \end{equation}

Our aim is 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 (4.1) and solve the resulting ordinary differential equations (ODEs).

Starting with our guess at a solution, and splitting the resultant system into real and imaginary parts, we have the equations

(4.11)\begin{gather} \sigma R+V\phi'R+\frac{\omega_0''}{2}(R''-(\nu+\phi')^2R)-\varGamma R^3+\varepsilon Q (\nu+\phi')R^3 = 0, \end{gather}

(4.12)\begin{gather}-VR'+\frac{\omega_0''}{2} (2R'\phi'+R\phi''+2 \nu R')-\varepsilon Q R^2R' = 0. \end{gather}

We start by taking the second equation, (4.12), and multiplying by $R$, which makes it an exact derivative. Integrating gives

(4.13)\begin{equation} -\frac{V}{2}R^2+\frac{\omega_0''}{2} (R^2\phi'-\nu R^2)-\frac{\varepsilon Q}{4}R^4 = I \end{equation}

for some constant $I$. It will be convenient 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 $\xi \to \pm \infty$. This yields that $I=0$ and allows one to manipulate (4.13) to show that the sweep rate $\mathscr {S}$ is

(4.14)\begin{equation} \mathscr{S} =-\frac{\varepsilon^3(c_g+\varepsilon V)^2 Q}{2\omega_0''}U'. \end{equation}

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 $\nu = 0$ to simplify the analysis, and it can be reintroduced under suitable mappings. In this case, we get a quartic potential,

(4.15)\begin{equation} (U')^2+B_2 U^2+B_3U^3+B_4U^4 \equiv (U')^2-\mathscr{V}(U) = 0 \end{equation}

with

(4.16a–c)\begin{equation} B_2 =\frac{4V^2}{\omega_0''^2}+\frac{8\sigma}{\omega_0''} , \quad B_3 =\frac{4}{\omega_0''}\left(\frac{\varepsilon VQ}{\omega_0''}-\varGamma\right), \quad B_4 = \left(\frac{\varepsilon Q}{\omega_0''}\right)^2. \end{equation}

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 $\mathscr {V}(U)>0$. As $B_4>0$, we need $B_2<0$ and so

(4.17)\begin{equation} \frac{4V^2}{\omega_0''^2}+\frac{8\sigma}{\omega_0''} <0. \end{equation}

Since $\sigma$ is a free real-valued parameter, this is always possible to satisfy. Following Kamchatnov et al. (Reference Kamchatnov, Kuo, Lin, Horng, Gou, Clift, El and Grimshaw2012) let us factorise this ODE as follows:

(4.18)\begin{equation} U' =\left\vert\frac{\varepsilon Q}{\omega_0''}\right\vert\sqrt{U^2(U-U_-)(U_+-U)}, \end{equation}

where the roots are given by

(4.19)\begin{equation} U_\pm=\frac{2\omega_0''}{\varepsilon^2Q^2}\left[\varGamma-\frac{\varepsilon V Q}{ \omega_0''} \pm \sqrt{\varGamma^2-\frac{2\varepsilon Q}{\omega_0''}(V\varGamma+\varepsilon Q \sigma)} \right]. \end{equation}

The negative subscript denotes the negative root, and plus the positive root, so that $U_-<0< U_+$. Thus, the expression in (4.18) 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\leq U< U_+$. The solution to this ODE is of the form of a Gardner/extended Korteweg–de Vries soliton (Grimshaw et al. Reference Grimshaw, Pelinovsky, Talipova and Kurkina2010; Kamchatnov et al. Reference Kamchatnov, Kuo, Lin, Horng, Gou, Clift, El and Grimshaw2012),

(4.20)\begin{equation} U = \frac{A}{(1+B)\cosh^2 \varTheta - 1} \end{equation}

with

(4.21a–c)\begin{equation} A =-U_-{>}0, \quad B =-\frac{U_-}{U_+}>0, \quad \varTheta = \sqrt{-\frac{2 \sigma \omega_0''+V^2}{\omega_0''^2}}(Z-VT), \end{equation}

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 wave trains 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.

5. 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 (4.1) using a time stepping procedure. For computational ease and to better identify the effects of system parameters on the sweep rate, we rescale (4.1), highlighting two cases dependant on the sign of $\omega _0''$. We employ the scalings

(5.1a–c)\begin{equation} X = {\rm sign}(\omega_0'') Z, \quad \tau = {\rm sign}(\omega_0'') \frac{\omega_0''}{2}T , \quad B = \begin{cases} \sqrt{\dfrac{|\omega_0''|}{2|\varGamma|}}A^* & \omega_0''<0, \\ \sqrt{\dfrac{|\omega_0''|}{2|\varGamma|}}A & \omega_0''>0, \end{cases} \end{equation}

to transform (4.1) into

(5.2)\begin{equation} {\rm i}A_\tau+A_{XX}-\sigma\left(|A|^2A+\frac{{\rm i}\varepsilon Q}{\varGamma}|A|^2A_X\right) = 0, \quad {\rm with} \quad \sigma = {\rm sign}(\omega_0''\varGamma). \end{equation}

This reduces the problem to a single tuneable parameter, which we can see from the non-dimensionalisations earlier that

(5.3)\begin{equation} \frac{\varepsilon Q}{\varGamma} = \frac{\varepsilon c}{\varOmega K}\frac{(2W-1)(2W^2-2W-\alpha^2)}{2W^3-4W^2+2W+\alpha^2}, \end{equation}

thus reducing the tuneable parameters to just $\varepsilon, W, \varOmega$ (or equivalently, $B_0$) and $\alpha$ 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\sim 6$. This corresponds to

(5.4a–d)\begin{equation} \alpha = 7.2, \quad B_0 =1.4\times 10^{-7} , \quad n_e = 1\times 10^7, \quad \varepsilon = O(10^{-5}). \end{equation}

In addition to this, we focus on wave frequencies lying in the lower frequency band range of $0.1 \varOmega$–$0.5\varOmega$.

We advance (5.2) in time using an exponential time differencing method with fourth-order Runge–Kutta time stepping (ETDRK4) (Cox & Matthews Reference Cox and Matthews2002; Kassam & Trefethen Reference Kassam and Trefethen2005), using periodic boundary conditions to take advantage of the speed and spectral accuracy of Fourier-based schemes. To initialise a multielement solution and remain close to the analytic solution found, we initialise the simulations with the wave packet

(5.5)\begin{equation} A(X,0) = \sqrt{\frac{1}{A_0\cosh\left(\dfrac{X}{\varLambda}\right)^2-1}} {\rm e}^{{\rm i} p X}, \end{equation}

which has amplitude ${1}/({A_0-1})$, width $\varLambda$ 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 $\varLambda$ 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 (5.2) 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 (4.14) would appear to be the source for the frequency sweeps that emerge numerically. The sweep rate decreases over each simulation, as the subelements grow farther 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 (5.3) enhance the onset of frequency sweeping, which further points to the role that the modified term in the NLS plays in generating frequency shifts in the WMC waves.

Figure 2. Examples of the power spectrum $|B_W|^2/\delta \omega$ generated by the parameter choices $(\varepsilon,n_e,\alpha,\varOmega ) = (4\times 10^{-5},1\times 10^7,7.2,2.49\times 10^4)$ and (a) $W=0.2, k_BT =8.02$ KeV, (b) $W= 0.35, k_BT =7.37$ KeV.

It is worth commenting on the contrast between the simulations here and the more conventional studies using particle-in-cell simulations. These have demonstrated chorus wave frequency chirping in a one-dimensional uniform field, but using reflective boundary conditions for particles and open boundary conditions for electromagnetic waves (Wu et al. Reference Wu, Tao, Zonca, Chen and Wang2020). In that work, the open boundary conditions for the waves facilitate the chirping behaviour in the sense that they allow electromagnetic perturbations to leave the domain and so not to return and interfere with the spatially localised wave–particle interaction that is creating the chirp. The mathematical analysis that is presented in this work is based on the assumption of an infinite uniform plasma in a uniform background magnetic field, and we employ periodic boundary conditions for all quantities in our numerical experiments. However, our experiments still permit the chirping mechanism to be fully local/isolated in the same way as the work by Wu et al. (Reference Wu, Tao, Zonca, Chen and Wang2020), since the spatial dimensions of the domain $[0\le x < L]$ are sufficiently large that the relevant phase and group velocities do not traverse a full box-length $L$ over the time scale of chirping events and indeed the whole numerical experiment. This is to say that $v_{g}\ll LT^{-1} , L \omega ^{-1}\,\textrm {d}\omega /\textrm{d}t$. Furthermore, the paper by Wu et al. (Reference Wu, Tao, Zonca, Chen and Wang2020) demonstrates that in their regime, it is the background non-uniformity that establishes a symmetry breaking in order to preferentially permit either rising or falling tones, based on the curvature of the background magnetic field. It is therefore interesting that our results demonstrate a preferential chirping direction (rising tone) in the case of a uniform field. Once again, we state again that our standalone mechanism is of a different origin and nature, and it will be important in future studies to investigate the relationship between this, and others.

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 3(a) akin to the Kuznetsov–Ma soliton found in the NLS (Ma Reference Ma1979; Akhmediev, Eleonskii & Kulagin Reference Akhmediev, Eleonskii and Kulagin1985), albeit with a zero background. The structure of maximum amplification for widths that generate two to three main packets also bears some semblance to higher-order rogue wave solutions (Chabchoub et al. Reference Chabchoub, Hoffmann, Onorato and Akhmediev2012; Slunyaev et al. Reference Slunyaev, Pelinovsky, Sergeeva, Chabchoub, Hoffmann, Onorato and Akhmediev2013), 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.

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

We may also use (3.22) to explore the effect on both the number density and the parallel electron 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:

(5.6)\begin{equation} E_{KE} = \frac{m_e}{2}n_e\|{\boldsymbol V}_e-v_{\|}\hat{\boldsymbol{z}}\|^2+\frac{k_BT}{2}(n_e-n_{e,0}), \end{equation}

which allows us to determine which stage of the WMC evolution energises the electrons. We visualise these quantities in figure 3(c), 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 Reference Tsurutani and Smith1974; Li et al. Reference Li, Thorne, Bortnik, Tao and Angelopoulos2012; Gao et al. Reference Gao, Li, Thorne, Bortnik, Angelopoulos, Lu, Tao and Wang2014), 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.

Figure 4. (a) Snapshots of the magnetic field wave versus the power spectra for the time series snapshot for the parameter choice $(\varepsilon,n_e,\alpha,\varOmega, W, k_BT) = (4\times 10^{-5},1\times 10^7,7.2,2.49\times 10^4,0.35, 7.373$ KeV). (b) 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.

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 modified NLS equation (4.1) appears to be unable to produce a single isolated rising tone event that has been successfully created via other theoretical modelling (e.g. Nunn et al. Reference Nunn, Omura, Matsumoto, Nagano and Yagitani1997; Omura et al. Reference Omura, Katoh and Summers2008; Tao Reference Tao2014). 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. Reference Agapitov, Krasnoselskikh, Zaliznyak, Angelopoulos, Le Contel and Rolland2011; Gao et al. Reference Gao, Chen, Lu, Chen, Chen and Wang2022) 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.

6. Concluding remarks

This paper has provided an overview of the formation and dynamics of parallel propagating WMC waves with rising tone. We have outlined that 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_g^2-c_s^2$. 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 here emphasise that the mechanisms proposed in this paper (namely chirping and potential contributions to the power-gap phenomenon as a result of the ponderomotive force and non-zero ion motion) are not proposed to be replacing other proposed mechanisms, such as those that may even include fixed ions in, for example, Katoh & Omura (Reference Katoh and Omura2007) and Tao et al. (Reference Tao, Zonca and Chen2021). Future work should investigate to what extent different viable mechanisms such as these and others interact and compete in different regimes.

Furthermore, we note that the power-gap phenomenon is not always observed in spacecraft measurements, for example, see Teng et al. (Reference Teng, Tao and Li2019). It is interesting to note that the occurrence rate of ‘no-gap’ whistler waves has been observed to peak off the equator ($|\textrm {MLAT}|\approx 8\unicode{x2013} 10^\circ$), i.e. in regimes where magnetic curvature cannot be ignored. It will be interesting in future work to see if we can obtain similar results. We emphasise once more, the mechanism that we propose is shown to be viable within the context of a uniform plasma with uniform background magnetic field, and it is not yet known how this mechanism will change in different geometry and regimes, and/or compete with other viable mechanisms.

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 it 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 $\varepsilon$ 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^\circ$ (close to parallel) or $70^\circ$ 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 (5.3). 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.