On reduced modelling of the modulational dynamics in magnetohydrodynamics

Abstract

This paper explores structure formation in two-dimensional magnetohydrodynamic (MHD) turbulence as a modulational instability (MI) of turbulent fluctuations. We focus on the early stages of structure formation and consider simple backgrounds that allow for a tractable model of the MI while retaining the full chain of modulational harmonics. This approach allows us to systematically examine the validity of popular closures such as the quasilinear approximation and other low-order truncations. We find that, although such simple closures can provide quantitatively accurate approximations of the MI growth rates in some regimes, they can fail to capture the modulational dynamics in adjacent regimes even qualitatively, falsely predicting MI when the system is actually stable. We find that this discrepancy is due to the excitation of propagating spectral waves (PSWs) which can ballistically transport energy along the modulational spectrum, unimpeded until dissipative scales, thereby breaking the feedback loops that would otherwise sustain MIs. The PSWs can be self-maintained as global modes with real frequencies and drain energy from the primary structure at a constant rate until the primary structure is depleted. To describe these waves within a reduced model, we propose an approximate spectral closure that captures them and MIs on the same footing. We also find that introducing corrections to ideal MHD, conservative or dissipative, can suppress PSWs and reinstate the accuracy of the quasilinear approximation. In this sense, ideal MHD is a ‘singular’ system that is particularly sensitive to the accuracy of the closure within mean-field models.

1 Introduction

1.1 Motivation

Coherent-structure formation from turbulence is ubiquitous in nature, intrinsically compelling and one of the precious few aspects of turbulence that are relatively yielding to analytical efforts (Hussain 1986; Smolyakov, Diamond & Malkov 2000; Krashennikov, D'Ippolito & Myra 2008; Karimabadi et al. 2013). For magnetohydrodynamic (MHD) turbulence (Biskamp 2003; Beresnyak 2019; Schekochihin 2022) and the turbulent-dynamo problem (Brandenburg, Sokoloff & Subramanian 2012; Rincon et al. 2016; Tobias 2021), this has spawned a vast body of work known as mean-field electrodynamics, or mean-field dynamo theory (Moffatt 1978; Krause & Rädler 1980; Gruzinov & Diamond 1994; Rädler 2007; Brandenburg 2018). In such approaches, the velocity and magnetic fields are split into fluctuations and mean, or average, fields (the definition of average depends on the problem); then, various closures are implemented to obtain the mean-field evolution in response to the fluctuations (Pouquet, Frisch & Léorat 1976; Kraichnan 1977; Nicklaus & Stix 1988; Blackman & Field 2002; Brandenburg & Subramanian 2005b).

These closures ultimately rely on a truncation of a formally infinite chain of correlation functions of increasingly high order.¹ This is a common approach to reduced modelling of nonlinear systems, and one can expect it to work when the fluctuations are small compared with the mean field. However, in turbulent structure formation, the usual roles of perturbation and background are switched, that is, mean fields emerge as small modulations of order-one correlation functions that characterize turbulent fluctuations. Hence, the question of the validity of low-order closures becomes non-trivial.

In this paper, we are concerned with the popular closure called the first-order smoothing (Krause & Rädler 1980) or second-order correlation approximation (Rädler 1982). It implies that the mean-field evolution is determined only by second-order correlations of the fluctuations, which, in a broader context, is also known as the quasilinear approximation (Dodin 2022). On the one hand, one can expect its validity domain to be extremely limited. It can be rigorously justified only for small magnetic Reynolds numbers or short correlation times (Brandenburg & Subramanian 2005a; Rincon 2019), and comparisons with direct numerical simulations (DNS) often corroborate this expectation (Schrinner et al. 2005; Pipin & Proctor 2008). On the other hand, the quasilinear approximation endures in its popularity as a workhorse that is too convenient to cast aside in analytical calculations (Masada & Sano 2014; Squire & Bhattacharjee 2015; Lingam & Bhattacharjee 2016; Gopalakrishnan & Singh 2023). Furthermore, in some cases, quasilinear calculations have been found to produce good agreement with DNS even outside of their formal validity domain (Käpylä et al. 2006).² It is therefore important to understand when, and how exactly, the quasilinear approximation fails. Answering this question requires exploring a mean-field model that retains high-order correlations. That is what the present paper is about.

1.2 Outline

Here, we explore structure formation in MHD turbulence as a modulational instability (MI) (Zakharov & Ostrovsky 2009) of the flow-velocity and magnetic-field fluctuations comprising the background turbulence. (For a recent overview of this approach in application to drift-wave and hydrodynamic turbulence, see Tsiolis, Zhou & Dodin (2020) and Zhu & Dodin (2021).) For simplicity, we restrict our consideration to two-dimensional (2-D) turbulence. Although bounded planar flows cannot sustain long-lasting magnetic fields in two dimensions (an antidynamo theorem by Zeldovich (1957)), 2-D turbulence nevertheless hosts rich physics and has long served as a popular testing ground for various aspects of MHD turbulence theory due to its (relative) tractability (Pouquet 1978; Biskamp, Schwarz & Drake 1996; Biskamp & Schwarz 2001; Nazarenko 2007; Giorgio, Servidio & Veltri 2017). The goal of this paper is to complement those works with a detailed study of MIs using both reduced models and DNS.

Specifically, we study how energy is transferred from primary fluctuations to mean magnetic field and mean velocity during the early stages of structure formation, when these mean fields are weak and the primary structure can be treated as fixed. To analyse this process in detail, we assume a maximally reduced model of the turbulent fluctuations, namely, a spatially monochromatic primary structure. The simplicity of this model has been detracting attention from it in the literature in favour of more realistic turbulence settings, but there is more to it that meets the eye. Even with a monochromatic primary structure, the modulational dynamics remains remarkably intricate and, depending on the system parameters, can change drastically (figure 1) in ways that are not captured by the common closures or obvious from the standard eigenmode analysis (Friedberg 1987). However, the simplicity of our model helps unravel this mystery.

Figure 1. The DNS of modulational dynamics of (2.3) seeded at $\tau =0$ with random noise. (This figure is discussed in detail in later sections. For the notation, see §§ 2.1 and 2.2.) Panels (a,c,e) correspond to a modulationally unstable primary mode (2.17) with $\theta ^\pm =\pm {\rm \pi}/6$, and (b,d,f) correspond to a modulationally stable primary mode with $\theta ^\pm =\pm {\rm \pi}/3$. Panels (a,b) show $z_x^+(t,y)$, (c,d) show $z_y^+(t,x)$ and (e,f) show the corresponding energy breakdown. Specifically, $\mathcal {E}_v$ is the normalized kinetic energy (2.21c), $\mathcal {E}_b = E_b/E_{\text {tot}}$ is the normalized magnetic energy (2.21d), $\mathcal {E}_{\text {mod}}$ is the normalized modulational-mode energy (2.21b) and $\mathcal {E}_{\text {pri}}$ is the normalized primary-mode energy (2.21a). The colour bar shows the field amplitudes normalized to $\mathcal {A}/p$. The spatial coordinates are normalized to $1/p$.

We find that, when the primary structure truly experiences a MI, this MI can be described well with typical truncated models that neglect high-order correlations. However, already in adjacent regimes, such truncated models can fail spectacularly and produce false positives for instability. To study this process in detail, we propose an ‘extended’ quasilinear theory (XQLT) that treats the primary structure as fixed but includes the entire spectrum of modulational harmonics (as opposed to just the low-order harmonics, as usual). From XQLT, also corroborated by DNS, we find that the difference between said regimes is due to a fundamental difference in the structure of the modulational spectrum. For unstable modes, the spectrum is exponentially localized at low harmonic numbers, so truncated models are justified. But this localization does not always occur. At other parameters, modulational modes turn into constant-amplitude waves propagating down the spectrum, unimpeded until dissipative scales. These spectral waves are self-maintained as global modes with real frequencies and cause ballistic energy transport along the spectrum, breaking the feedback loops that could otherwise sustain MI. The ballistic transport drains energy from the primary structure at a constant rate until the primary structure is depleted. Notably, these effects are overlooked if one assumes from the start that the modulation spectrum is confined and that, accordingly, dissipation entirely vanishes in the limit of zero viscosity and resistivity.

A comment is due here regarding how the established truncated MHD closures fit into this picture. Given the indispensability of such simple closures, a higher-resolution understanding of their limitations can help guide and interpret their inevitable application. For example, we show that both conservative and dissipative corrections to ideal MHD tend to restore the validity of the quasilinear approximation. In this sense, ideal incompressible MHD is a ‘singular’ theory that is particularly sensitive to the accuracy of the closure within mean-field models. That said, it is an important conclusion of our work that, unless deviations from ideal incompressible MHD are substantial, changing the turbulence parameters even slightly can produce qualitative inaccuracies in an otherwise workable reduced model. Therefore, efforts to extrapolate closures or infer them from general considerations should be supplemented with first-principle calculations based on the underlying modulational dynamics. As a step in this direction, we propose a spectral closure that captures both the propagating spectral waves (PSWs) and MIs within our model. (But, of course, more work remains to be done for generic MHD turbulence.)

The paper is organized as follows. In § 2, we present the basic equations and derive the XQLT. In § 3, we examine the ability of truncated models of low-order to predict the MI rate. In § 4, we introduce spectral waves and describe their properties. In § 5, we propose a closure that can simultaneously capture MIs and spectral waves. In § 6, we demonstrate the impact of corrections to ideal MHD on the modulational dynamics. In § 7, we summarize our main results. Auxiliary calculations are presented in appendices.

2 Extended quasilinear theory

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 1950), 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.³ 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.⁴ (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. 1993; Kanekar et al. 2015; Parker 2016; Adkins & Schekochihin 2018; Nastac et al. 2023). 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 2005; Bender et al. 2019). 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 Truncated models

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 2018) (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 energy⁵ , 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 Spectral waves

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 amplitude⁶ , 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. 2000).

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 1992).

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 (2005)) 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.

5 Unified closure

As discussed in § 3.1, naive $Q$MTs assume the perfectly reflecting boundary conditions (3.1) for spectral waves, thus precluding PSWs and ignoring potential dissipation at large $|n|$. One can expect that, instead of (3.1), it would be more accurate to adopt the following closure in anticipation of spectral waves:

(5.1)\begin{equation} \psi_{N+1}^\sigma = \mathrm{e}^{\mathrm{i} K^\sigma} \psi_N^\sigma \end{equation}

(and similarly for $\psi _{-(N+1)}^\sigma$). The value of $K^\sigma$ is found, for given $\omega$, from (4.3a,b),

(5.2)\begin{equation} \mathrm{e}^{\mathrm{i} K^\sigma} = \exp({-\mathrm{i}\sigma\theta}) [\mathrm{i}\varOmega^\sigma + \sqrt{1 - (\varOmega^\sigma)^2}], \end{equation}

where $\varOmega ^\sigma \doteq \omega /\omega _0^\sigma$, and we restrict our attention to the range $|\operatorname {Re}\varOmega |\leq 1$ to avoid discontinuities, a choice which is consistent with the range of observed solutions. The sign of the square root is chosen to enforce outward propagation at $\operatorname {Im}\omega =0$, which also enforces that unstable solutions exponentially decay in the direction of propagation $|\exp (\mathrm {i} K^\sigma _+)|<1$.⁷ (Conversely, damped solutions grow in the direction of propagation.) As seen from figure 12, the resulting closure

(5.3)\begin{equation} \psi_{N+1}^\sigma = \exp({-\mathrm{i}\sigma\theta}) [\mathrm{i}\varOmega^\sigma + \sqrt{1 - (\varOmega^\sigma)^2}]\psi_N^\sigma \end{equation}

leads to reasonably accurate results already $N = 2$, both for MIs and PSWs. In particular, note that $\operatorname {Im} \omega$ is substantially suppressed relative to the predictions of $Q$MTs, as the PSW closure captures effective dissipation to higher harmonics.⁸

Figure 12. The $\theta$-dependence of the global-mode frequencies derived from truncated models with the closure (5.3) for various $N$: (a–c) $\operatorname {Re}\omega$ and (d–f) $\operatorname {Im}\omega$. Here (a,d) $N = 2$; (b,e) $N = 3$; (c,f) $N = 4$. The red and blue curves indicate the branches that are the best match to the, respectively, unstable and stable solutions obtained through XQLS (black dashed lines). The grey curves show the remaining spurious solutions due to finite $N$. Note that MIs exist only for $\theta \in (0, {\rm \pi}/4)$, while PSWs exist in the entire range $\theta \in (0, {\rm \pi}/2)$.

Because the closure becomes exact only in the limit $N \rightarrow \infty$, it also yields spurious solutions (grey curves in figure 12) at finite $N$. One can understand this from the analogy with the Case–van Kampen modes mentioned in § 4.2. These solutions become increasingly stable, and thus less of an issue, as $N$ increases. Also note that the closure adequately describes true MIs. This is explained by the fact that $\psi _n^\sigma$ at large $|n|$ are exponentially small and do not matter for unstable modes anyway.

6 Quasilinear approximation beyond ideal MHD

As discussed in the previous sections, the existence of PSWs undermines the standard quasilinear approximation in application to ideal incompressible MHD. Given the ubiquity of quasilinear modelling in the literature, it may seem concerning that the quasilinear approximation can fail so spectacularly. Interestingly, though, the quasilinear approximation is somewhat more robust beyond the ideal-MHD limit, both due to conservative corrections and dissipation.

Let us discuss the former first. Without attempting to describe any particular physical system, let us consider the following modified version of (2.3):

(6.1)\begin{equation} \partial_t \boldsymbol{z}^\pm = - (\boldsymbol{z}^\mp\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{z}^\pm - \boldsymbol{\nabla} P + \hat{\chi}\boldsymbol{z}^\pm,\quad \hat{\chi} \doteq \lambda \partial_x\nabla^2, \end{equation}

where $\lambda$ is a constant parameter. The last term is intended as a simple ad hoc correction that, while causing deviation from ideal MHD, leaves the primary-mode evolution unaffected and preserves MHD's key invariants, specifically, the energy and cross-helicity (Appendix B). Perhaps notably, this term introduces rotational asymmetry in the $(x, y)$ plane. Similar terms can appear due to a background magnetic field or differential rotation (Heinonen et al. 2023).

Taking $\phi = {\rm \pi}/4$ for simplicity, one arrives at the following corresponding modification of the linear (2.23):

(6.2)\begin{equation} \sqrt{2}(\partial_\tau+ \mathrm{i} \delta_n)\psi_n^\pm = \mathrm{e}^{{\mp}\mathrm{i}\theta} \alpha_n^- \psi_{n-1}^\pm + \mathrm{e}^{{\pm} \mathrm{i}\theta}\beta_n^- \psi_{n-1}^\mp + \mathrm{e}^{{\pm}\mathrm{i}\theta}\alpha_n^+ \psi_{n+1}^\pm + \mathrm{e}^{{\mp} \mathrm{i}\theta}\beta_n^+ \psi_{n+1}^\mp, \end{equation}

where $\delta _n^\pm \doteq \varLambda r(r^2+n^2)$, $\varLambda \doteq \lambda /\mathcal {A}p^3$ and $\alpha _n^\pm$, $\beta _n^\pm$ are as in (2.24). In the large-$|n|$ limit, one has $\alpha _n^ \pm \sim 1$, $\beta _n^\pm \to 0$, $\delta _n \sim \varLambda n^2$, so one obtains the following scaling:

(6.3)\begin{equation} |\psi_{n+1}^\pm| \sim \frac{|\psi_n^\pm|}{\varLambda n^2}. \end{equation}

This shows that the harmonic magnitude $|\psi _n|$ decreases rapidly (superexponentially) with $|n|$, and thus low-order truncations may be justified. Note that although the opposite scaling, $|\psi _{n-1}^\pm | \sim |\psi _n^\pm |/\varLambda n^2$ is also formally possible, such modes cannot be excited, as is the case with inward propagating PSWs. Figure 13(a) shows that, indeed, the agreement between analytical QL growth rates and XQLT improves as the parameter $\varLambda$ increases. The agreement becomes nearly perfect for $\varLambda \gtrsim 0.5$. Figure 13(b) shows that the same effect can be achieved if one introduces sufficiently strong viscosity. In this case, one also has a modified (2.23) with the exact form of (6.2), but with $\hat {\chi }\doteq \nu _+\nabla ^2$ such that $\delta _n^\pm \doteq \mu (r^2+n^2)$, where $\mu \doteq \nu _+/\mathcal {A}p^2$ (with $\nu _-=0$ for simplicity). Again, the agreement becomes nearly exact for $\mu \gtrsim 0.5$.

Figure 13. The growth rate $\varGamma$, at $\nu _- = 0$, versus (a) $\varLambda \doteq \lambda /\mathcal {A}p^3$ and (b) $\mu \doteq \nu _+/\mathcal {A}p^2$. The corresponding ‘centre of energy’ of the unstable eigenmode, $\bar {n}\doteq \sqrt {\sum _n \mathcal {E}_n n^2/\mathcal {E}}$, is shown in (c,d), respectively. In all figures, the colour markers indicate the results obtained through XQLS, while the black solid curves indicate solutions obtained from the 4MT truncation of (6.2). The results are presented for the representative cases $\theta = 0$ and $\theta = {\rm \pi}/2$, both at $r = 0.5$ and $\phi = {\rm \pi}/4$.

More generally, note that any other linear operator $\hat {\chi }$ can be encapsulated as some $\delta _n$. The left-hand side of (6.2) therefore determines the linear response of the system, while the right-hand side corresponds to nonlinear (modulational) dynamics, and such a decomposition into linear and nonlinear can be done for any generic system. If the linear waves of the system are dispersive, i.e. $\hat {\chi }$ is a differential operator of the second or higher order, then $\delta _n$ grows with $|n|$. Then, unless the coefficients on the right-hand side of (6.2) also grow proportionally, the amplitudes $|\psi _n^\pm |$ decrease with $|n|$ exponentially or superexponentially and PSW are not possible.

7 Summary

In this paper, we explore structure formation in 2-D MHD turbulence as a MI of turbulent fluctuations. We focus on the early stages of structure formation and consider simple backgrounds that allow for a tractable model of the MI while retaining the full chain of modulational harmonics. This approach allows for a systematic examination of the importance of high-order correlations that are typically ignored in mean-field theories.

We find that, when the primary structure truly experiences a MI, this MI can be described well with typical truncated models that neglect high-order correlations. However, already in adjacent regimes, such truncated models can fail spectacularly and produce false positives for instability. To study this process in detail, we propose an XQLT that treats the primary structure as fixed but includes the entire spectrum of modulational harmonics (as opposed to just the low-order harmonics, as usual). From XQLT, also corroborated by DNS, we find that the difference between said regimes is due to a fundamental difference in the structures of the modulational spectra. For unstable modes, the spectrum is exponentially localized at low harmonic numbers, so truncated models are justified. But this localization does not always occur. At other parameters, modulational modes turn into constant-amplitude waves propagating down the spectrum, unimpeded until dissipative scales. These spectral waves are self-maintained as global modes with real frequencies and cause ballistic energy transport along the spectrum, breaking the feedback loops that could otherwise sustain MI.

The ballistic transport by PSWs drains energy from the primary structure at a constant rate until the primary structure is depleted. Because global PSWs exist at almost all parameters, this means that almost any primary structure in ideal incompressible MHD will eventually be depleted. This means, in particular, that sustainability of MHD structures is not entirely limited to the issue of exponentially growing linear instabilities, as PSWs must also be taken into consideration. To describe them within a reduced model, we propose a closure that successfully captures both the PSWs and MIs for the assumed primary structure.

Finally, we find that departures from ideal MHD constrains the form of modulational eigenmodes, in turn suppressing the amplitude and impact of high harmonics. This allows us to end on an informed yet optimistic note regarding the applicability of the quasilinear approximation for structure formation in dispersive forms of MHD. That said, it is an important conclusion of our work that, unless deviations from ideal incompressible MHD are substantial, changing the turbulence parameters even slightly can utterly destroy the applicability of an otherwise workable reduced model. An understanding of the complex modulational dynamics supported by (nearly) incompressible MHD provides important context in the interpretation of existing simple closures and potentially opens the path to building more reliable alternatives.