2.1 Base model

2.1.1 Incompressible resistive MHD

As a base model, we assume incompressible resistive MHD with homogeneous mass density $\rho = \text {const.}$ The corresponding governing equations are

(2.1a)\begin{gather} \partial_t \boldsymbol{v} + (\boldsymbol{v} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{v}= (\boldsymbol{b} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{b} - \boldsymbol{\nabla} P + \nu \nabla^2 \boldsymbol{v}, \end{gather}

(2.1b)\begin{gather}\partial_t \boldsymbol{b} = \boldsymbol{\nabla} \times (\boldsymbol{v} \times \boldsymbol{b}) + \eta \nabla^2 \boldsymbol{b}. \end{gather}

Here $\boldsymbol {v}$ is the fluid velocity; $\boldsymbol {b} \doteq \boldsymbol {B}/\sqrt {4{\rm \pi} \rho }$ is the local Alfvén velocity (the symbol $\doteq$ denotes definitions), i.e. rescaled magnetic field $\boldsymbol {B}$; $P \doteq (P_{\text {kin}} + B^2/8{\rm \pi} )/\rho$ is the normalized total pressure, with $P_{\text {kin}}$ being the kinetic pressure. Also, $\nu$ is the viscosity and $\eta$ is the resistivity, both of which are assumed either constant or negligible (except at small enough scales). Due to incompressibility and magnetic Gauss's law, one also has

(2.1c)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} = 0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{b} = 0. \end{equation}

In order to put (2.1) in a more symmetric form, let us rewrite them in terms of the two Elsässer fields $\boldsymbol {z}^\pm$ (Elsasser Reference Elsasser1950), which are also solenoidal,

(2.2a,b)\begin{equation} \boldsymbol{z}^\pm \doteq \boldsymbol{v} \pm \boldsymbol{b},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{z}^\pm = 0. \end{equation}

This leads to two coupled equations for $\boldsymbol {z}^\pm$,

(2.3)\begin{equation} \partial_t \boldsymbol{z}^\pm = - (\boldsymbol{z}^\mp\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{z}^\pm - \boldsymbol{\nabla} P +\nu_+ \nabla^2 \boldsymbol{z}^\pm + \nu_- \nabla^2 \boldsymbol{z}^\mp, \end{equation}

where $\nu _+ \doteq (\nu + \eta )/2$ and $\nu _- \doteq (\nu - \eta )/2$. The normalized total pressure $P$ can be found as follows. By taking the divergence of (2.1a) and using (2.2a,b), one obtains

(2.4)\begin{equation} \nabla^2 P = - \boldsymbol{\nabla}\boldsymbol{\cdot}[(\boldsymbol{z}^\mp \boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{z}^\pm]. \end{equation}

Let us introduce the wavevector operator $\hat {\boldsymbol {k}} \doteq - \mathrm {i} \boldsymbol {\nabla }$, so $\hat {k}^2 \doteq \smash {\hat {\boldsymbol {k}}}^2 = - \nabla ^2$. Let us also assume some appropriate (say, periodic) boundary conditions. Then, (2.4) yields

(2.5)\begin{equation} P = - \hat{k}^{ - 2}\hat{\boldsymbol{k}}\boldsymbol{\cdot}[(\boldsymbol{z}^\mp \boldsymbol{\cdot} \hat{\boldsymbol{k}})\boldsymbol{z}^\pm] + \text{const.}, \end{equation}

whence $\boldsymbol {\nabla } P$ in (2.3) can be expressed through $\boldsymbol {z}^\pm$. Alternatively, $\boldsymbol {\nabla } P$ can be eliminated from (2.3) by taking the curl of this equation and rewriting the result in terms of the Elsässer vorticities

(2.6)\begin{equation} \boldsymbol{w}^\pm \doteq \boldsymbol{\nabla} \times \boldsymbol{z}^\pm. \end{equation}

Indeed, due to (2.2a,b), one can express $\boldsymbol {z}^\pm$ using a vector potential $\boldsymbol {a}^\pm$ such that $\boldsymbol {z}^\pm = \boldsymbol {\nabla } \times \boldsymbol {a}^\pm$. Let us assume the gauge such that $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {a}^\pm = 0$. Then,

(2.7)\begin{equation} \boldsymbol{w}^\pm = \boldsymbol{\nabla} \times (\boldsymbol{\nabla} \times \boldsymbol{a}^\pm) = - \nabla^2 \boldsymbol{a}^\pm \equiv \hat{k}^2 \boldsymbol{a}^\pm, \end{equation}

whence

(2.8)\begin{equation} \boldsymbol{z}^\pm = \mathrm{i} \hat{k}^{ - 2} (\hat{\boldsymbol{k}} \times \boldsymbol{w}^\pm). \end{equation}

(Here, we assume that $\boldsymbol {z}^\pm$ have zero spatial average, so $\hat {k}^2$ is invertible.) Then, (2.3) can be expressed through $\boldsymbol {w}^\pm$ alone,

(2.9)\begin{equation} \partial_t \boldsymbol{w}^\pm = - \hat{\boldsymbol{k}} \times\{[(\hat{\boldsymbol{k}} \times \hat{k}^{ - 2}\boldsymbol{w}^\mp)\boldsymbol{\cdot} \hat{\boldsymbol{k}}] (\hat{\boldsymbol{k}}\times\hat{k}^{ - 2}\boldsymbol{w}^\pm)\}-\hat{k}^2(\nu_+ \boldsymbol{w}^\pm + \nu_- \boldsymbol{w}^\mp). \end{equation}

2.1.2 The 2-D collisionless model

For simplicity, below we adopt a 2-D model, meaning that $\boldsymbol {z}^\pm$ lie in the $(x, y)$ plane and $\partial _z = 0$. In this case, the only potentially non-zero component of $\boldsymbol {w}^\pm$ is the $z$-component, $w^\pm \doteq (\boldsymbol {\nabla }\times \boldsymbol {z}^\pm )_z$. We also forgo an explicit treatment of the viscosities $\nu _\pm$ below. (However small, though, these viscosities remain non-negligible in that they determine absorbing boundary conditions at infinity in the spectral space; see § 4.) Then, the vector equation (2.9) can be replaced with a scalar equation for $w^\pm$,

(2.10)\begin{equation} \partial_t w^\pm = \epsilon_{lmz} \left[\left(\frac{\hat{k}_m}{\hat{k}^2}w^\mp \right) (\hat{k}_l w^\pm)- \left(\frac{\hat{k}_n\hat{k}_l}{\hat{k}^2} w^\mp\right) \left(\frac{\hat{k}_n\hat{k}_m}{\hat{k}^2} w^\pm\right)\right], \end{equation}

where $\epsilon _{lmz}$ is the Levi–Civita symbol with the third index fixed to $z$. Note that the right-hand side of (2.10) vanishes whenever $w^+$ or $w^-$ is zero. This means that any stationary $w^+$ is a solution as long as $w^-$ is zero and vice versa.

2.1.3 Fourier representation

Assuming the Fourier representation

(2.11)\begin{equation} w^\pm = \sum_{\boldsymbol{k}} w_{\boldsymbol{k}}^\pm(t) \exp({\mathrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{x}}), \end{equation}

where $\boldsymbol {x} \doteq (x, y)$, one arrives at the following equation for the Fourier coefficients $w_{\boldsymbol {k}}^\pm$:

(2.12)\begin{equation} \partial_t w^\pm_{\boldsymbol{k}}= \sum_{\boldsymbol{k}_1,\boldsymbol{k}_2} T(\boldsymbol{k}_1,\boldsymbol{k}_2) w^\mp_{\boldsymbol{k}_1} w^\pm_{\boldsymbol{k}_2} \delta_{\boldsymbol{k},\boldsymbol{k}_2+\boldsymbol{k}_1}, \end{equation}

where the dot is the time derivative, $\delta$ is the Kronecker delta, $T(\boldsymbol {k}_1,\boldsymbol {k}_2)$ are the coupling coefficients given by

(2.13)\begin{equation} T(\boldsymbol{k}_1,\boldsymbol{k}_2)= [\boldsymbol{e}_{z} \boldsymbol{\cdot} (\boldsymbol{k}_2\times\boldsymbol{k}_1)]\frac{(\boldsymbol{k}_2 + \boldsymbol{k}_1)\boldsymbol{\cdot}\boldsymbol{k}_2}{k_1^2 k_2^2}, \end{equation}

and $\boldsymbol {e}_{j}$ is the unit vector along the $j$th axis. Accordingly, the kinetic and magnetic energy are given by

(2.14a)\begin{gather} E_v = \frac{\rho}{8}\sum_{\boldsymbol{k}}\frac{1}{k^2}|w_{\boldsymbol{k}}^+ + w_{\boldsymbol{k}}^-|^2 \equiv \sum_{\boldsymbol{k}} E_{{v},\boldsymbol{k}}, \end{gather}

(2.14b)\begin{gather}E_b = \frac{\rho}{8}\sum_{\boldsymbol{k}}\frac{1}{k^2}|w_{\boldsymbol{k}}^+ - w_{\boldsymbol{k}}^-|^2 \equiv \sum_{\boldsymbol{k}} E_{{b},\boldsymbol{k}}, \end{gather}

where $E_{{v},\boldsymbol {k}}$, and $E_{{b},\boldsymbol {k}}$ are the corresponding spectral energy densities. Then, the total energy density is

(2.15)\begin{equation} E = \frac{\rho}{4}\sum_{\boldsymbol{k}}\frac{1}{k^2}(|w_{\boldsymbol{k}}^+|^2 + |w_{\boldsymbol{k}}^-|^2) \equiv \sum_{\boldsymbol{k}} E_{\boldsymbol{k}}, \end{equation}

where $E_{\boldsymbol {k}}$ is the spectral total-energy density.

2.2 Harmonic coupling mediated by a primary structure

2.2.1 Equation for the modulational spectrum

Let us explore modulational stability of a Fourier harmonic with a given wavevector $\boldsymbol {k} = \boldsymbol {p}$ and the corresponding order-one Fourier amplitudes $w_{\boldsymbol {p}}^\pm$. We will call it the primary harmonic or, primary structure. Let us choose axes such that $\boldsymbol {p} = p\boldsymbol {e}_{y}$ and assume a perturbation with a wavevector $\boldsymbol {q}$. As can be seen from (2.13), non-trivial coupling requires a component of $\boldsymbol {q}$ perpendicular to $\boldsymbol {p}$, and the component of $\boldsymbol {q}$ parallel to $\boldsymbol {p}$ does not qualitatively affect the modulational dynamics.Footnote ³ For simplicity then, let us assume $\boldsymbol {q} = q\boldsymbol {e}_{x}$. Let us also assume initial conditions such that at $t=0$, $w_{\boldsymbol {p}}^\pm = O(1)$, and $w_{\boldsymbol {q}+n\boldsymbol {p}}^\pm = O(\epsilon )$, where $\epsilon \ll 1$ is a small parameter.Footnote ⁴ (Remember that $O$($\epsilon$) means, loosely speaking, ‘scales as $\epsilon$ or smaller’ but not necessarily ‘of order $\epsilon$’.) This ordering is not guaranteed to hold as the system evolves (for example, a MI corresponds to exponential growth of $w_{\boldsymbol {q}+n\boldsymbol {p}}^\pm$), but it will hold at least for some initial linear stage; i.e. equations resulting from this ordering will describe the linear growth of a MI but not its saturation.

For early times such that primary mode remains dominant, the largest contributions to the sum in (2.12) have $\boldsymbol {k}_1=\boldsymbol {p}$ or $\boldsymbol {k}_2=\boldsymbol {p}$. To first order in $\epsilon$ we then have

(2.16)\begin{align} \partial_t w^\pm_{\boldsymbol{k}} & \approx T(\boldsymbol{p},\boldsymbol{k}_-)w^\mp_{\boldsymbol{p}} w^\pm_{\boldsymbol{k}_-} + T(-\boldsymbol{p},\boldsymbol{k}_+)w^\mp_{-\boldsymbol{p}} w^\pm_{\boldsymbol{k}_+} \nonumber\\ & \quad + T(\boldsymbol{k}_-,\boldsymbol{p})w^\mp_{\boldsymbol{k}_-} w^\pm_{\boldsymbol{p}} + T(\boldsymbol{k}_+,-\boldsymbol{p})w^\mp_{\boldsymbol{k}_+} w^\pm_{-\boldsymbol{p}}, \end{align}

where $\boldsymbol {k}_{\pm }\doteq \boldsymbol {k} \pm \boldsymbol {p}$ and $\boldsymbol {k}= \boldsymbol {q}+n\boldsymbol {p}$ with integer $n$. (All the five terms in (2.16) are of order $\epsilon$, so, whether or not one chooses to introduce $\epsilon$ explicitly in this equation, the small parameter cancels out in any case. For the same reason, (2.16) is invariant with respect to rescaling $w^\pm _{\boldsymbol {k}} \to Cw^\pm _{\boldsymbol {k}}$, where $C$ is any non-zero constant.) For all other $\boldsymbol {k}$, one obtains $\partial _t w^\pm _{\boldsymbol {k}} = O(\epsilon ^2)$ and therefore negligible within the assumed approximation. In particular, this means that $w^\pm _{\boldsymbol {p}}$ can be considered fixed; i.e. $A^\pm \doteq w^\pm _{\boldsymbol {p}} \approx \text {const}$. In terms of the Elsässer fields, this corresponds to a primary structure

(2.17)\begin{equation} \boldsymbol{z}^\pm = - \frac{2}{p}\mathcal{A}^\pm \sin(py + \theta^\pm)\boldsymbol{e}_{x}, \end{equation}

where $\mathcal {A^\pm } \doteq |A^\pm |$ and $\theta ^\pm \doteq \arg A^\pm$.

For the modulation energy density and the primary-wave energy density,

(2.18a)\begin{gather} E_{\text{mod}} \doteq \sum_{\boldsymbol{k}=\boldsymbol{q}+n\boldsymbol{p}} E_{\boldsymbol{k}}=\sum_n \frac{\rho}{4}\sum_{\sigma=\pm1}\frac{|w_{\boldsymbol{q}+n\boldsymbol{p}}^{\sigma}|^2}{q^2 + n^2p^2}, \end{gather}

(2.18b)\begin{gather}E_{\text{pri}} \doteq E_{\boldsymbol{p}}=\frac{\rho}{4}\sum_{\sigma=\pm1} \frac{|w_{\boldsymbol{p}}^{\sigma} |^2}{p^2}, \end{gather}

one has

(2.19a)\begin{gather} \frac{\mathrm{d} E_{\text{mod}}}{\mathrm{d} t} = \sum_n \frac{\rho}{4}\sum_{\sigma=\pm1} \frac{w_{\boldsymbol{q} +n\boldsymbol{p}}^{\sigma*}\partial_t w_{\boldsymbol{q} + n\boldsymbol{p}}^{\sigma}}{q^2 + n^2p^2}+ \text{c.c.}, \end{gather}

(2.19b)\begin{gather}\frac{\mathrm{d} E_{\text{pri}}}{\mathrm{d} t} = \frac{\rho}{4}\sum_{\sigma=\pm1} \frac{w_{\boldsymbol{p}}^{\sigma*}\dot w_{\boldsymbol{p}}^{\sigma}}{p^2} + \text{c.c.}, \end{gather}

where ‘c.c.’ stands for ‘complex conjugate’. Using

(2.20)\begin{equation} \partial_t w_{\boldsymbol{p}}^{\sigma}=\sum_{\boldsymbol{k}=\boldsymbol{q}+n\boldsymbol{p}} T(\boldsymbol{k},\boldsymbol{p-k}) w_{\boldsymbol{k}}^{-\sigma} w_{\boldsymbol{p-k}}^\sigma, \end{equation}

for $\boldsymbol {k}=\boldsymbol {q}+n\boldsymbol {p}$, one finds that our approximate equations conserve the total energy density $E_{\text {tot}} = E_{\text {pri}} + E_{\text {mod}}$ exactly within this approximation. (For brevity, from now on we will refer to the energy densities simply as energies.) We will assume that the initial modulational energy is negligible, so $E_{\text {tot}} = \rho \mathcal {A}^2/4p^2$, where $\mathcal {A}^2\doteq (\mathcal {A}^+)^2+(\mathcal {A}^-)^2$, and we will also use the following dimensionless energies:

(2.21a)\begin{gather} \mathcal{E}_{\text{pri}}=E_{\text{pri}}/E_{\text{tot}}, \end{gather}

(2.21b)\begin{gather}\mathcal{E}_{\text{mod}}=E_{\text{mod}}/E_{\text{tot}}, \end{gather}

(2.21c)\begin{gather}\mathcal{E}_v=E_v/E_{\text{tot}}, \end{gather}

(2.21d)\begin{gather}\mathcal{E}_v=E_b/E_{\text{tot}}. \end{gather}

2.2.2 Dimensionless equations

In order to reduce the number of free parameters, let us perform a variable transformation $w_{\boldsymbol {q}+n\boldsymbol {p}}^\pm \mapsto \psi _n^\pm$,

(2.22)\begin{equation} w_{\boldsymbol{q}+n\boldsymbol{p}}^\pm = \sqrt{\frac{\mathcal{A}^3(q^2+n^2p^2)}{\mathcal{A}^\mp p^2}} \exp\left(\mathrm{i} n \frac{\theta^+ + \theta^-}{2}\right) \psi_n^\pm, \end{equation}

where $\psi _n^\pm = O(\epsilon )$. Then, (2.16) becomes

(2.23)\begin{align} \partial_t \psi_n^\pm & = \mathcal{A}^\mp \mathrm{e}^{{\mp}\mathrm{i}\theta} \alpha_n^- \psi_{n-1}^\pm + \sqrt{\mathcal{A}^+\mathcal{A}^-}\mathrm{e}^{{\pm} \mathrm{i}\theta}\beta_n^- \psi_{n-1}^\mp \nonumber\\ & \quad + \mathcal{A^\mp}\mathrm{e}^{{\pm}\mathrm{i}\theta}\alpha_n^+ \psi_{n+1}^\pm + \sqrt{\mathcal{A}^+\mathcal{A}^-}\mathrm{e}^{{\mp} \mathrm{i}\theta}\beta_n^+ \psi_{n+1}^\mp, \end{align}

where $\theta \doteq (\theta ^+-\theta ^-)/2$ and $\alpha _n^\pm$ and $\beta _n^\pm$ are given by

(2.24a)\begin{gather} \alpha_n^\pm \doteq \mp r\frac{r^2 + n(n \pm 1)}{\sqrt{(r^2+n^2)(r^2 + (n \pm 1)^2})}, \end{gather}

(2.24b)\begin{gather}\beta_n^\pm \doteq -r\frac{n}{\sqrt{(r^2+n^2)(r^2 + (n \pm 1)^2)}}. \end{gather}

Note that these coefficients depend only on the ratio $r \doteq q/p$ rather than on $q$ and $p$ separately. From (2.23), one can also see that the phases $\theta ^+$ and $\theta ^-$ per se are, expectedly, unimportant and the dynamics is affected by these phases only in the combination $\theta ^+ - \theta ^-$.

It is also convenient for our purposes to introduce

(2.25)\begin{equation} \tau \doteq \mathcal{A} t \end{equation}

as the new dimensionless time, as well as $\phi \doteq \arctan (\mathcal {A}^-/\mathcal {A}^+) \in [0, {\rm \pi}/2]$, so that

(2.26a,b)\begin{equation} \mathcal{A}^+ = \mathcal{A}\cos\phi,\quad \mathcal{A}^- = \mathcal{A}\sin\phi. \end{equation}

In particular, $\phi = 0$ and $\phi = {\rm \pi}/2$ correspond to structures consisting only of $w_{\boldsymbol {p}}^+$ and $w_{\boldsymbol {p}}^-$, respectively, and $\phi = {\rm \pi}/4$ corresponds to a ‘balanced’ background with $\mathcal {A}^+ = \mathcal {A}^-$. For simplicity, $\phi = {\rm \pi}/4$ will be assumed from now on throughout the paper (except where specified otherwise). In this case (for $\phi = {\rm \pi}/4$), one has

(2.27)\begin{equation} \frac{E_{\text{pri}, b}}{E_{\text{pri}, v}} = \frac{|w_{\boldsymbol{p}}^+ - w_{\boldsymbol{p}}^-|^2}{|w_{\boldsymbol{p}}^+ + w_{\boldsymbol{p}}^-|^2} = \tan^2\theta, \end{equation}

where $E_{\text {pri}, v}$ and $E_{\text {pri}, b}$ are the kinetic and magnetic energy densities in the primary mode, respectively. Hence, $\theta$ can be understood as the parameter that determines the relative weights of the kinetic and magnetic-field energy in the primary-mode energy. In particular, $\theta = 0$ corresponds to a primary structure that involves only velocity perturbations, while $\theta = {\rm \pi}/2$ corresponds to a primary structure that involves only magnetic-field perturbations.

Using this notation, one can express (2.23) in the following compact form:

(2.28)\begin{equation} \dot{\boldsymbol{\psi}}_n = \sum_{\sigma ={\pm} 1}\boldsymbol{G}_{n}^\sigma\boldsymbol{\psi}_{n+\sigma} \end{equation}

(the dot denotes $\partial _\tau$), or even more compactly as

(2.29)\begin{equation} \dot{\boldsymbol{\psi}} = \boldsymbol{M}\boldsymbol{\psi}. \end{equation}

Here, $\boldsymbol {\psi }_n \doteq (\psi _n^-, \psi _n^+)^\intercal$ is a two-component column vector (the symbol $^\intercal$ denotes transposition), $\boldsymbol {\psi } \doteq (\unicode{x22EF}, \boldsymbol {\psi }_{-1}, \boldsymbol {\psi }_0,\boldsymbol {\psi }_1 \unicode{x22EF} )^\intercal$ is an infinite-dimensional block vector consisting of $\boldsymbol {\psi }_n$, the matrices $\boldsymbol {G}_{n}^\sigma$ are given by

(2.30)\begin{equation} \boldsymbol{G}_{n}^\sigma \doteq \begin{pmatrix} \alpha_n^\sigma\cos\phi \exp({-\mathrm{i}\sigma\theta}) & \beta_n^\sigma\sqrt{\cos\phi\sin\phi}\,\mathrm{e}^{\mathrm{i}\sigma\theta}\\ \beta_n^\sigma \sqrt{\cos\phi\sin\phi}\exp({-\mathrm{i}\sigma\theta}) & \alpha_n^\sigma\sin\phi \,\mathrm{e}^{\mathrm{i}\sigma\theta} \end{pmatrix} \end{equation}

and $\boldsymbol {M}$ is the following block matrix:

(2.31)\begin{equation} \boldsymbol{M}\doteq\begin{pmatrix} \ddots & & & \vdots & & & \unicode{x22F0}\\ & 0 & \boldsymbol{G}_{ - 2}^+ & 0 & 0 & 0 & \\ & \boldsymbol{G}_{ - 1}^- & 0 & \boldsymbol{G}_{ - 1}^+ & 0 & 0 & \\ \unicode{x22EF} & 0 & \boldsymbol{G}_0^- & 0 & \boldsymbol{G}_0^+ & 0 & \unicode{x22EF} \\ & 0 & 0 & \boldsymbol{G}_1^- & 0 & \boldsymbol{G}_1^+\\ & 0 & 0 & 0 & \boldsymbol{G}_2^- & 0 & \\ \unicode{x22F0} & & & \vdots & & & \ddots \end{pmatrix}. \end{equation}

It is readily seen then that the system's dynamics is controlled only by the dimensionless parameters $r$, $\theta$ and $\phi$, while $\mathcal {A}$ and $p$ determine only its characteristic temporal and spatial scales, respectively.

The dimensionless modulational energy (2.21b), or

(2.32)\begin{equation} \mathcal{E}_{\text{mod}} \doteq\sum_\sigma \frac{\mathcal{A}}{\mathcal{A}^{-\sigma}} |\psi_n^\sigma|^2\equiv \sum_n \mathcal{E}_n, \end{equation}

is governed by

(2.33)\begin{equation} \dot{\mathcal{E}}_n=\sum_\sigma (I_n^\sigma+F_n^\sigma) +\text{c.c.}, \end{equation}

where $I_n$ and $F_n$ are given by

(2.34)\begin{gather} I_n^\sigma = \frac{\mathcal{A}\sqrt{\mathcal{A}^\sigma}}{(\mathcal{A}^{-\sigma})^{3/2}} (\beta_n^- \mathrm{e}^{\mathrm{i}\sigma\theta} \psi_{n-1}^{-\sigma} + \beta_n^+\exp({-\mathrm{i}\sigma\theta})\psi_{n+1}^{-\sigma}) \psi_n^{\sigma *}, \end{gather}

(2.35)\begin{gather}F_n^\sigma = \frac{\mathcal{A}}{\mathcal{A}^{-\sigma}}(\alpha_n^- \exp({-\mathrm{i}\sigma\theta}) \psi_{n-1}^\sigma + \alpha_n^+ \mathrm{e}^{\mathrm{i}\sigma\theta} \psi_{n+1}^\sigma) \psi_n^{\sigma*}. \end{gather}

Because $\sum _n F_n^\sigma = 0$, the terms $F_n^\sigma$ represent the energy flux that is carried along the modulation spectrum and is conserved within each subchannel $\sigma$. In contrast, $I_n$ can be understood as injection terms in that they also appear, with the opposite signs, in the equation for the primary wave,

(2.36)\begin{equation} \dot{\mathcal{E}}_{\text{pri}}= - \sum_n \sum_{\sigma=\pm1}I_n^\sigma+\text{c.c.}, \end{equation}

which follows from (2.20). Note that the energy in each subchannel $\sigma$ is separately conserved, i.e. both energy and cross-helicity are conserved.

Equations (2.28) and (2.29) are similar to the equations that describe a one-dimensional chain of coupled linear oscillators. The $n$th elementary cell of the chain consists of two different oscillators characterized by $\psi _n^+$ and $\psi _n^-$, respectively. Similar equations emerge in studies of phase-space turbulence for the coupling between Hermite moments of the particle phase-space probability distribution (Hammett et al. Reference Hammett, Beer, Dorland, Cowley and Smith1993; Kanekar et al. Reference Kanekar, Schekochihin, Dorland and Loureiro2015; Parker Reference Parker2016; Adkins & Schekochihin Reference Adkins and Schekochihin2018; Nastac et al. Reference Nastac, Ewart, Sengupta, Schekochihin, Barnes and Dorland2023). Throughout this work, we use (2.28) and (2.29) to study the nature of collective oscillations in the chain of modulational harmonics, which can be understood as Floquet modes of the linearized system. Unstable oscillations of this system (i.e. MIs) lead to structure formation on top of the primary structure.

One can also consider (2.29) as a Schrödinger equation with a Hamiltonian $\mathrm {i}\boldsymbol {M}$. This Hamiltonian is not Hermitian, because modulations are parametrically coupled with the primary mode (through which energy can be either gained or lost) and dissipate at $n\to \infty$. At the same time, this Hamiltonian is invariant under the time-reversal transformation ($\mathrm {i}\to -\mathrm {i}$) and the parity transformation in the spectral space ($n\rightarrow -n$, $\sigma \rightarrow -\sigma$). This makes the system (2.28) $\mathcal {PT}$-symmetric (Bender Reference Bender2005; Bender et al. Reference Bender, Dorey, Dunning, Fring, Hook, Jones, Kuzhel, Lévai and Tateo2019). Depending on the balance of sources and sinks, such systems can support modes with entirely real frequencies (unbroken $\mathcal {PT}$ symmetry) and pairs of modes whose frequencies are mutually complex-conjugate (broken $\mathcal {PT}$ symmetry). This will be discussed further in § 4.

As long as the underlying ordering assumption [$|w_{\boldsymbol {q}+n\boldsymbol {p}}^\pm |/|w_{\boldsymbol {p}}^\pm | \sim O(\epsilon )$] holds, XQLT (2.28) provides an excellent approximation of the modulational dynamics of (2.3) (figure 2). However, practical applications of XQLT require truncations of this system. In the remainder of this work, we discuss to what extent such truncations can be adequate and compare our analytical results with XQLS (numerical simulations of (2.28) with the spectral boundary sufficiently far such that it is never reached in the time simulated), in lieu of the full (2.3).

Figure 2.

Comparison of the time evolution of various $\psi _n^+$ as predicted by DNS of (2.3) (red dots) and XQLS (numerical simulations of XQLT) (2.28) (black curves) for two representative cases of (a,c,e) the MI ($\theta = {\rm \pi}/6$) and (b,d,f) stable modulational dynamics ($\theta = {\rm \pi}/3$). Nonlinear DNS are seeded with random noise, resulting in a broadband modulational dynamics, but the results shown are for a modulation corresponding to $r = 0.5$, and the XQLS are initialized accordingly. For the $\theta = {\rm \pi}/6$ case, the harmonics shown are (a) $n = 0$, (b) $n = 1$ and (c) $n = 2$. In this case, XQLT adequately approximates the nonlinear dynamics until approximately $\tau \sim 30$. For the $\theta = {\rm \pi}/3$ case, the harmonics shown are (d) $n = 0$, (e) $n = 25$ and (f) $n = 50$. In this case, XQLT adequately approximates the nonlinear dynamics indefinitely.

3.1 Basic equations

A common way to truncate a problem like (2.29) is to postulate that all $\psi _n$ with $|n| > N$ are negligible. One can also understand this as imposing reflective boundary conditions in the oscillator chain $\{\boldsymbol {\psi }_n^\sigma \}$,

(3.1)\begin{equation} \boldsymbol{\psi}_{N+1}^\sigma = \boldsymbol{\psi}_{-(N+1)}^\sigma = 0. \end{equation}

Retained in this case are harmonics with wavevectors $\boldsymbol {p}$ and $\boldsymbol {q} + n\boldsymbol {p}$, with $n = 0, \pm 1, \ldots, \pm N$. This means that the resulting system has $Q = 1 + (1 + 2N)$ degrees of freedom. We call this procedure $Q$-mode truncation ($Q$MT). The corresponding truncation of $\boldsymbol {\psi }$ will be denoted as $\boldsymbol {\psi }_{Q\text {MT}}$, and the corresponding truncation of $\boldsymbol {M}$ will be denoted as $\boldsymbol {M}_{Q\text {MT}}$, so (2.29) becomes

(3.2)\begin{equation} \dot{\boldsymbol{\psi}}_{Q\text{MT}} = \boldsymbol{M}_{Q\text{MT}}\boldsymbol{\psi}_{Q\text{MT}}. \end{equation}

Because of (3.1), eigenmodes of this system can be understood as standing waves in the oscillator chain $\{\boldsymbol {\psi }_n^\sigma \}$ and have the form $\boldsymbol {\psi }_{Q\text {MT}} = \boldsymbol {Y}\exp (-\mathrm {i}\omega \tau )$, where $\omega$ are global frequencies and $\boldsymbol {Y}$ are constant ‘polarization vectors’. According to (2.29), these vectors satisfy $\boldsymbol {D}\boldsymbol {Y} = 0$, where $\boldsymbol {D} = \mathrm {i} \omega \mathbb {1} + \boldsymbol {M}_{Q\text {MT}}$ and $\mathbb {1}$ is a unit $Q\times Q$ matrix. Then, $\omega$ can be found from

(3.3)\begin{equation} \det \boldsymbol{D} = 0 \end{equation}

and MIs correspond to modes with $\operatorname {Im} \omega > 0$.

The success of this approach hinges on the assumption that higher harmonics truly remain negligible. In other words, a $Q$MT should be able to adequately describe the evolution of a modulational mode if its eigenvector $\boldsymbol {\psi }_{Q\text {MT}}$ is heavily weighted towards the low-order harmonics. Below, we explore to what extent this assumption is satisfied in various regimes.

3.2 Growth rates

The lowest-order $Q$MT is 4MT, which corresponds to $N = 1$. Within this model, the primary harmonic with $\boldsymbol {k} = \boldsymbol {p}$ is assumed to interact with the modulational mode at $\boldsymbol {k} = \boldsymbol {q}$ only through two sidebands $\boldsymbol {k} = \boldsymbol {q} \pm \boldsymbol {p}$; then, $\boldsymbol {\psi }_\text {4MT} = (\boldsymbol {\psi }_{-1}, \boldsymbol {\psi }_0, \boldsymbol {\psi }_1)^\intercal$ and

(3.4)\begin{equation} \boldsymbol{M}_{\text{4MT}} = \begin{pmatrix} 0 & \boldsymbol{G}_{ - 1}^+ & 0 \\ \boldsymbol{G}_0^- & 0 & \boldsymbol{G}_0^+\\ 0 & \boldsymbol{G}_1^- & 0 \end{pmatrix}. \end{equation}

(This model can be understood as quasilinear, because, although the modulational dynamics remains nonlinear, the second- and higher-order harmonics of the perturbation are neglected.) The 4MT equations yield

(3.5)\begin{equation} \omega^2 = \frac{r^4}{1+r^2}(1 \pm \sqrt{1 - (1 - r^{ - 4} \cos^22\theta) \sin^2 2\phi}). \end{equation}

(Note that these are normalized frequencies corresponding to the normalized time $\tau$. To obtain the actual physical frequencies, one must multiply the right-hand side of (3.5) by $\mathcal {A}$.) Assuming $\phi = {\rm \pi}/4$, one finds that two out of four branches of this dispersion relation are unstable at $|r| < 1$. The corresponding growth rate $\varGamma \doteq \operatorname {Im}\omega$ is maximized at $\theta = 0$ and $\theta = {\rm \pi}/2$ at the value

(3.6)\begin{equation} \varGamma = r \sqrt{\frac{1 - r^2}{1 + r^2}} \end{equation}

(from now on, we assume $r > 0$ for clarity of notation), and the corresponding polarization vectors are as follows:

(3.7)\begin{equation} \boldsymbol{Y}_0 = \begin{cases} (1, +1)^\intercal, & \theta = 0,\\ (1, -1)^\intercal, & \theta = {\rm \pi}/2. \end{cases} \end{equation}

The case $\theta = 0$ can be recognized as a purely hydrodynamic Kelvin–Helmholtz instability, i.e. one that does not involve magnetic field. The case $0 < \theta < {\rm \pi}/4$ can be understood as the MHD generalization of the Kelvin–Helmholtz instability. Note that the case $\theta = {\rm \pi}/2$ has the form of the tearing instability (Boldyrev & Loureiro Reference Boldyrev and Loureiro2018) (and similarly we may understand ${\rm \pi} /4 < \theta < {\rm \pi}/2$ as generalizations with shear flow), although the tearing instability should be stable in the ideal limit. However, the stability problem that we are studying here is somewhat different from the stability problem of tearing eigenmodes in that we also address transient dynamics, as will become clearer in § 4.

Potentially more accurate is a model with $N = 2$, or 6MT, which accounts for harmonics with $\boldsymbol {k} = \boldsymbol {q} \pm 2\boldsymbol {p}$; then, $\boldsymbol {\psi }_\text {4MT} = (\boldsymbol {\psi }_{-2}, \boldsymbol {\psi }_{-1}, \boldsymbol {\psi }_0, \boldsymbol {\psi }_1, \boldsymbol {\psi }_{2})^\intercal$ and

(3.8)\begin{equation} \boldsymbol{M}_{\text{6MT}} = \begin{pmatrix} 0 & \boldsymbol{G}_{ - 2}^+ & 0 & 0 & 0\\ \boldsymbol{G}_{ - 1}^- & 0 & \boldsymbol{G}_{ - 1}^+ & 0 & 0\\ 0 & \boldsymbol{G}_0^- & 0 & \boldsymbol{G}_0^+ & 0\\ 0 & 0 & \boldsymbol{G}_1^- & 0 & \boldsymbol{G}_1^+\\ 0 & 0 & 0 & \boldsymbol{G}_2^- & 0 \end{pmatrix}. \end{equation}

One can derive an analytical expression for $\omega$ from here just like we derived (3.5) from (3.4). However, this expression is cumbersome and not particularly instructive, so it is not presented explicitly. Truncations with larger $Q$, albeit also explored by us, will not be presented either for the same reason.

Figure 3 shows a comparison of the 4MT and 6MT predictions and the MI growth rates inferred numerically from XQLS. Note that the truncated models quantitatively agree with each other and with XQLS at some parameters but drastically disagree at other parameters. In particular, the 4MT predicts that $\varGamma$ is an even function of ${\rm \pi} /2 - \theta$, while the 6MT predicts that this is not the case. Furthermore, XQLT predicts that, for example, at $r = 0.5$, the system is stable at all $\theta > {\rm \pi}/2$ (figure 3b). Let us discuss this in detail.

Figure 3.

(a) The growth rate $\varGamma$ (of the most unstable mode) versus the normalized wavenumber $r$ at $\theta = 0$: 4MT (blue); 6MT (red); XQLS (black markers). (b) The same for $\varGamma$ versus $\theta$ at $r = 0.5$. The 4MT predicts identical growth rates at $\theta = 0$ and $\theta = {\rm \pi}$, while the 6MT and XQLS predict that the system is unstable at $\theta = 0$ and stable at $\theta = {\rm \pi}$. (c) The same for $\varGamma$ versus $\phi$ at $r = 0.5$ and $\theta =0$.

As expected from the analytical formula (3.6) based on the 4MT, and as we have also found numerically, MI is typically suppressed at $r \gtrsim 1$; hence, below we focus on the regime $r \in (0,1)$. A comparison of the growth rate predicted by XQLT and truncated models in this range is presented figure 4. The corresponding primary modes with $\theta \lesssim {\rm \pi}/4$ are modulationally unstable, and the growth rates predicted by 4MT and 6MT are in good agreement with XQLT. Such primary structures are dominated by the velocity field, $\mathcal {E}_{\text {pri}, v} > \mathcal {E}_{\text {pri}, b}$ (see (2.27)), and thus will be called $\boldsymbol {v}$-dominated primary modes (VDPMs). At $\theta \gtrsim {\rm \pi}/4$, primary structures are dominated by the magnetic field, $\mathcal {E}_{\text {pri}, b} > \mathcal {E}_{\text {pri}, v}$, and thus will be called $\boldsymbol {b}$-dominated primary modes (BDPMs). As seen in figure 4, BDPMs are modulationally stable, but truncated models predict otherwise. Figure 4(a,b) shows the discrepancy between the predictions of XQLT, 4MT and 6MT for representative parameters, specifically, $\theta \in (0, {\rm \pi}/2)$ and $\phi = {\rm \pi}/4$. (Stability is primarily determined by $\theta$ and $r$. Deviations of $\phi$ from ${\rm \pi} /4$ result only in quantitative adjustments, unless $\phi$ is close to $0$ or ${\rm \pi} /2$.) In other words, truncated modes tend to overestimate the growth rate and systematically produce false positives for instability. Figure 5 shows that, although the magnitude of disagreement decreases with $N$, the systematic overprediction of instability persists.

Figure 4.

The difference between the growth rate predicted by a $Q$MT and that inferred from XQLS, $\varGamma _{Q\text {MT}}-\varGamma _{\text {XQLS}}$, versus $(\theta, r)$ for $\phi = {\rm \pi}/4$: (a) 4MT; (b) 6MT. The dashed curves mark the stability boundary as determined by a parameter scan with XQLS. The preponderance of the red colour, especially outside the instability domain, indicates that truncated models tend to produce false positives for instability.

Figure 5.

The growth rates $\varGamma$ obtained through a $Q$MT for $\theta = 0$ (blue markers) and $\theta = {\rm \pi}/2$ (orange markers) versus the number of modes retained, $N$. The solid-coloured lines indicate the corresponding XQLS growth rates, to which the rates predicted by the truncated models converge at large $N$.

This can be understood as follows. Consider two possible scenarios for modulational stability. The first possibility is that there is no tendency for the primary mode to inject energy into such a modulation in the first place. (Injection is localized to small $|n|$, so this corresponds to regimes where both low-order $Q$MTs and XQLT predict stability.) Alternatively, the primary mode tends to inject energy but the energy escapes to higher harmonics instead of entering positive feedback loops at low $|n|$. (This corresponds to the false positive case where $Q$MTs predict instability when the system is actually modulationally stable.) This overprediction can therefore be understood as a consequence of simple truncations constituting perfectly reflecting boundary conditions for the energyFootnote ⁵ , i.e. the outward energy flux at the cutoff harmonic $n=N$ is enforced to be zero if $|\psi _{N+1}^\pm |$ are held to be zero (2.35). By trapping energy that might otherwise escape to higher $|n|$ at lower harmonics, $Q$MTs can generate false positive feedback loops with injection (2.34). True MIs then correspond to modulational dynamics that naturally remain confined to small $|n|$.

This can be seen by analysing the eigenmode structure. (See also Appendix A for an alternative explanation.) As to be discussed in § 4, unstable modes are evanescent spectral waves, whose $|Y_n|$ exponentially decrease with $|n|$ (figure 6a),

(3.9)\begin{equation} |Y_n| \propto \exp(-|n|/l). \end{equation}

(We assume $\phi = {\rm \pi}/4$ for simplicity, so $|Y_n^+| = |Y_n^-|$ and the upper index in $|Y_n^\sigma |$ can be omitted.) That is why truncated models correctly predict MI for modes that are actually unstable. The spectral scale $l$ depends on the system parameters, and the smaller $l$ the more accurate $Q$MT is. At some parameters, though, particularly when $\theta \rightarrow {\rm \pi}/4$, $l$ becomes large or even infinite, such that $|Y_n|$ asymptotes to a non-zero constant at $|n| \to \infty$ (figure 6b). In this case, $Q$MT is bound to fail. The corresponding modes are PSWs. While they also receive energy from the primary structure at low $n$, PSWs transport this energy along the spectrum to $|n| \to \infty$. There, the energy is unavoidably dissipated by viscosity or resistivity, however small those are. Below, we discuss these effects in more detail.

Figure 6.

Typical structures of eigenmodes, specifically, $|Y_n|$ versus $n$, at $r=0.5$ and $\phi ={\rm \pi} /4$ for various $\theta$: (a) unstable modes supported by VDPMs; (b) oscillatory solutions supported by BDPMs. (The upper index in $|Y_n^\sigma |$ is omitted because $|Y_n^+| = |Y_n^-|$ for the assumed parameters.) Here $|Y_n|$ decreases exponentially with $n$ for the former (notice the logarithmic scale) but asymptote to non-zero constants for the latter. The asymptotes are indicated by the dotted lines. These results are obtained through XQLS using the initial conditions of the form (4.6) and normalized such that $|Y_0| = 1$.

4.1 Dynamics at large $|n|$

First, let us consider spectral waves at $|n| \gg 1$. In this limit, one has $\alpha _n^\pm \rightarrow \mp r$ and $\beta _n^\pm \rightarrow 0$, so (2.29) yields two decoupled wave equations for $\psi _n^+$ and $\psi _n^-$,

(4.1a)\begin{gather} \dot{\psi}_n^- = r\cos\phi (\mathrm{e}^{\mathrm{i}\theta} \psi_{n-1}^- - \mathrm{e}^{-\mathrm{i}\theta}\psi_{n+1}^-), \end{gather}

(4.1b)\begin{gather}\dot{\psi}_n^+ = r\sin\phi (\mathrm{e}^{-\mathrm{i}\theta} \psi_{n-1}^+ - \mathrm{e}^{\mathrm{i}\theta}\psi_{n+1}^+). \end{gather}

These equations allow solutions in the form of monochromatic waves,

(4.2)\begin{equation} \psi_n^\sigma = \varPsi^\sigma \exp(- \mathrm{i}\omega \tau + \mathrm{i} K^\sigma n), \end{equation}

with constant amplitude $\varPsi ^\sigma$, frequency $\omega$ and ‘wavenumber’ $K^\sigma$. The quotation marks are added as a reminder that the corresponding waves propagate in the spectral space rather than in the physical space. Also note that $K^\pm$ are defined unambiguously only within the first Brillouin zone, $K^\pm \in (-{\rm \pi}, {\rm \pi})$.

From (4.1), one readily finds that spectral waves obey the following dispersion relations:

(4.3a,b)\begin{equation} \omega = \omega_0^\sigma \sin \kappa^\sigma,\quad \kappa^\pm \doteq K^\pm{\pm}\theta, \end{equation}

where $\omega _0^+ = 2r\sin \phi$ and $\omega _0^- = 2r\cos \phi$. Accordingly, spectral waves in the $n$ space along the $\psi ^\sigma$ channel have the group velocity

(4.4)\begin{equation} v_g^\sigma = \omega_0^\sigma \cos \kappa^\sigma \end{equation}

and can propagate either up or down the spectrum (figure 7). The dependence of $\omega$ and $v_g^\pm$ on $\theta$ can be removed by a gauge transformation. Specifically, $\kappa ^\pm$ becomes the true wavenumber if one adopts the variables $\bar {\psi }_n^\pm \doteq \psi _n^\pm \exp ({\mp \mathrm {i} n\theta })$ instead of $\psi _n^\pm$. Also notably, to the extent that the (large) index $n$ can be approximately considered a continuous variable, (4.1) can be written as usual non-dispersive-wave equations for $\zeta ^\pm (\tau, n) \doteq \bar {\psi }_n^\pm (\tau )$,

(4.5)\begin{equation} \partial_\tau \zeta^\pm = - \omega_0^\pm \partial_n \zeta^\pm. \end{equation}

This model corresponds to the small-$\kappa ^\sigma$ limit of (4.3a,b), that is, $\omega \approx \omega _0^\sigma \kappa ^\sigma$.

Figure 7.

Numerical simulations showing PSW packets propagating (a) down the spectrum ($|n|$ of the packet centre increases with time) and (b) up the spectrum ($|n|$ of the packet centre decreases with time). In both cases, $r = 0.5$ and $\theta = {\rm \pi}/3$. Both packets are initialized using $\psi _n^+(\tau = 0) = \exp [-(n - n_0)^2/\varsigma + \mathrm {i} K^+ n]$ with $n_0 = 200$, $\varsigma = 150$, and (a) $K^+ = -{\rm \pi} /6$, (b) $K^+ = {\rm \pi}/2$.

The above equations can be used to describe PSW packets localized at large $|n|$ but not the global eigenmodes (because the latter involve dynamics at small $|n|$). In particular, (4.3a,b) are not enough to find the global-mode frequencies. Yet they allow one to determine the mode structure at $|n| \gg 1$ if one knows the mode frequency from other considerations (or, vice versa, to find the mode frequency if $\kappa ^\pm$ are known). In particular, unstable modes have complex $\omega$, so, according to (4.3a,b), their asymptotic wavenumber cannot be real. This explains why unstable modes are evanescent. In contrast, for real $\omega$, (4.3a,b) allows for real wavenumbers (provided that $|\omega | \le \omega _0^\pm$), i.e. predicts PSWs. These results are corroborated by XQLS. For example, figure 6 shows the mode structures inferred from the asymptotic dynamics of $\psi _{n}^\sigma (\tau )$ at large enough $\tau$ for the initial conditions

(4.6)\begin{equation} \psi_{n}^\sigma(\tau=0) = \begin{cases} \epsilon \exp(\mathrm{i}\xi^\sigma_n), & n ={\pm} 1,\\ 0, & n \neq{\pm} 1. \end{cases} \end{equation}

Here, $\epsilon$ is a constant small amplitudeFootnote ⁶ , and $\xi _n$ are parameters that change the polarization of the initial conditions while keeping the initial modulation energy fixed. (The specific values of $\xi _n^\sigma$ are given in the captions of figures in which initial condition dependent results of XQLS are presented.) This form of the initial conditions is chosen to emulate a simple modulation of the primary mode.

The fact that the properties of spectral waves are determined only by the asymptotic properties of $\alpha _n^\pm$ and $\beta _n^\pm$ makes these waves particularly robust. In particular, note that the asymptotic form of $\alpha _n^\pm$ and $\beta _n^\pm$ at $|n| \to \infty$ remains the same even when $\boldsymbol {p}$ and $\boldsymbol {q}$ are not orthogonal, so the above results readily extend to general $\boldsymbol {q}$. Likewise, a similar picture holds also for three-dimensional interactions, although the mode polarization can be more complicated in this case. Notably, PSWs provide ballistic rather than diffusive (or superdiffusive) energy transport along the spectrum. This distinguishes them from the cascades that are commonly discussed in turbulence theories and result from eddy–eddy interactions, or wave–wave collisions (Galtier et al. Reference Galtier, Nazarenko, Newell and Pouquet2000).

4.2 Global modes in the form of PSWs

In realistic settings, spectral waves are excited at finite $|n|$ and propagate down the spectrum. Evanescent waves reach $|n| \sim l$ (see (3.9)) in finite time $\tau _0$, get reflected, and eventually (asymptotically) settle as stationary standing waves on time scales of the order of a few $\tau _0$. Such waves can be adequately described by truncated models. In contrast, PSWs continue to propagate towards infinity indefinitely (until they dissipate at scales where viscosity or resistivity is no longer negligible). Thus, they never develop components propagating up the spectrum at large $|n|$ and so they never become quite like the standing-wave eigenmodes predicted by truncated models (§ 3.1). This means that actual PSWs cannot be adequately described by analysing $Q$MT in principle. Instead, they should be considered in the context of the initial-value problem, much like Langmuir waves appear in plasma kinetic theory within Landau's initial-value approach as opposed to the Case–van Kampen true-eigenmode approach (Stix Reference Stix1992).

Consequently, even though (2.29) has infinitely many degrees of freedom, the number of relevant global modes, i.e. those that propagate towards infinity, is finite. We find that there are only two of them, which can be attributed to the fact that $\dim \boldsymbol {\psi }_n = 2$. (Likewise, at a given wavenumber, only one, Langmuir, branch of relevant electrostatic waves in electron plasma remains relevant after transients are gone.) For example, for $\phi = {\rm \pi}/4$, we find from XQLS that

(4.7)\begin{equation} \omega ={\pm} \sqrt{2}r \cos \theta. \end{equation}

This dependence and the corresponding mode structures are illustrated in figure 8. Each of the two modes has two wavenumbers associated with it: one determines $\psi ^\sigma _n$ at $n \to +\infty$, and the other one determines $\psi ^{-\sigma }_n$ at $n \to -\infty$ (figure 9). Thus, there are four $\kappa ^\sigma _d$ overall, where $\sigma$ denotes the corresponding component of $\boldsymbol {\psi }_n$, and $d$ denotes the direction of propagation ($d = \pm 1$ is the sign of the group velocity at $|n| \to \infty$). They can be found by combining (4.7) with (4.3a,b) and applying $\operatorname {sgn} v_g = d$, i.e. $\omega _0^\sigma d\cos \kappa ^\sigma _d > 0$. This yields

(4.8a)\begin{gather} \omega > 0: \quad \kappa^+_+ = \frac{\rm \pi}{2}-\theta,\quad \kappa^-_- = -\frac{3{\rm \pi}}{2}+\theta, \end{gather}

(4.8b)\begin{gather}\omega < 0: \quad \kappa^+_- = \frac{3{\rm \pi}}{2}-\theta,\quad \kappa^-_+ = -\frac{\rm \pi}{2}+\theta, \end{gather}

where $\kappa ^\sigma _d=K^\sigma _d+\sigma \theta$ and the values of $K^\sigma _d$ are defined modulo the Brillouin zone (i.e. $K^-_-=-3{\rm \pi} /2$ should be understood as ${\rm \pi} /2$ and so on). Remember, though, that this applies only to global eigenmodes. Transient waves propagating in the region $|n| \gg 1$ can have any $\kappa ^\sigma$ (figure 7).

Figure 8.

(a) The dependence of the global-mode frequency $\omega > 0$ on $\theta$ for various $r$. The dashed lines are the inferred solutions (4.7). (b) Here $|Y_n^+|$ (blue) and $|Y_n^-|$ (orange) for the mode with $\omega > 0$. The eigenmode amplitudes for the $\omega < 0$ mode are identical with the roles of $|Y_n^+|$ and $|Y_n^-|$ switched.

Figure 9.

The XQLS showing global-mode PSWs for $r=0.5$, $\theta = {\rm \pi}/3$ and $\phi = {\rm \pi}/4$: (a) $\operatorname {Re} \psi _n^+/\epsilon$ for a PSW seeded by the initial conditions (4.6) with $\xi ^\sigma _d = d\exp (\sigma \mathrm {i} {\rm \pi}/4)$ (colour bar). The dashed lines indicate the fronts propagating at (i) the maximum spectral speed $\sqrt {2}r$ and (ii) the actual group velocity of the mode. The field between these dashed lines consists of transients, which are negligible behind the second front. (b) The total normalized energy of the modulation, $\mathcal {E}_{\text {mod}}/\epsilon ^2$, versus $\tau$ (black), along with its kinetic (red) and magnetic (blue) components. (c) The profiles of the spectral energy density at $\tau =100, 150, 200$, clearly exhibiting left- and right-propagating fronts. The horizontal dashed line indicates the average spectral energy density $\overline {\mathcal {E}_n}$, where the average is taken over the PSW period and over $n$ between the energy fronts.

The fact that the frequencies of these modes are real indicates that the modes are self-sustained notwithstanding the dissipation that they unavoidably experience at $|n| \to \infty$. The corresponding evolution of the wave spectrum is illustrated in figure 9(a) and the corresponding evolution in the physical space is shown in figure 10. This remarkable dynamics is understood from the fact that the energy drain via this dissipation is balanced by the energy injection from the primary mode at small $|n|$ (see figure 9c).

Figure 10.

The global-mode PSW mode with $\omega > 0$, $r = 0.5$ and $\theta = {\rm \pi}/3$; the same as figure 9, but in the real space as opposed to the spectral space. Specifically shown are (a) $z_y^+$ as a function of $(\tau, x)$ at $y = 0$; (b) $z_x^+$ as a function of $(\tau, y)$ at $x = 0$.

The linear spectral waves discussed so far correspond to the regime when the change of $\mathcal {E}_{\text {pri}}$ is negligible. Eventually, though, as higher harmonics are populated and thus absorb more energy, the primary wave will be depleted. This means, in particular, that no primary wave can be truly stationary. Even an arbitrarily small perturbation will generally launch a PSW, and this PSW will eventually deplete the primary wave.

4.3 Anomalous dissipation

Since a PSW mode carries energy towards $|n| \rightarrow \infty$, eventually, a large enough $|n|$ is reached where the energy is dissipated regardless of how small the viscosity and resistivity are. Thus, a PSW exhibits an effective, or ‘anomalous’, dissipation rate $\gamma$ that is independent of $\nu$ and $\eta$ in the limit $\nu, \eta \to 0$. This effect is different from the anomalous transport caused by eddy–eddy collisions in turbulence (for example, see Donzis, Sreenivasan & Yeung (Reference Donzis, Sreenivasan and Yeung2005)) in that the energy transport caused by a PSW is ballistic. As a result, $\gamma$ is straightforward to calculate, which is done as follows.

As discussed in § 4.2, a global PSW at large $|n|$ has a flat mode structure, i.e. a structure with $\mathcal {E}_n$ independent of $n$. This structure establishes itself as an expanding ‘shelf’ whose edges (wavefronts) propagate across the modulational spectrum to $n \to \pm \infty$ at the group velocity $v_g$ (figure 11). Since the height of the shelf, $\overline {\mathcal {E}_n}$ (the overbar here denotes averaging over the PSW temporal period and over all $n$ within the shelf), remains constant, this process drains energy from the primary mode linearly in time,

(4.9)\begin{equation} \dot{\mathcal{E}}_{\text{pri}} = - 2|v_g|\overline{\mathcal{E}_n} = \text{const.}, \end{equation}

where we consider the dynamics average over the PSW temporal period. Hence,

(4.10)\begin{equation} \gamma \doteq - \frac{\dot{\mathcal{E}}_{\text{pri}}}{\mathcal{E}_{\text{pri}}} = 2|v_g|\frac{\overline{\mathcal{E}_n}}{\mathcal{E}_{\text{pri}}} \sim \omega_0^\pm \epsilon^2. \end{equation}

(The rate relative to the physical time $t$, as opposed to $\tau$, is $\mathcal {A}$ times this $\gamma$.)

Figure 11.

(a) The magnitude of the group velocity $|v_g|$ of the global PSW mode (i.e. the speed at which the energy front propagates along the spectrum) versus $\theta$ for various $r$. (b) The same for the average spectral energy density $\overline {\mathcal {E}_n}/\epsilon ^2$. (c) The same for the resulting average drain rate $\gamma /\epsilon ^2$, where $\gamma$ is defined in (4.10). The initial conditions used throughout these figures are given by (4.6), with $\xi _{1}^\sigma =-\xi _{-1}^\sigma =\exp (\sigma \mathrm {i} {\rm \pi}/2)$.

Let us again assume the initial conditions (4.6) and $\phi = {\rm \pi}/4$. Through XQLS, we find that the global PSW-mode amplitude is maximized when $|\arg (\xi ^\sigma _{\pm 1}/\xi ^{-\sigma }_{\pm 1})| = |\arg (\xi ^\sigma _{\pm 1} / \xi ^\sigma _{\mp 1})| = {\rm \pi}$, irrespective of $\theta$. (Any particular choice of $\xi ^\sigma _n$ that satisfies this condition affects only phases of the modulational dynamics and does not impact averaged quantities that determine $\gamma$.) This corresponds to the case when all modulational-seed energy is in the magnetic field; i.e. $\mathcal {E}_b = \mathcal {E}$ and $\mathcal {E}_v = 0$. Figure 11 shows the corresponding anomalous dissipation rate normalized to the seed energy, along with its determining factors – the system-dependent PSW group velocity $v_g$ and the average spectral energy density $\overline {\mathcal {E}_n}$, which is determined by the initial conditions. It can be seen that, for the MI unstable VDPMs ($\theta < {\rm \pi}/4$), spectral waves are relatively slow and have a small amplitude. The transition to modulational stability occurs at $|v_g| \sim \varGamma _{4MT}$, where $\varGamma _{4\mathrm {MT}}$ is the MI growth rate predicted by $4$MT. This condition can be understood as a threshold beyond which (i.e. at $|v_g| \gtrsim \varGamma$) PSWs provide a sufficiently fast escape route for the energy injected at small $|n|$, such that the positive feedback loops (see (2.34)) supporting MIs can no longer be sustained. In the context of $\mathcal {PT}$ symmetry (§ 2.2.2), the spectral group velocity can be understood as the effective coupling between the links in the oscillator chain, connecting the energy injection from the primary mode at small $n$ to the energy sink at $n\to \infty$. As $|v_g|$ increases with $\theta$, the system transitions from broken to unbroken $\mathcal {PT}$ symmetry.

Also, notice the following. If modulations at multiple $\boldsymbol {q}$s are present, they launch independent cascades along $\boldsymbol {k} = \boldsymbol {q} + n\boldsymbol {p}$ and thus the drain on the shared primary mode will be additive. This regime is illustrated by figure 1(d–f), which show the modulational dynamics for a BDPM ($\theta = {\rm \pi}/3$) seeded with noise. The rather nondescript modulational dynamics seen in figure 1(d,e) can be understood as the superposition of the many structures like that in figure 10 for various $\boldsymbol {q}$. Also, the linear-in-time depletion of the primary mode seen in figure 1(f) can be understood as the sum of the PSW-driven energy fluxes along the constituent modulational channels.