Formation of solitary zonal structures via the modulational instability of drift waves

TL;DR

This paper investigates how modulational instability of drift waves can lead to the formation of solitary or stationary zonal structures, using a nonlinear Schrödinger equation and phase-space modeling, highlighting quantum-like effects.

Contribution

It demonstrates the formation of solitary and stationary zonal structures via modulational instability in a modified Hasegawa-Mima framework, incorporating quantum-like phase-space effects.

Findings

Small amplitude drift waves produce solitary structures.

Large amplitude drift waves produce stationary structures.

Incoherent drift waves can also form zonal structures depending on beam intensity.

Abstract

The dynamics of the radial envelope of a weak coherent drift wave is approximately governed by a nonlinear Schr\"odinger equation, which emerges as a limit of the modified Hasegawa-Mima equation. The nonlinear Schr\"odinger equation has well-known soliton solutions, and its modulational instability can naturally generate solitary structures. In this paper, we demonstrate that this simple model can adequately describe the formation of solitary zonal structures in the modified Hasegawa-Mima equation, but only when the amplitude of the coherent drift wave is relatively small. At larger amplitudes, the modulational instability produces stationary zonal structures instead. Furthermore, we find that incoherent drift waves with beam-like spectra can also be modulationally unstable to the formation of solitary or stationary zonal structures, depending on the beam intensity. Notably, we show that these drift waves can be modeled as quantumlike particles ("driftons") within a recently developed phase-space (Wigner-Moyal) formulation, which intuitively depicts the solitary zonal structures as quasi-monochromatic drifton condensates. Quantumlike effects, such as diffraction, are essential to these condensates; hence, the latter cannot be described by wave-kinetic models that are based on the ray approximation.

Equations

$${\partial }_{t}w+\mathbf{v}\cdot \nabla w$$

$$w\doteq {\nabla }^2\phi$$

$${\partial }_t\tilde{w}+⟨\mathbf{v}⟩\cdot \nabla \tilde{w}+\tilde{\mathbf{v}}\cdot \nabla ⟨w⟩-\beta {\partial }_y\tilde{\phi }$$

$${\partial }_t⟨w⟩+⟨\tilde{\mathbf{v}}\cdot \nabla \tilde{w}⟩$$

$${\partial }_t\tilde{w}+U{\partial }_y\tilde{w}-(\beta +U^{\prime \prime }){\partial }_y\tilde{\phi }$$

$${\partial }_{t}U-{\partial }_x⟨{\partial }_x\tilde{\phi }{\partial }_y\tilde{\phi }⟩$$

$$i{\partial }_t\tilde{w}=\hat{H}\tilde{w}.$$

$$\hat{H}\doteq (\beta +{\hat{U}}^{\prime \prime }){\hat{p}}_y{\hat{\bar{p}}}^{-2}+\hat{U}{\hat{p}}_y,$$

$${\hat{H}}_0\approx \Omega +v_{\text{g}}{\hat{p}}_x+(\chi /2){\hat{p}}_x^2,{\hat{H}}_U\approx k_{y}U,$$

$$v_{\text{g}}=-2\beta k_x{k}_y/{\bar{k}}^4,\chi =(2\beta k_y/{\bar{k}}^6)(4k_x^2-{\bar{k}}^2).$$

$$i({\partial }_t+v_{\text{g}}{\partial }_x)\psi \approx -(\chi /2){\partial }_x^2\psi +k_{y}U\psi .$$

$${\partial }_{t}U\approx {\partial }_x⟨{\partial }_x\tilde{w}{\partial }_y\tilde{w}⟩/{\bar{k}}^4\approx k_x{k}_y{\partial }_x|\psi {|}^2/(2{\bar{k}}^4),$$

$$U\approx |\psi {|}^2/(4\beta ).$$

$$i({\partial }_t+v_{\text{g}}{\partial }_x)\psi \approx -(\chi /2){\partial }_x^2\psi +k_y|\psi {|}^2\psi /(4\beta ).$$

$$(\omega -q{v}_{\text{g}}{)}^2=\frac{{\chi }^2{q}^4}{4}\left(1+\frac{k_y|{\psi }_0{|}^2}{\beta \chi q^2}\right).$$

$${\partial }_{t}W$$

$${\partial }_{t}U$$

$$W\doteq ⟨\int {\mathrm{d}}^{2}s{e}^{-i\mathbf{p}\cdot \mathbf{s}}\tilde{w}\left(t,\mathbf{x}+\frac{\mathbf{s}}{2}\right)\tilde{w}\left(t,\mathbf{x}-\frac{\mathbf{s}}{2}\right)⟩.$$

$$\mathcal{H}$$

$$\Gamma$$

$$\omega =\sum _{\pm }\int \frac{{\mathrm{d}}^{2}p}{(2\pi {)}^2}\mathcal{W}(\mathbf{p})\frac{\mp q{p}_y^2(p_x\pm q/2)(1-q^2/{\bar{p}}^2)}{{\bar{p}}^2{\bar{p}}_{\pm }^2\omega +2\beta q{p}_y(p_x\pm q/2)},$$

$$\mathcal{W}(\mathbf{p})={\pi }^2|{\psi }_0{|}^2[\delta (\mathbf{p}+\mathbf{k})+\delta (\mathbf{p}-\mathbf{k})].$$

$$\prod _{\pm }$$

$$|{\psi }_0{|}^2{q}^2{k}_y^2({\bar{k}}^2-q^2)(4k_x^2-{\bar{k}}^2-q^2)/2.$$

$$\psi =2\eta \sqrt{-\frac{\beta \chi }{k_y}}\frac{\exp (i\chi {\eta }^{2}t/2)}{\cosh [\eta (x-v_{\text{g}}t)]}.$$

$$\tilde{w}$$

$$U$$

$$\mathcal{W}(\mathbf{p})={\pi }^2|{\psi }_0{|}^2\sum _{\pm }\delta (p_y\pm k_y)f_{\sigma }(p_x\pm k_x),$$

$$\mathcal{W}(\mathbf{p})={\pi }^2|{\psi }_0{|}^2\sum _{\pm }\delta (p_y\pm k_y)[f_{\sigma }(p_x\pm k_x)+f_{\sigma }(p_x\mp k_x)].$$

$$\prod _{\pm }$$

$$=$$

$$\times [{\bar{k}}^4{\bar{k}}_{+}^2{\bar{k}}_{-}^2{\omega }^2+{\beta }^2{q}^2{k}_y^2(4k_x^2-q^2)].$$

$${\partial }_{t}W$$

$${\partial }_{t}U$$

$$\mathcal{H}$$

Full text

I Introduction

In magnetically confined plasmas, radially propagating coherent structures are of great interest, as they can induce transport that is ballistic rather than diffusive. Examples include turbulence spreading Garbet et al. (1994); Lin and Hahm (2004); Hahm et al. (2004); Naulin et al. (2005); Guo et al. (2009) and avalanching Politzer (2000); Candy and Waltz (2003); McMillan et al. (2009), as well as the density “blobs” in edge plasmas Zweben et al. (2007); Krasheninnikov et al. (2008); Krasheninnikov (2016); Zhang and Krasheninnikov (2017). Recently, such structures have also been identified in gyrokinetic simulations of subcritical plasmas van Wyk et al. (2016, 2017); McMillan et al. (2018).

An arguably simplest model of radially propagating coherent structures considers the radial envelope dynamics of a weak coherent drift wave (DW) without forcing and dissipation. It has been shown that the envelope approximately follows a nonlinear Schrödinger equation (NLSE) Champeaux and Diamond (2001); Dewar and Abdullatif (2007), where the (cubic) nonlinearity originates from the quasilinear interaction between the primary DW and a secondary zonal flow (ZF). The well-known soliton solution to the NLSE corresponds to a zonal structure that propagates radially at the DW group velocity. In particular, Guo et al. Guo et al. (2009) first studied this DW–ZF soliton in the context of turbulence spreading from linearly unstable regions to stable regions. Accordingly, they examined soliton formation due to the inhomogeneity of the linear growth rate. However, an intrinsic mechanism for the formation of solitary zonal structures is still needed, since events such as turbulence avalanching seem to be dominated by local physics McMillan et al. (2009). One natural candidate is the modulational instability (MI) Champeaux and Diamond (2001); Dewar and Abdullatif (2007), which is known to generate solitary structures in the NLSE. [Here, “solitary” means propagating at a (roughly) constant speed while maintaining a (roughly) constant shape.] Nevertheless, the relevance of this mechanism in the modified Hasegawa–Mima equation (mHME) Dewar and Abdullatif (2007); Chandre et al. (2014), the parent model of the NLSE, has remained unexplored.

In this paper, we explicitly demonstrate that the NLSE can adequately describe the formation of solitary zonal structures via the MI in the mHME, using both quasilinear and nonlinear simulations of the latter. However, these structures only emerge from primary DWs with relatively small amplitudes; at larger amplitudes, the MI produces stationary zonal structures instead. Then, as a generalization, we simulate the MI of incoherent DWs with beam-like spectra, using a recently developed Wigner–Moyal (WM) model of DW–ZF dynamics Ruiz et al. (2016). While the finite beam width has a stabilizing effect on the MI, the results are similar to those from coherent DWs. That is, with moderate beam intensity, solitary zonal structures are formed; as the intensity increases, the zonal structures become stationary.

One novelty of the WM model is that it treats DWs as quantumlike particles (“driftons”) and facilitates analysis of the zonal structures from an instructive phase-space perspective. The Wigner function of the DW–ZF solitons show concentration of DW quanta in both position and momentum, depicting them as self-trapping drifton condensates, akin to Bose–Einstein condensates in quantum mechanics. The MI-induced solitary zonal structures exhibit similar features, which suggests that they are essentially the DW–ZF solitons. In turn, it also implies that these structures are not the same as those obtained from wave-kinetic models Smolyakov and Diamond (1999); Trines et al. (2005), which are the conventional tools for treating incoherent DWs. The reason is that diffraction, among other quantumlike effects, is fundamental to these condensates (in fact, it determines their sizes), but neglected by the ray approximation that the wave-kinetic models invoke Diamond et al. (2005).

This paper is organized as follows. In Sec. II, we introduce the mHME and its quasilinear approximation. In Sec. III, we consider coherent DWs by deriving the NLSE and its dispersion relation of the MI. In Sec. IV, we describe the WM model and the MI of general DW spectra. Most of our new results are presented in Sec. V, where we study solitary zonal structures and their formation from coherent DWs and incoherent DW spectra. Our results are summarized and discussed in Sec. VI.

II Basic model

II.1 Modified Hasegawa–Mima equation

In this paper, we study DWs within the modified Hasegawa–Mima equation (mHME) Dewar and Abdullatif (2007); Chandre et al. (2014), which is the simplest yet useful model that captures many basic effects of interest. (The mHME has been referred to as the generalized Krommes and Kim (2000) or extended Connaughton et al. (2015) Hasegawa–Mima equation as well.) In a dimensionless form, the mHME can be written as

[TABLE]

It is a 2D model in slab geometry, with coordinates $\mathbf{x}\equiv(x,y)$ and a uniform magnetic field normal to the plane. The gradient of the plasma density $n_{0}$ is in the radial ($x$) direction, and is parameterized by a (positive) constant $\beta\doteq a/L_{n}$, where $a$ is the system scale length and $L_{n}\doteq(-\mathrm{d}\ln n_{0}/\mathrm{d}x)^{-1}$ is the local scale length of the gradient. (The symbol $\doteq$ denotes definitions.) The ZF is in the poloidal ($y$) direction. Time $t$ is normalized by the transit time $a/c_{\text{s}}$, where $c_{\text{s}}$ is the sound speed. Space is normalized by the ion sound radius $\rho_{\text{s}}\doteq c_{\text{s}}/\Omega_{\text{ci}}$, where $\Omega_{\text{ci}}$ is the ion gyro-frequency. The electrostatic potential $\varphi(t,\mathbf{x})$ is normalized by $T_{\text{e}}\rho_{\text{s}}/(ea)$, where $e$ is the unit charge and $T_{\text{e}}$ is the electron temperature. Accordingly, $\mathbf{v}\doteq(-\partial_{y}\varphi,\partial_{x}\varphi)$ is the $\mathbf{E}\times\mathbf{B}$ velocity.

In the mHME, the definition of the generalized vorticity $w$ (1b) involves separating the total $\varphi$ into the zonal component $\langle{\varphi}\rangle$ and non-zonal component $\tilde{\varphi}$. The former is the “zonal average” of $\varphi$, $\langle{\varphi}\rangle\doteq\int\mathrm{d}y\,\varphi/L_{y}$ (where $L_{y}$ is the domain length in $y$), and corresponds to the ZF. The latter is the fluctuating component, $\tilde{\varphi}\doteq\varphi-\langle{\varphi}\rangle$, and corresponds to DWs. The same notations apply to $w$ and $\mathbf{v}$ as well.

In contrast, in the original Hasegawa–Mima equation (oHME) Hasegawa and Mima (1978), the generalized vorticity is defined as $w\doteq\nabla^{2}\varphi-{\varphi}$. The modification in the mHME is due to the finding that the zonal potential $\langle{\varphi}\rangle$ does not contribute to the adiabatic electron response Dorland and Hammett (1993); Hammett et al. (1993). The oHME is also called the Charney–Hasegawa–Mima equation for its equivalence to the Charney equation Charney (1971). Meanwhile, with $w\doteq\nabla^{2}\varphi$, Eq. (1a) becomes equivalent to the barotropic vorticity equation Charney et al. (1950). Both the Charney equation and the barotropic vorticity equation are widely used in studies of geophysical fluids. The similarity between the mHME and these equations suggests that our results can, to an extent, be relevant to Rossby-wave turbulence in geophysics.

The mHME does not have a primary instability, i.e., an instability that generates DWs. Thus, external forcing is sometimes introduced as a proxy, and ad hoc dissipation must also be added to balance the energy input. However, due to the existence of the drift term (the linear term proportional to $\beta$), the mHME can support finite-amplitude DWs even in the absence of forcing. Then, ZFs can emerge from these DWs through a secondary instability, which is also known as the MI Champeaux and Diamond (2001); Dewar and Abdullatif (2007) or the zonostrophic instability Srinivasan and Young (2012). In this paper, we focus on this particular process and will not be concerned with the origin of the DWs. Instead, we will introduce finite-amplitude DWs via initial conditions, and exclude forcing and dissipation (except briefly in the end of Sec. V.3).

II.2 Quasilinear approximation

To proceed, let us separate Eq. (1a) into the non-zonal and zonal components, respectively:

[TABLE]

In the studies of ZF formation, it is common to assume that the nonlinearity on the right-hand side (RHS) of Eq. (2a) is weak, which can be physically interpreted as neglecting the eddy–eddy interactions between the DWs. This is called the quasilinear approximation, for it leads to two linear equations that are nonlinearly coupled:

[TABLE]

For convenience, we introduce the ZF velocity $U(t,x)\doteq\langle v_{y}\rangle=\partial_{x}\langle\varphi\rangle$, with $U^{\prime\prime}\doteq\partial^{2}_{x}U$.

The quasilinear mHME (3) has been shown to reproduce many of the basic features of the original nonlinear system (2), at least qualitatively Srinivasan and Young (2012). Hence, we consider the quasilinear approximation sufficient for our purposes, and adopt it throughout the rest of the paper (except in Figs. 4 and 5, where we briefly present some nonlinear simulation results as verifications for the corresponding quasilinear simulations).

The non-zonal component of the quasilinear mHME (3a) can also be written in the form of a Schrödinger equation for DW quanta (driftons),

[TABLE]

Unlike the truly quantum Schrödinger equation, Eq. (4) does not contain $\hbar$, so it is purely classical. (Likewise, $|\tilde{w}|^{2}$ is an action density rather than number density.) Also, the Hamiltonian operator,

[TABLE]

is not entirely Hermitian, since $\hat{U}^{\prime\prime}$ does not commute with $\hat{\bar{p}}^{-2}$. Here, $\hat{\mathbf{p}}\doteq-i\nabla$ can be understood as the momentum (wave-vector) operator, and $\hat{\bar{p}}^{2}\doteq 1+\hat{p}_{x}^{2}+\hat{p}_{y}^{2}$ such that $\tilde{w}=-\hat{\bar{p}}^{2}\tilde{\varphi}$. Also, $\hat{U}\doteq U(t,\hat{x})$ and $\hat{U}^{\prime\prime}\doteq U^{\prime\prime}(t,\hat{x})$, with $\hat{\mathbf{x}}$ being the position operator. This quantumlike formalism proves useful in deriving the nonlinear Schrödinger equation that governs the envelope dynamics of coherent DWs (Sec. III), as well as the WM formulation that can describe incoherent DWs (Sec. IV).

III Coherent drift waves

The quasilinear mHME (3) has an exact plane-wave solution with finite amplitude, $\tilde{w}=\text{Re}(\psi_{0}e^{i\mathbf{k}\cdot\mathbf{x}-i\Omega t})$, where $\mathbf{k}\doteq(k_{x},k_{y})$ is the wave-vector and $\Omega\doteq\beta k_{y}/\bar{k}^{2}$ is the DW frequency, with $\bar{k}^{2}\doteq 1+k_{x}^{2}+k_{y}^{2}$ and $\psi_{0}$ being a complex constant denoting the amplitude.

This primary wave, when subject to large-scale modulations, can become unstable. One simplified way to study this MI is to consider the envelope dynamics of a coherent DW. The first study of such kind appears to be Ref. Nozaki et al. (1979), which is based on the oHME and only considers poloidal modulations (“streamers”). In Ref. Champeaux and Diamond (2001), the mHME is employed, and radial and poloidal modulations are treated on the same footing. More comprehensive reviews of the envelope formalism can be found in Refs. Diamond et al. (2005); Dewar and Abdullatif (2007). Also notably, related equations were later rediscovered independently in Refs. Guo et al. (2009) (with over-simplified coefficients) and Jovanovic et al. (2010).

In all of these studies, it is noticed that the envelope equation can be approximated as a NLSE, which is well-known to have a MI. In Sec. III.1, we show how the NLSE follows naturally from our quantumlike formalism. In Sec. III.2, we rederive the corresponding dispersion relation of the MI.

III.1 Nonlinear Schrödinger equation

Let us represent the Hamiltonian operator (5) as $\hat{H}=\hat{H}_{0}+\hat{H}_{U}$, where $\hat{H}_{0}$ is the $U$-independent part and $\hat{H}_{U}$ scales linearly with $U$. Since we focus on zonal structures, here we consider a coherent DW with radial modulation only, $\tilde{w}=\text{Re}[\psi(t,x)e^{i\mathbf{k}\cdot\mathbf{x}-i\Omega t}]$. We assume that the envelope $\psi$ is slow, i.e., $|\partial_{x}\ln\psi|\ll|k_{x}|$, and also that $|U|$ is small. Then, the Hamiltonian operator (5) can be approximated as

[TABLE]

where $v_{\text{g}}\doteq\partial\Omega/\partial k_{x}$ is the radial group velocity and $\chi\doteq\partial^{2}\Omega/\partial k_{x}^{2}$, or explicitly,

[TABLE]

Equation (6) can be viewed as a special case of the Weyl expansion derived in Ref. Dodin et al. (2019) for an inhomogeneous medium. Additionally, the term proportional to $U^{\prime\prime}$ has been neglected because both $|U|$ and $\partial^{2}_{x}$ are assumed small. Then, the resulting equation for $\psi$ is

[TABLE]

Meanwhile, using $\tilde{\varphi}\approx-\tilde{w}/\bar{k}^{2}$, we can approximate the ZF equation (3b) as

[TABLE]

where the factor $1/2$ originates from zonal averaging. From Eq. (8), we can see that, to the leading order, the modulation propagates at the group velocity $v_{\text{g}}$. Hence, we can assume that $\partial_{t}\approx-v_{\text{g}}\partial_{x}$ in Eq. (9), and integrate in $x$ to obtain $U\approx-k_{x}k_{y}|\psi|^{2}/(2\bar{k}^{4}v_{\text{g}})$, or more explicitly,

[TABLE]

Here, vanishing boundary condition in $x$ is implied. For other boundary conditions (e.g., $\langle\varphi\rangle$ periodic in $x$), an additional integration constant may be needed on the RHS of Eq. (10) for consistency. (This constant can be easily removed by a Galilean transformation, however.) Also, with $|U|$ assumed small, it is implied that the DW amplitude $|\psi|$ should be small too. Substituting Eq. (10) into Eq. (8), we obtain

[TABLE]

Equation (11) has previously been derived (using somewhat different approaches) in Refs. Champeaux and Diamond (2001); Dewar and Abdullatif (2007). It has the form of a NLSE, or the Gross–Pitaevskii equation, so the structures it describes can be viewed as “drifton condensates” (by analogy with the Bose–Einstein condensate). Namely, Eq. (11) shows that it is energetically favorable for driftons to be in a correlated state rather than have random phases. Also, Eq. (10) can be interpreted as an equation of state of the condensates, as it provides a local relation between the drifton density $|\psi|^{2}/2$ and another “thermodynamic” property of the condensates, $U$.

In Fig. 1, we compare the evolution of an initially Gaussian envelope in numerical simulations of the NLSE and the quasilinear mHME. The good agreement between the solutions confirms that the former is a reasonable approximation of the latter. All of our simulations using configuration-space models (the NLSE, the quasilinear mHME, and the nonlinear mHME) are pseudo-spectral, dealiased, and performed on periodic domains.

III.2 Modulational instability

The NLSE (11) has an exact homogeneous solution $\psi=\psi_{0}\exp[-ik_{y}|\psi_{0}|^{2}t/(4\beta)]$. The frequency $\omega$ and wavenumber $q$ of a linear perturbation on this solution satisfies the following dispersion relation Dewar and Abdullatif (2007),

[TABLE]

When $\beta\chi/k_{y}\propto 4k_{x}^{2}-\bar{k}^{2}<0$, the frequency is complex for small $q$ and the solution is linearly unstable, with the wavenumber of the fastest-growing mode given by $q_{\text{max}}=|\psi_{0}|\sqrt{-k_{y}/(2\beta\chi)}$. This is the well-known MI of the NLSE, arising from the interplay of diffraction and self-focusing [the first and second terms on the RHS of Eq. (11), respectively].

The NLSE (11) offers an intuitive perspective on the MI of coherent DWs. However, as an approximate model, it is restricted to slow modulations and small $|\psi_{0}|$. When $|\psi_{0}|$ is large, $q_{\text{max}}$ can be comparable or larger than $k_{x}$, which is inconsistent with the underlying assumption of the NLSE. In addition, the NLSE only applies to primary waves with non-zero radial group velocity $v_{\text{g}}\propto k_{x}k_{y}$. While in this paper we focus on such waves for this very feature, primary waves with $k_{x}=0$ are also of interest, since they correspond to the fastest growing modes in some primary instabilities, particularly, ion-tempatature-gradient modes Chen et al. (2000); Rogers et al. (2000).

In fact, there are more general approaches to deriving the dispersion relation of the MI. One way is to employ the four-mode truncation (4MT) method. As the name suggests, the 4MT is a truncation of the mHME in spectral representation by only keeping four modes: a primary wave with wave-vector $\mathbf{k}$, a modulation with wave-vector $\mathbf{q}$, and two sidebands with wave-vectors $\mathbf{k}\pm\mathbf{q}$. In general, the modulation does not have to be purely radial. For example, the MI with a purely radial primary wave (a ZF) and a purely poloidal modulation is a tertiary instability of the ZF Zhu et al. (2018a, b). Detailed discussions on the 4MT can be found in Refs. Connaughton et al. (2010); Gallagher et al. (2012); Zhu et al. (2019). Meanwhile, for purely radial modulations with $\mathbf{q}=(q,0)$, which we focus on in this paper, the MI of coherent DWs can be considered as a special case of the MI of general DW spectra, which is discussed in Sec. IV.

IV Drift-wave ensembles

Equation (4), along with Eq. (3b), governs the quasilinear dynamics of a single realization of DWs (in quantum mechanical terms, a “pure state”). However, due to the incoherent nature of DW turbulence, it is useful to consider the dynamics of an ensemble of DWs statistically. This is equivalent to studying the von Neumann equation that follows from the Schrödinger equation (4), which can describe the dynamics of “mixed states”. In double-configuration-space representation, this leads to the theory of second-order cumulant expansion (CE2), which has been widely used in geophysical fluids (e.g., Farrell and Ioannou (2003); Marston et al. (2008); Srinivasan and Young (2012); Parker and Krommes (2014)) and subsequently introduced to plasma physics Farrell and Ioannou (2009); Parker and Krommes (2013). A mathematically equivalent yet physically more intuitive alternative to the CE2 model is the phase-space representation of the von Neumann equation. This leads to the Wigner–Moyal (WM) model of DW–ZF dynamics, which we introduce in Sec. IV.1. The WM model can describe the MI of general DW spectra, which is presented in Sec. IV.2.

IV.1 Wigner–Moyal formulation

The WM model of DW–ZF dynamics was first derived in Ref. Ruiz et al. (2016). The derivation starts from the quasilinear mHME (3) and leads to the following equations:

[TABLE]

Here, $\mathbf{p}\doteq(p_{x},p_{y})$ is the coordinate in the DW momentum (wave-vector) space, and $W(t,x,\mathbf{p})$ is the zonal-averaged Wigner function Wigner (1932). For a single realization of $\tilde{w}(t,\mathbf{x})$, $W$ can be written as

[TABLE]

For an ensemble of realizations, the zonal average (again, denoted by the angle bracket) can be regarded as an ensemble average. The Wigner function $W$ is a quasi-probability distribution of driftons, and the ZF velocity $U$ serves as a collective field through which they interact. Since $\tilde{w}$ is real (unlike in quantum mechanics where the wave functions are complex), the DW Wigner function has a unique parity in $\mathbf{p}$ that $W(t,x,\mathbf{p})=W(t,x,\mathbf{-p})$. This implies that driftons come in pairs, i.e., each drifton with wave-vector $\mathbf{p}$ has a twin with wave-vector $-\mathbf{p}$.

The specific dynamics of the driftons is governed by

[TABLE]

which are the Hermitian and anti-Hermitian parts of the Hamiltonian, respectively. Here, the Moyal star product $A(\mathbf{x},\mathbf{p})\star B(\mathbf{x},\mathbf{p})\doteq A\exp(i\hat{\mathcal{L}}/2)B$ Moyal (1949), where the Janus operator $\hat{\mathcal{L}}$ is defined as $A\hat{\mathcal{L}}B\doteq(\partial_{\mathbf{x}}A)\cdot(\partial_{\mathbf{p}}B)-(\partial_{\mathbf{p}}A)\cdot(\partial_{\mathbf{x}}B)$. The Moyal sine and cosine brackets are given by $\{\!\{A,B\}\!\}\doteq 2A\sin(\hat{\mathcal{L}}/2)B$ and $[[A,B]]\doteq 2A\cos(\hat{\mathcal{L}}/2)B$, respectively. A detailed derivation of Eq. (13) and a review of the Weyl calculus can be found in Ref. Ruiz et al. (2016).

The WM equation (13) can be understood as a kinetic model treating driftons as quantumlike particles with finite “de Broglie” wavelengths. As such, it captures “full-wave” effects missing in wave-kinetic models of DW–ZF dynamics based on the ray approximation Smolyakov and Diamond (1999); Parker (2016), which treat driftons as point particles with zero wavelengths. While the wave-kinetic models prove useful in some scenarios Parker (2018); Zhu et al. (2018c); Ruiz et al. (2019), they are insufficient for the problems that we study in this paper. A detailed discussion on the limitations of the wave-kinetic models is presented in Appendix A.

IV.2 Modulational instability

In the WM model, a statistically homogeneous equilibrium can be described by a Wigner function $W=\mathcal{W}(\mathbf{p})$, which can be interpreted as a DW spectrum. Linearizing Eq. (13) about $\mathcal{W}(\mathbf{p})$ leads to the following dispersion relation of the MI Ruiz et al. (2016); Zhu et al. (2018a),

[TABLE]

where $\bar{p}_{\pm}^{2}\doteq 1+(p_{x}\pm q)^{2}+p_{y}^{2}$. The MI of a monochromatic DW $\tilde{w}=\text{Re}(\psi_{0}e^{i\mathbf{k}\cdot\mathbf{x}-i\Omega t})$ can be considered as a special case. Using Eq. (14), we can obtain the corresponding spectrum of this DW

[TABLE]

Substituting Eq. (17) into Eq. (16) leads to the dispersion relation of the MI of this monochromatic DW,

[TABLE]

Note that this dispersion relation agrees with that given by the 4MT method Connaughton et al. (2010); Gallagher et al. (2012), and is the exact dispersion relation of the quasilinear mHME (3).

It is instructive to compare the exact dispersion relation (18) with the approximate one (12) obtained from the NLSE. In Fig. 2, we show both dispersion relations (for a given $\mathbf{k}$) with multiple values of $|\psi_{0}|$. It can be seen that the agreement between the dispersion relations is better for small $|\psi_{0}|$ and small $q$, which is consistent with the fact that the NLSE is derived based on the assumptions of slow modulation and small DW amplitude. Still, in Fig. 2, the NLSE appears to be a reasonable approximation of the quasilinear mHME even for relatively large $|\psi_{0}|$, in terms of capturing the linear MI. (This is not the case in the nonlinear stage of the MI, as we will discuss in Sec. V.) For comparison, Fig. 2 also shows the dispersion relations obtained from a wave-kinetic model based on the ray approximation (overviewed in Appendix A), which do not accurately approximate the exact ones. The reason for the discrepancy is that the wave-kinetic model misses essential full-wave (quantumlike) effects, particularly diffraction.

From Fig. 2, we can see that the unstable modulations have finite real frequencies due to the DW group velocity, and hence propagate while growing. This feature has largely been overlooked in the past, possibly because zonal structures are typically perceived as (quasi-) stationary. Nevertheless, as we will show in Sec. V, the propagation of these zonal modes can have consequences in the nonlinear stage of the MI, leading to the formation of solitary zonal structures.

V Solitary zonal structures

V.1 Drift-wave–zonal-flow soliton

Since Eq. (11) is of the same form as the well-known NLSE, it has the usual soliton solution

[TABLE]

Here, $\eta$ is a free parameter that can be regarded as the soliton inverse width. With $\tilde{w}=\text{Re}[\psi(t,x)e^{i\mathbf{k}\cdot\mathbf{x}-i\Omega t}]$ and the equation of state (10), Eq. (19) translates to an approximate soliton solution to the quasilinear mHME (3):

[TABLE]

To our knowledge, this DW–ZF soliton was first explicitly discussed in Ref. Guo et al. (2009) and then in Ref. Jovanovic et al. (2010), even though the NLSE (11) that governs the DW envelope dynamics had been derived earlier Champeaux and Diamond (2001); Dewar and Abdullatif (2007). It is fundamentally different from the vortex-pair solution called “modon” Meiss (1983); Horton and Hasegawa (1994), which is a 2D structure that propagates poloidally, whereas the DW–ZF soliton is an essentially 1D structure that propagates radially.

In Fig. 3, we show snapshots of a DW–ZF soliton in both configuration space (1-a) and phase space (1-b). The Wigner function $W$, obtained using Eq. (14), is concentrated in both position and momentum, illustrating that the soliton is a quasi-monochromatic drifton condensate. It is also shown that the contours of $W$ do not coincide with those of the Hamiltonian in the moving frame $\mathcal{H}_{\text{m}}\doteq\mathcal{H}-v_{\text{g}}p_{x}$. This distinguishes the DW–ZF soliton from the BGK-type structures obtained from wave-kinetic models based on the ray approximation (Appendix A), where $W$ is a function of $\mathcal{H}_{\text{m}}$ only Smolyakov et al. (2000); Kaw et al. (2002); Itoh et al. (2005); Singh et al. (2014). After all, the DW–ZF soliton is a result of the balance between self-focusing and diffraction, and the latter is a quantumlike effect missing in wave-kinetic models Diamond et al. (2005). Therefore, wave-kinetic models cannot describe the DW–ZF soliton (20), even though the requirement that its envelope be slow may seem consistent with the ray approximation.

Also shown in Fig. 3 (row 2) are snapshots of two superposing solitons with wave-vectors $(k_{x},k_{y})$ and $(-k_{x},k_{y})$, respectively. Accordingly, these solitons have opposing group velocities and hence counter-propagate. Upon colliding, they tunnel through each other. The striations between the drifton condensates in (2-b) are signatures of the quantum superposition of macroscopically distinct states, i.e., “cat states” Kenfack and Yczkowski (2004); Weinbub and Ferry (2018).

In principle, $\eta\ll|k_{x}|$ is required for the soliton solution (20) to stand in the quasilinear mHME, since the NLSE is derived under the assumptions of slow envelope modulation and small DW amplitude. In practice, however, we find that the solitary behavior of this solution is quite robust, and extends even to $\eta\sim|k_{x}|$. In Fig. 4 (row 1), we show the spatial-temporal evolution of the DW envelope from quasilinear mHME simulations initialized with Eq. (20) for various values of $\eta$. As $\eta$ is doubled from (1-a) to (1-b), the zonal structure keeps propagating much like a soliton, while some small-amplitude structures emerge and the speed decreases slightly. As $\eta$ is tripled and quadrupled in (1-c) and (1-d), respectively, the zonal structure gradually breaks down and eventually stops propagating.

In Fig. 4 (row 2), we also show corresponding results from nonlinear mHME simulations to verify the quasilinear simulations. The same qualitative features can be observed: the solitary behavior is robust at relatively small $\eta$, while the propagation eventually stops as $\eta$ keeps increasing. Admittedly, at larger $\eta$, the propagation stops more quickly in the nonlinear simulations. Nonetheless, for solitary structures with smaller $\eta$, which we focus on in this paper, the quasilinear approximation proves reasonable and acceptable.

Having confirmed that the mHME indeed allows for solitary zonal structures, next, we demonstrate how they can spontaneously form via the MI, from either coherent DWs (Sec. V.2) or incoherent DW spectra (Sec V.3). These are the main results of this paper.

V.2 Zonal structures from coherent drift waves

In Fig. 5, we present numerical simulations of the MI, through the nonlinear stage, using three different models: the NLSE, the quasilinear mHME, and the nonlinear mHME. The initial states are chosen to be primary waves with various amplitudes ($|\psi_{0}|$) and random perturbations on top to seed the instability. The amplitudes of the perturbations are the same in all simulations.

In the NLSE simulations (row 1), for all values of $|\psi_{0}|$, the MI leads to the formation of zonal structures that behave approximately as solitons, notwithstanding the increasingly apparent nonlinear oscillations and interactions with increasing $|\psi_{0}|$. The wavenumbers of the structures are consistent with the fastest-growing wavenumbers in Fig. 2. The amplitudes of the structures increase with $|\psi_{0}|$ while the widths decrease, which is also consistent with the properties of the soliton solution (19). Note that in (1-c) and (1-d), the amplitudes of the structures are already comparable to those that stop propagating in Figs. 4(c) and 4(d), respectively. The implication is that in the corresponding mHME simulations, we may not be able to observe these solitary zonal structures.

Indeed, this can be seen in both the quasilinear (row 2) and nonlinear (row 3) simulations in Fig. 5. At small $|\psi_{0}|$ [column (a)], the spatial-temporal evolutions of the DW envelope closely resemble their NLSE counterpart (1-a) and shows solitary zonal structures. As $|\psi_{0}|$ increases [column (b)], the nonlinear interactions between the zonal structures begin to disrupt their propagation. At even larger $|\psi_{0}|$ [columns (c) and (d)] the zonal structures cease to propagate quickly after they develop, in contrast to the NLSE cases [(1-c) and (1-d)]. Again, we notice that the quasilinear and nonlinear simulation results are qualitatively similar, demonstrating the sufficiency of the quasilinear approximation for our purposes.

In summary, in the nonlinear stage of the MI of coherent DWs, the NLSE does not properly approximate the (quasilinear) mHME when the primary-wave amplitude is relatively large. Only when the primary-wave amplitude is relatively small can propagating zonal structures similar to the DW–ZF solitons form.

V.3 Zonal structures from incoherent drift waves

As discussed in Sec. IV.2, the MI of a monochromatic DW, which is studied in Sec. V.2, is equivalent to that of the delta-shaped DW spectrum (17) in the WM model. Now, let us consider a slightly more general case, where the DW spectrum has a finite width in $p_{x}$. Specifically, we adopt

[TABLE]

with $f_{\sigma}(p)\doteq{\exp[{-{p}^{2}/(2\sigma^{2})}}]/\sqrt{2\pi\sigma}$. Here, we keep the distribution in $p_{y}$ as delta functions for simplicity. The justification is that within the quasilinear approximation, $p_{y}$ is a constant of motion, and the coupling in $p_{y}$ is weak. For convenience, we refer to the spectrum (21) as “2-beam” for it consists of a pair of quasi-monochromatic drifton beams. (As mentioned in Sec. IV.1, driftons come in pairs.) Equation (17) is reproduced as the limit of Eq. (21) as the beam width $\sigma\rightarrow 0$.

In Fig. 6 (row 1), we present numerical simulations of the MI of the 2-beam spectrum (21) with random perturbations. These simulations employ the spectral representation of the WM equation (13) derived in the appendix of Ref. Ruiz et al. (2016). The simulation domain is periodic in $x$. In contrast to Fig. 5(2-a), there is no instability at small $|\psi_{0}|$ (1-a), which demonstrates the stabilizing effect of the finite beam width $\sigma$. Still, as the effective amplitude $|\psi_{0}|$ increases, the system becomes modulationally unstable, and the corresponding features are qualitatively similar to the quasilinear mHME simulations in Fig. 5 (row 2). With moderate $|\psi_{0}|$ (1-b), solitary zonal structures emerge. When $|\psi_{0}|$ is increased further [(1-c) and (1-d)], the zonal structures stop propagating and eventually become stationary.

As the next generalization, we consider the MI of a DW spectrum that has two pairs of quasi-monochromatic drifton beams,

[TABLE]

For convenience, we refer to the spectrum (22) as “4-beam”. In the limit as $\sigma\rightarrow 0$, this spectrum corresponds to the mixed state of two plane waves with amplitude $|\psi_{0}|$ and wave-vectors $(k_{x},k_{y})$ and $(-k_{x},k_{y})$, respectively. In this limit, the dispersion relation of the MI reads Ruiz et al. (2016)

[TABLE]

In Fig. 6 (row 2), we show numerical simulations of the MI of the 4-beam spectrum (22) with random perturbations. On one hand, we observe many features common with the 2-beam case (row 1): the system is stable at small $|\psi_{0}|$ (2-a); solitary zonal structures form as $|\psi_{0}|$ increases (2-b); when $|\psi_{0}|$ is even larger, the zonal structures become stationary [(2-c) and (2-d)]. On the other hand, unlike in the 2-beam case, counter-propagating solitary zonal structures emerge. This is due to the fact that the group velocities of the two drifton beams have opposite signs. The implication of the 4-beam case is that solitary zonal structures can also form via the MI of multiple (pairs of) quasi-monochromatic drifton beams.

Furthermore, we model the 4-beam case with some external forcing $F$ added to the RHS of Eq. (13a). To balance the energy input, we also add frictional damping $-2\mu W$ and $-\mu U$ to the RHS of Eqs. (13a) and (13b), respectively Ruiz et al. (2016). In order to compare directly with the unforced 4-beam case above, we choose $F=2\mu\mathcal{W}(\mathbf{p})$, with $\mathcal{W}$ given by the 4-beam spectrum (22), such that the corresponding homogeneous equilibrium is $W=F/(2\mu)=\mathcal{W}$. The simulations are initialized with such equilibria in place and random perturbations on top, and the results are shown in Fig. 6 (row 3). Many qualitative features of the unforced case (row 2) are reproduced here: stability at small $|\psi_{0}|$ (3-a), solitary zonal structures at moderate $|\psi_{0}|$ (3-b), and stationary zonal structures at large $|\psi_{0}|$ (3-d). Hence, we conclude that the formation of solitary zonal structures is still possible even with forcing and dissipation, provided that the forcing spectrum consists of quasi-monochromatic peaks, producing quasi-monochromatic distributions of driftons. That being said, the forced case visibly differs from the unforced case in that the former has larger amplitudes of DWs between the zonal structures, which is maintained by the homogeneous production of driftons by the external forcing. Another distinction is that the solitary zonal structures in the forced case (3-b) are more coherent than those in the unforced case (2-b), which also owes to the fact that the external forcing tends to keep the DWs quasi-monochromatic.

V.4 Phase-space structures

It is instructive to examine these zonal structures in phase space by studying snapshots of the Wigner function. Figure 3(1-b) illustrates that a DW–ZF soliton is a quasi-monochromatic drifton condensate that is concentrated in both space and momentum (at some nonzero $p_{x}$). For comparison, the phase-space snapshots in Fig. 7 (row 1) are taken from solitary zonal structures in various simulations. They all reveal the presence of trains of such condensates located at $p_{x}\approx(\pm)k_{x}$. In the quasilinear mHME simulation (1-a), the striations between the drifton condensates are signatures of “cat states”, akin to those in Fig. 3(2-b). In the 4-beam cases [(1-c) and (1-d)], the phase-space snapshots show two trains of drifton condensates with opposing group velocities, also similar to Fig. 3(2-b), such that the zonal structures counter-propagate.

In contrast, the phase-space snapshots in Fig. 7 (row 2) are taken from stationary zonal structures in different simulations. Accordingly, in these cases, the DW quanta are mostly localized at $p_{x}=0$, which is consistent with the stationarity of the zonal structures. Finally, we note that the phase-space structures are more coherent in the forced cases [column (d)] than in the unforced cases [column (c)], due to the homogeneous quasi-monochromatic external forcing applied.

VI Summary and discussion

In this paper, we first consider coherent DWs, by rederiving the NLSE that approximately governs their radial envelope dynamics, and subsequently, its dispersion relation of the MI and soliton solution. Using mHME simulations, both quasilinear and nonlinear, we validate the NLSE and the DW–ZF soliton solution. Then, we demonstrate that the NLSE can adequately describe the spontaneous generation of solitary zonal structures in the mHME, which takes place in the nonlinear stage of the MI, but only when the amplitude of the primary DW is relatively small. Otherwise, stationary zonal structures are formed instead.

Next, we consider the MI of incoherent DW spectra, using the recently developed WM model. We show that DW spectra that consist of quasi-monochromatic drifton beams can also be modulationally unstable to the formation of solitary zonal structures, but only when the beam intensity is moderate. At higher intensity, the zonal structures become stationary, similar to the case with coherent DWs. Meanwhile, due to the stabilizing effect of the finite beam width, the system becomes stable to modulations at lower beam intensity.

In addition, using the WM formulation, we compare the solitary zonal structures formed via the MI with the DW–ZF solitons in phase space. This approach enables extraction of information that can be obscure in configuration space, especially when the DW spectrum has multiple quasi-monochromatic peaks. (As a data analysis tool alone, it can be straightforwardly applied to systems that are more complicated than the mHME). The phase-space distributions of DW quanta show common features of quasi-monochromatic drifton condensates, which suggests that the MI-induced solitary zonal structures are essentially the DW–ZF solitons. These structures cannot be described by wave-kinetic models that are based on the ray approximation, which neglect critical quantumlike effects such as diffraction. In contrast, the WM model retains these effects, subsumes both the NLSE and the wave kinetic models, and hence can support these solitary zonal structures.

It is worthwhile to comment on the relevance of our results to Rossby-wave turbulence in geophysics. In Ref. Dewar and Abdullatif (2007), the envelope dynamics within the oHME is also discussed, and the equation of state is more complicated than our Eq. (10); namely, it is not a local but an integral equation. Hence, the applicability of our results to the Charney equation remains to be further investigated. Meanwhile, our results are readily applicable to the barotropic vorticity equation, and the only adaption needed is to replace $\hat{\bar{p}}^{2}$ with $\hat{{p}}^{2}\doteq\hat{p}_{x}^{2}+\hat{p}_{y}^{2}$, and similarly for ${\bar{p}}^{2}$ and ${\bar{k}}^{2}$.

In the future, we plan to investigate the effect of background shear flows on the DW–ZF solitons, motivated by the radially propagating coherent structures recently identified in gyrokinetic simulations of subcritical plasmas van Wyk et al. (2016, 2017); McMillan et al. (2018). In particular, the structures in Ref. McMillan et al. (2018) [Fig. 3(b) therein] seem close to quasi-monochromatic, much similar to our Fig. 3(1-a). The (normalized) sizes and amplitudes of those structures are also comparable to those of the DW–ZF solitons discussed in this paper. However, to properly account for subcriticality would be a challenge, and it is possible that one needs to resort to more sophisticated models with primary instabilities, such as the (modified) Hasegawa–Wakatani equations Hasegawa and Wakatani (1983); Numata et al. (2007).

Acknowledgements.

We thank L. Chen and F. Zonca for bringing Ref. Guo et al. (2009) to our attention, which led to this work. The research was supported by the U.S. Department of Energy under Contract No. DE-AC02-09CH11466.

Appendix A Limitations of wave-kinetic models

Wave-kinetic models of DW–ZF dynamics invoke the ray (geometrical-optics) approximation, i.e., assume that the DW wavelength is negligible compared with the ZF wavelength. In this case, the WM model (13) reduces to

[TABLE]

with the Poisson bracket $\{A,B\}=(\partial_{\mathbf{x}}A)\cdot(\partial_{\mathbf{p}}B)-(\partial_{\mathbf{p}}A)\cdot(\partial_{\mathbf{x}}B)$. The Hermitian and anti-Hermitian parts of the Hamiltonian are given by, respectively Ruiz et al. (2016),

[TABLE]

This model was first derived as the geometrical-optics limit of the CE2 equations Parker (2016). Following Refs. Zhu et al. (2018a, c), we refer to it as the improved wave-kinetic equation (iWKE). The dispersion relation of the MI in the iWKE reads Parker (2016); Zhu et al. (2018a)

[TABLE]

By substituting the Wigner function (17) into Eq. (26), we obtain the iWKE dispersion relation of the MI of a monochromatic DW,

[TABLE]

The iWKE dispersion relation (27) is plotted in Fig. 2 for comparison with the NLSE dispersion relation (12) and the WM dispersion relation (18). To better understand the discrepancy between the them, let us also consider the traditional wave-kinetic equation (tWKE) Smolyakov and Diamond (1999). The tWKE differs from the iWKE by further neglecting the high-order derivatives of $U$ in the Hamiltonian (25), such that

[TABLE]

Accordingly, the tWKE dispersion relation for the MI can be obtained by simply neglecting the factor ${q^{2}}/{\bar{k}^{2}}$ in Eq. (27):

[TABLE]

This dispersion relation diverges at large $q$, which is typical for the tWKE. Such divergence can give rise to unphysical grid-scale ZFs in numerical simulations, as shown in Refs. Parker (2016); Ruiz et al. (2016), where more detailed discussions on the differences between the tWKE, the iWKE, and the WM (CE2) model can be found.

Despite this caveat, the tWKE dispersion relation (29) warrants direct comparison with the NLSE dispersion relation (12), since both models neglect the same high-oder derivatives of $U$. The difference is the $\chi^{2}q^{4}/4$ term. This term is due to diffraction, which is not included in the ray approximation that the tWKE (and iWKE) is based on, but retained in the quasi-optical approximation that the NLSE invokes. Likewise, it is the balance of diffraction and self-focusing that determines the size of the solitary zonal structures discussed in Sec. V. In summary, because of the absence of full-wave (quantumlike) effects such as diffraction, wave-kinetic models based on the ray approximation do not properly capture the MI or the solitary zonal structures discussed in the main text.