Nonlinear saturation and oscillations of collisionless zonal flows

TL;DR

This paper develops a semi-analytic theory explaining the transition between oscillations and saturation of collisionless zonal flows in drift-wave turbulence, based on phase-space analysis of drift-wave quanta within the geometrical-optics approximation.

Contribution

It introduces a control parameter N that predicts whether ZFs oscillate or saturate, supported by analytical derivations and simulations, offering new insights into collisionless ZF dynamics.

Findings

N determines transition between oscillations and saturation

Passing drifton trajectories cause oscillations at N<<1

Phase mixing leads to saturation at N≥1

Abstract

In homogeneous drift-wave (DW) turbulence, zonal flows (ZFs) can be generated via a modulational instability (MI) that either saturates monotonically or leads to oscillations of the ZF energy at the nonlinear stage. This dynamics is often attributed as the predator--prey oscillations induced by ZF collisional damping; however, similar dynamics is also observed in collisionless ZFs, in which case a different mechanism must be involved. Here, we propose a semi-analytic theory that explains the transition between the oscillations and saturation of collisionless ZFs within the quasilinear Hasegawa--Mima model. By analyzing phase-space trajectories of DW quanta (driftons) within the geometrical-optics (GO) approximation, we argue that the parameter that controls this transition is $N\sim\gamma_{{\rm MI}}/\omega_{{\rm DW}}$, where $\gamma_{{\rm MI}}$ is the MI growth rate and $\omega_{{\rm DW}}$ is the linear DW frequency. We argue that at $N\ll1$, ZFs oscillate due to the presence of so-called passing drifton trajectories, and we derive an approximate formula for the ZF amplitude as a function of time in this regime. We also show that at $N\gtrsim1$, the passing trajectories vanish and ZFs saturate monotonically, which can be attributed to phase mixing of higher-order sidebands. A modification of $N$ that accounts for effects beyond the GO limit is also proposed. These analytic results are tested against both quasilinear and fully-nonlinear simulations. They also explain the earlier numerical results by Connaughton $\textit{et al}$. [J. Fluid Mech. $\textbf{654}$, 207 (2010)] and Gallagher $\textit{et al}$. [Phys. Plasmas $\textbf{19}$, 122115 (2012)] and offer a revised perspective on what the control parameter is that determines the transition from the oscillations to saturation of collisionless ZFs.

Equations

$$\frac{\partial w}{\partial t}+(\hat{\mathbf{z}}\times \nabla \phi )\cdot \nabla w+\beta \frac{\partial \phi }{\partial x}=0,$$

$$w=({\nabla }^2-L_{\mathrm{D}}^{-2}\hat{\alpha })\phi .$$

$$\frac{\partial {\phi }_{\mathbf{k}}}{\partial t}=i{\omega }_{\mathbf{k}}{\phi }_{\mathbf{k}}+\frac{1}{2}\sum _{{\mathbf{k}}_1,{\mathbf{k}}_2}T(\mathbf{k},{\mathbf{k}}_1,{\mathbf{k}}_2){\phi }_{{\mathbf{k}}_1}{\phi }_{{\mathbf{k}}_2}{\delta }_{\mathbf{k},{\mathbf{k}}_1+{\mathbf{k}}_2}.$$

$${\omega }_{\mathbf{k}}\doteq -\frac{\beta k_x}{k_{\mathrm{D}}^2}$$

$$T(\mathbf{k},{\mathbf{k}}_1,{\mathbf{k}}_2)\doteq -\frac{k_{1,\mathrm{D}}^2-k_{2,\mathrm{D}}^2}{k_{\mathrm{D}}^2}({\mathbf{k}}_1\times {\mathbf{k}}_2)\cdot \hat{\mathbf{z}}$$

$$k_{\mathrm{D}}^2\doteq |\mathbf{k}{|}^2+L_{\mathrm{D}}^{-2}{\alpha }_{\mathbf{k}},$$

$$k_{n,\mathrm{D}}^2\doteq |{\mathbf{k}}_n{|}^2+L_{\mathrm{D}}^{-2}{\alpha }_{{\mathbf{k}}_n}.$$

$${\partial }_t{\Phi }_{\mathbf{p}}=T(\mathbf{p},\mathbf{q},{\mathbf{p}}_{-}){\Phi }_{\mathbf{q}}{\Phi }_{{\mathbf{p}}_{-}}{e}^{i{\Delta }_{-}t}$$

$$+T(\mathbf{p},-\mathbf{q},{\mathbf{p}}_{+}){\Phi }_{\mathbf{q}}^{\ast }{\Phi }_{{\mathbf{p}}_{+}}{e}^{i{\Delta }_{+}t},$$

$${\partial }_t{\Phi }_{\mathbf{q}}=T(\mathbf{q},\mathbf{p},-{\mathbf{p}}_{-}){\Phi }_{\mathbf{p}}{\Phi }_{{\mathbf{p}}_{-}}^{\ast }{e}^{-i{\Delta }_{-}t}$$

$$+T(\mathbf{q},-\mathbf{p},{\mathbf{p}}_{+}){\Phi }_{\mathbf{p}}^{\ast }{\Phi }_{{\mathbf{p}}_{+}}{e}^{i{\Delta }_{+}t},$$

$${\partial }_t{\Phi }_{{\mathbf{p}}_{-}}=T({\mathbf{p}}_{-},\mathbf{p},-\mathbf{q}){\Phi }_{\mathbf{p}}{\Phi }_{\mathbf{q}}^{\ast }{e}^{-i{\Delta }_{-}t},$$

$${\partial }_t{\Phi }_{{\mathbf{p}}_{+}}=T({\mathbf{p}}_{+},\mathbf{p},\mathbf{q}){\Phi }_{\mathbf{p}}{\Phi }_{\mathbf{q}}{e}^{-i{\Delta }_{+}t},$$

$$\displaystyle\left(\begin{array}[]{c}\Phi_{\boldsymbol{p}}\\ \Phi_{\boldsymbol{q}}\\ \Phi_{\boldsymbol{p}_{+}}\\ \Phi_{\boldsymbol{p}_{-}}\end{array}\right)=\left(\begin{array}[]{c}\Phi_{0}\\ 0\\ 0\\ 0\end{array}\right)+\varepsilon\left(\begin{array}[]{c}0\\ A_{\boldsymbol{q}}e^{-i\Omega_{\boldsymbol{q}}t}\\ A_{\boldsymbol{p}_{+}}e^{-i\Omega_{\boldsymbol{p}_{+}}t}\\ A_{\boldsymbol{p}_{-}}e^{-i\Omega_{\boldsymbol{p}_{-}}t}\end{array}\right),$$

$$\begin{gathered} {\Omega }_{\mathbf{q}}+\frac{|{\Phi }_0{|}^{2}T(\mathbf{q},-\mathbf{p},{\mathbf{p}}_{+})T({\mathbf{p}}_{+},\mathbf{q},\mathbf{p})}{{\Delta }_{+}+{\Omega }_{\mathbf{q}}} \\ -\frac{|{\Phi }_0{|}^{2}T(\mathbf{q},-{\mathbf{p}}_{-},\mathbf{p})T({\mathbf{p}}_{-},-\mathbf{q},\mathbf{p})}{{\Delta }_{-}-{\Omega }_{\mathbf{q}}}=0. \end{gathered}$$

$${\Delta }_{\pm }\equiv \Delta \doteq \frac{\beta p{q}^2}{p_{\mathrm{D}}^2(p_{\mathrm{D}}^2+q^2)},$$

$$p_{\mathrm{D}}^2\doteq L_{\mathrm{D}}^{-2}+p^2.$$

$${\gamma }_{\mathrm{MI}}^2=\frac{2|{\Phi }_0{|}^2{p}^2{q}^2}{(1+\alpha L_{\mathrm{D}}^{-2}{q}^{-2})}\frac{{\delta }^{\prime }-1}{\delta +1}-{\Delta }^2.$$

$$\delta \doteq \frac{p_{\mathrm{D}}^2}{q^2},{\delta }^{\prime }\doteq \frac{p_{\mathrm{D}}^2-\alpha L_{\mathrm{D}}^{-2}}{q^2},$$

$$q^2<p_{\mathrm{D}}^2-\alpha L_{\mathrm{D}}^{-2}.$$

$$\frac{\partial \tilde{w}}{\partial t}+U\frac{\partial \tilde{w}}{\partial x}+\left[\beta -\left(\frac{{\partial }^2}{\partial y^2}-\alpha L_{\mathrm{D}}^{-2}\right)U\right]\frac{\partial \tilde{\phi }}{\partial x}=f_{\mathrm{NL}},$$

$$\left[1-\alpha L_{\mathrm{D}}^{-2}{\left(\frac{\partial }{\partial y}\right)}^{-2}\right]\frac{\partial U}{\partial t}=-\frac{\partial }{\partial y}⟨{\tilde{v}}_x{\tilde{v}}_y⟩.$$

$$U=u(t)\cos qy,$$

$$\begin{gathered} \frac{\partial {\phi }_m}{\partial t}=i{\omega }_m{\phi }_m+\frac{ipu}{2}\left(\frac{q^2+\alpha L_{\mathrm{D}}^{-2}-k_{m+1,\mathrm{D}}^2}{k_{m,\mathrm{D}}^2}\right){\phi }_{m+1} \\ +\frac{ipu}{2}\left(\frac{q^2+\alpha L_{\mathrm{D}}^{-2}-k_{m-1,\mathrm{D}}^2}{k_{m,\mathrm{D}}^2}\right){\phi }_{m-1}, \end{gathered}$$

$${\omega }_m\doteq -\frac{\beta p}{k_{m,\mathrm{D}}^2},k_{m,\mathrm{D}}^2\doteq p_{\mathrm{D}}^2+(mq{)}^2.$$

$$W(y,\mathbf{k},t)=\int d^{2}s{e}^{-i\mathbf{k}\cdot \mathbf{s}}⟨\tilde{w}(\mathbf{x}+\frac{\mathbf{s}}{2},t)\tilde{w}(\mathbf{x}-\frac{\mathbf{s}}{2},t)⟩.$$

$${\partial }_{t}W(y,\mathbf{k},t)=\{\mathcal{H},W\}+2\Gamma W.$$

$$\{A,B\}\doteq \frac{\partial A}{\partial \mathbf{x}}\cdot \frac{\partial B}{\partial \mathbf{k}}-\frac{\partial A}{\partial \mathbf{k}}\cdot \frac{\partial B}{\partial \mathbf{x}},$$

$$\mathcal{H}(y,\mathbf{k},t)=k_x\left[u(t)\cos qy-\frac{\beta +(q^2+\alpha L_{\mathrm{D}}^{-2})u(t)\cos qy}{k_{\mathrm{D}}^2}\right],$$

$$\Gamma (y,\mathbf{k},t)=-\frac{k_x{k}_{y}q(q^2+\alpha L_{\mathrm{D}}^{-2})u(t)\sin qy}{k_{\mathrm{D}}^4}.$$

$$W_{\lambda }(\mathbf{k},t)\doteq \int W(y,\mathbf{k},t)\exp (-i\lambda y)dy,$$

$$\begin{gathered} {\partial }_t{W}_{\lambda }(p,k_y,t)=i\beta p\left(\frac{1}{{\kappa }_{+\lambda }^2}-\frac{1}{{\kappa }_{-\lambda }^2}\right)W_{\lambda ,0} \\ -\frac{iup}{2}{Q}_{\lambda -2q}{W}_{\lambda -q,-q}+\frac{iup}{2}{Q}_{2q-\lambda }{W}_{\lambda -q,+q} \\ -\frac{iup}{2}{Q}_{\lambda +2q}{W}_{\lambda +q,+q}+\frac{iup}{2}{Q}_{-2q-\lambda }{W}_{\lambda +q,-q}, \end{gathered}$$

$$W_{a,b}\doteq W_a\left(p,k_y+\frac{b}{2},t\right),{\kappa }_a^2\doteq p_{\mathrm{D}}^2+{\left(k_y+\frac{a}{2}\right)}^2,$$

$$Q_a\doteq 1-\frac{\alpha L_{\mathrm{D}}^{-2}+q^2}{{\kappa }_a^2}.$$

Full text

I Introduction

Zonal flows (ZFs) are banded sheared flows that can spontaneously emerge from drift-wave (DW) turbulence in magnetized plasmas (Lin98, ; Diamond05, ; Rogers00, ; Jenko00, ) and, similarly, from Rossby-wave turbulence in the atmospheres of rotating planets (Vasavada05, ). They are considered important as regulators of turbulent transport and thus have been studied actively by many researchers. One mechanism of the ZF generation is the secondary (or zonostrophic) instability (Rogers00, ; Smolyakov00a, ; Smolyakov00b, ; Strintzi07, ; Srinivasan2012, ; Parker2013, ; Parker2014, ; Parker14Thesis, ; St-Onge17a, ; Marston16, ), and the modulational instability (MI) as a special case when the DW is monochromatic (Qi19, ; Chen00, ; Connaughton10, ; Gallagher12, ; Connaughton15, ; Champeaux01, ; Koshkarov16, ; Lashkin08, ; Jenko06, ; Shukla02, ; Manfredi01, ; Lashmore-Davies01, ; Gill74, ; Mahanti81, ; St-Onge17b, ). The linear stage of the MI is generally understood, but the dynamics of ZFs at the nonlinear stage is not sufficiently explored. Some progress in this area has been made by applying quasilinear (QL) models (Sec. II.4), such as the second-order cumulant expansion theory (CE2) (Srinivasan2012, ; Parker2013, ; Parker2014, ; Parker14Thesis, ; St-Onge17a, ; Marston16, ), or the stochastic structural stability theory (SSST) (Constantinou18, ; Farrell07, ; Bakas15, ); however, those are not particularly intuitive. A more intuitive paradigm was proposed based on a simpler QL model known as the wave-kinetic equation (WKE) (Smolyakov00a, ; Smolyakov00b, ; Ruiz16, ; Ruiz18, ; Zhu18a, ; Zhu18b, ; Parker16, ; Smolyakov99, ; Kaw2001, ). The WKE treats DW turbulence as a collection of DW quanta (“driftons”), for which the ZF velocity serves as a collective field. Then, the ZF–DW interactions can be viewed as the predator–prey dynamics, which can result in either monotonic or oscillatory energy exchange between the ZFs and DWs (Malkov01, ; Miki12, ; Diamond94, ; Kim03, ; Berionni11, ). Such dynamics is indeed observed in simulations (Kobayashi15, ). However, the existing paradigm substantially relies on collisional damping of ZFs, whereas simulations indicate that nonlinear saturation and oscillations are also possible when the ZF is collisionless (Connaughton10, ; Gallagher12, ; Manfredi01, ; Lashmore-Davies01, ; Mahanti81, ). This leads to the question of how this dynamics can be explained within a collisionless theory.

Here, we propose a semi-analytic theory that explains the transition between the oscillations and saturation of collisionless ZFs within the QL Hasegawa–Mima model. By analyzing phase-space trajectories of DW quanta (driftons) within the geometrical-optics (GO) approximation, we argue that the parameter that controls this transition is $N\sim\gamma_{{\rm MI}}/\omega_{{\rm DW}}$, where $\gamma_{{\rm MI}}$ is the MI growth rate and $\omega_{{\rm DW}}$ is the linear DW frequency. We argue that at $N\ll 1$, ZFs oscillate due to the presence of so-called passing drifton trajectories, and we derive an approximate formula for the ZF amplitude as a function of time in this regime. In doing so, we also extend the applicability of the popular “four-mode truncation” (4MT) model (Chen00, ; Connaughton10, ; Gallagher12, ; Connaughton15, ; Champeaux01, ; Koshkarov16, ; Lashkin08, ; Jenko06, ; Shukla02, ; Manfredi01, ; Lashmore-Davies01, ; Mahanti81, ; Gill74, ; St-Onge17b, ), which is commonly used for the linear stage, to nonlinear ZF–DW interactions. We also show that at $N\gtrsim 1$, the passing trajectories vanish and ZFs saturate monotonically, which can be attributed to phase mixing of higher-order sidebands. A modification of $N$ that accounts for effects beyond the GO limit is also proposed. These analytic results are tested against both QL and fully-nonlinear (NL) simulations. They also explain the earlier numerical results by Connaughton et al. (Connaughton10, ) and Gallagher et al. (Gallagher12, ) and offer a revised perspective on what the control parameter is that determines the transition from oscillations to saturation of collisionless ZFs.

Our paper is organized as follows. In Sec. II, we introduce the basic equations, including the Hasegawa–Mima equation, the 4MT, the QL approximation, the Wigner function, and the WKE that we use later on. In Sec. III, we discuss the nonlinear stage of the MI using both the 4MT description and the WKE, and we also propose our control parameter $N$. In Sec. IV, we compare our results with those in Refs (Connaughton10, ; Gallagher12, ). In Sec. VI we summarize our main conclusions.

II Basic equations

II.1 Hasegawa–Mima equation

Both electrostatic DW turbulence in plasmas and Rossby-wave turbulence in the atmospheres of rotating planets are often modeled by the Hasegawa–Mima equation (HME) (Hasegawa77, ):

[TABLE]

(In geophysics, it is also known as the Obukhov–Charney equation (DolzhanskyBook, ).) Here, the geophysics coordinate convention is used to simplify comparisons with the earlier relevant studies (Sec. V). The HME describes wave turbulence on a two-dimensional plane $(x,y)$, where ZFs develop along the $x$-direction, and $\nabla^{2}\doteq\partial_{x}^{2}+\partial_{y}^{2}$ is the Laplacian. In the plasma-physics context, the system is assumed to be immersed in a uniform magnetic field along the $z$ axis, the plasma is assumed to have a constant density gradient along the $y$ axis, $\beta$ is a scalar constant proportional to this gradient, $L_{{\rm D}}$ is the ion-sound radius, and $\varphi$ is the perturbation of the electrostatic potential. In the geophysical context, the constant $\beta$ is proportional to the latitudinal gradient of the vertical rotation frequency, $L_{{\rm D}}$ is the deformation radius, and $\varphi$ is the stream function.

The operator $\hat{\alpha}$ is an identity operator in the original HME (oHME). The so-called modified HME (mHME), which is also known as “generalized” (Ruiz16, ; Krommes00, ) or “extended” HME (Connaughton15, ; Gallagher12, ), uses $\hat{\alpha}f\doteq f-\langle f\rangle$, where $f=f(x,y,t)$ is any field quantity and $\langle f\rangle\doteq\int_{0}^{L_{x}}fdx/L_{x}$ denotes the zonal average, where $L_{x}$ is the system length in the $x$-direction. (For more details, see, for example, Ref. (St-Onge17b, ; Hammett93, ).) Below, we consider both the oHME and mHME and treat them on the same footing.

II.2 Fourier decomposition and the 4MT

It is common to approach the dynamics of $\varphi$ in the Fourier representation, $\varphi=\sum_{\boldsymbol{k}}\varphi_{\boldsymbol{k}}(t)\exp(i\boldsymbol{k}\cdot\boldsymbol{x})$, where $\boldsymbol{k}=(k_{x},k_{y})$ and $\boldsymbol{x}=(x,y)$. This leads to the following equation for $\varphi_{\boldsymbol{k}}$:

[TABLE]

Here, $\delta_{\boldsymbol{k},\boldsymbol{k}_{1}+\boldsymbol{k}_{2}}$ equals one only if $\boldsymbol{k}=\boldsymbol{k}_{1}+\boldsymbol{k}_{2}$ and is zero otherwise,

[TABLE]

is the linear DW frequency, and

[TABLE]

are the coefficients that govern the nonlinear mode coupling. Also,

[TABLE]

and $k_{n,{\rm D}}^{2}$ (where $n=1,2$) are similarly defined as

[TABLE]

And finally, $\alpha_{\boldsymbol{k}}$ is the Fourier representation of the operator $\hat{\alpha}$; namely, for the oHME $\alpha_{\boldsymbol{k}}$ is unity, and for the mHME $\alpha_{\boldsymbol{k}}$ equals zero if $k_{x}=0$ and equals unity if $k_{x}\neq 0$.

It is also common to introduce the 4MT, which is a truncation of the system (3) that retains only four Fourier harmonics, namely, those with wave vectors $\boldsymbol{p}$, $\boldsymbol{q}$, and $\boldsymbol{p}_{\pm}\doteq\boldsymbol{p}\pm\boldsymbol{q}$. (Since $\varphi$ is real, one has $\varphi_{-\boldsymbol{k}}=\varphi_{\boldsymbol{k}}^{*}$; hence, the harmonics with wave vectors $-\boldsymbol{p}$, $-\boldsymbol{q}$, and $-\boldsymbol{p}_{\pm}$ are included too.) Then, for $\Phi_{\boldsymbol{k}}\doteq\varphi_{\boldsymbol{k}}\exp(-i\omega_{\boldsymbol{k}}t)$, Eq. (3) gives

[TABLE]

where $\Delta_{\pm}=\omega_{\boldsymbol{p}}\pm\omega_{\boldsymbol{q}}-\omega_{\boldsymbol{p}_{\pm}}$.

II.3 Modulational instability

Suppose a perturbation on a primary wave $\Phi_{\boldsymbol{p}}=\Phi_{0}$,

[TABLE]

where $\varepsilon$ is small. Then, the linearized Eq. (8) gives $\Omega_{\boldsymbol{p}\pm}=\Delta_{\pm}\pm\Omega_{\boldsymbol{q}}$ together with the following dispersion relation:

[TABLE]

(The derivation can be found, for example, in Ref. (Gallagher12, ).) This equation can have a complex solution for $\Omega_{\boldsymbol{q}}$ with $\text{Im}\,\Omega_{\boldsymbol{q}}>0$, which signifies the presence of the MI. As shown in Refs. (Connaughton10, ; Gill74, ), the 4MT is indeed often a good approximation at the linear stage of the MI.

In the following, we restrict our discussions to the case when $\boldsymbol{p}=(p,0)$ and $\boldsymbol{q}=(0,q)$. In this case, $\omega_{\boldsymbol{q}}=0$ and

[TABLE]

where

[TABLE]

Then, one finds that $\Omega_{\boldsymbol{q}}^{2}$ is real, and hence the MI has a positive growth rate $\gamma_{{\rm MI}}$ if $\Omega_{\boldsymbol{q}}^{2}$ is negative; namely, $\gamma_{{\rm MI}}^{2}\doteq-\Omega_{\boldsymbol{q}}^{2}$, which can be explicitly written as

[TABLE]

Here,

[TABLE]

with $\alpha=1$ for the oHME and $\alpha=0$ for the mHME. Notably, this implies that the mHME gives much larger growth rates than the oHME does. Also, a necessary condition for the modulation to be unstable is $\delta^{\prime}>1$, which requires

[TABLE]

II.4 Quasilinear approximation

In order to describe the nonlinear stage of the MI analytically, we proceed as follows. Let us decompose the field quantities into the zonal-averaged part and the fluctuation part, $f=\langle f\rangle+\tilde{f}$. Then, Eqs. (1) and (2) become (Parker14Thesis, )

[TABLE]

Here, $U(y,t)\doteq-\partial_{y}\langle\varphi\rangle$ is the ZF velocity, $\tilde{\boldsymbol{v}}\doteq\hat{\boldsymbol{z}}\times\nabla\tilde{\varphi}$ is the fluctuation velocity, $\partial_{y}^{-2}$ is an operator that in the wave-vector (Fourier) representation is simply a multiplication by $-k_{y}^{-2}$, and $f_{{\rm NL}}\doteq\tilde{\boldsymbol{v}}\cdot\nabla\tilde{w}-\langle\tilde{\boldsymbol{v}}\cdot\nabla\tilde{w}\rangle$ describes self-interactions of DWs, or eddy–eddy interactions. We shall simplify the problem by ignoring these interactions, i.e., by adopting $f_{{\rm NL}}=0$. This is the commonly-used QL approximation, which often yields an adequate description of ZF–DW interactions (Parker14Thesis, ). (We shall also discuss the applicability of this approximation in Sec. IV.2). Notably, once the QL approximation is adopted, the 4MT model becomes the exact description of the linear MI.

We also assume, for simplicity, that the ZF is sinusoidal and non-propagating (assuming $\gamma_{{\rm MI}}$ is real),

[TABLE]

where for clarity we choose the origin on the $y$ axis such that $y=0$ corresponds to the maximum of $U$. (Propagating zonal structures are studied in detail in our Ref. (soliton, ).) The assumption of spatially-monochromatic ZF holds approximately if the ZF’s second and higher harmonics do not outpace the fundamental harmonic during the linear MI. Hence, the ansatz (30) implies that the value of $q$ is close to the one that maximizes $\gamma_{{\rm MI}}$. Then, the DW field can be represented as $\tilde{\varphi}={\rm Re}\sum_{m=-\infty}^{\infty}\varphi_{m}(t)\exp(ipx+imqy)$, where $\varphi_{0}(t=0)=\Phi_{0}$ is the primary-wave amplitude, and the DW equation (28) becomes

[TABLE]

where

[TABLE]

Note that the definition of $k_{m,D}$ is consistent with that given by Eq. (7).

II.5 Wigner function and WKE

Equation (31) describes the mode-coupling among different Fourier modes of DWs due to the ZF. Now, we seek to interpret this equation in terms of the drifton phase-space dynamics. Consider the zonal-averaged Wigner function of DWs,

[TABLE]

It can be understood as the spectral representation of the DW two-point correlation function, such as the one used in the CE2. It can also be viewed as the quasiprobability of the drifton distribution in the $(y,\boldsymbol{k})$ space, where the prefix “quasi” reflects the fact that $W$ is not necessarily positive-definite. That said, it becomes positive-definite in the GO limit defined as (i) $\alpha L_{{\rm D}}^{-2}+q^{2}\ll L_{{\rm D}}^{-2}+p^{2}$ and (ii) $\partial_{k_{y}}W\ll q^{-1}W$, when it can be considered as the true distribution function of driftons (Cartwright76, ). Then, the Wigner function can be shown to satisfy the following partial differential equation (Ruiz16, ; Parker16, ):

[TABLE]

Here, $\{\cdot,\cdot\}$ is the Poisson bracket,

[TABLE]

and the functions $\mathcal{H}$ and $\Gamma$ can be interpreted as the drifton Hamiltonian and the drifton dissipation rate, respectively. Specifically, they are given by (Ruiz16, ; Zhu18b, )

[TABLE]

Equation (34) is the WKE that describes the phase-space dynamics of driftons (DW quanta). They obey Hamilton’s equations and move on constant-$\mathcal{H}$ surfaces in the phase-space $(y,k_{y})$ if the ZF is stationary. (The $x$-momentum $k_{x}=p$ serves only as a parameter.) In the GO limit, $\Gamma$ is small and unimportant for the discussions below. However, in general, keeping $\Gamma$ is necessary to ensure that the WKE preserves the conservation of the fundamental integrals of the HME (Parker16, ).

II.6 Evolution of $\boldsymbol{W}$ beyond the GO approximation

Beyond the GO approximation, the Wigner function satisfies a pseudo-differential equation known as the Wigner–Moyal equation (Ruiz16, ). But here, we use a somewhat different formulation. Following the procedure in Ref. (Ruiz16, ), we consider the spectrum of $W$,

[TABLE]

which is a function in the double-momentum space $(\lambda,k_{y})$. Then, Eq. (31) is transformed into

[TABLE]

where we introduced

[TABLE]

and

[TABLE]

Equation (38) is equivalent to Eq. (31), but describes the dynamics in the double-momentum space. The DW momentum flux can also be expressed through $W_{\lambda}$, whose gradient drives the Fourier component of the ZF, $U_{\lambda=q}\doteq\pi u(t)$. Specifically, from Eq. (29), the following equation is obtained:

[TABLE]

III Nonlinear stage of the modulation instability

III.1 4MT description

In order to describe the nonlinear stage of the MI, let us first consider this instability within the 4MT model. We shall show that the 4MT has an exact analytic solution that not only gives the linear growth rate (25), but also predicts the reversal of the ZF growth. We shall also propose a toy-model modification of the 4MT that qualitatively explains the transition from ZF oscillations to saturation.

For the MI, the initial condition is $\tilde{\varphi}(\boldsymbol{x},t=0)=\Phi_{0}\exp(ipx)+{\rm c.c.}$, which corresponds to a delta function in the double-momentum space [Eq. (33)],

[TABLE]

According to Eq. (38), $W_{\lambda}(\boldsymbol{k},t)$ remains delta-shaped also at $t>0$, so we search it in the form

[TABLE]

where the factor $\delta(k_{x}-p)$ is omitted, and $W_{\lambda}(-\boldsymbol{k},t)=W_{\lambda}(\boldsymbol{k},t)$. The 4MT model (Sec. II.2) corresponds to keeping the following nine peaks:

[TABLE]

Here, $a$, $b$, $c$, and $d$ are real, and the initial value of $a(t)$ is $a_{0}$. Using Eq. (38), one obtains

[TABLE]

Here, $\Delta$, $\delta$, and $\delta^{\prime}$ are defined by Eqs. (23) and (26). Also, the ZF amplitude is governed by [Eq. (40)]

[TABLE]

where we introduced

[TABLE]

Note that DWs with $k_{x}=p$ and $k_{x}=-p$ contribute equally to $\dot{u}$.

In the following, we simplify the above equations and obtain the time-evolution of $u$. First, notice that

[TABLE]

Thus, $d$ can be expressed through $a$:

[TABLE]

where $a_{0}$ is given by Eq. (42). Next, we introduce dimensionless variables

[TABLE]

where we choose $B=C=\beta\epsilon V$ and

[TABLE]

Then, Eqs. (44) and (45) become

[TABLE]

where $\bar{a}_{0}=a_{0}/A$. This leads to

[TABLE]

where

[TABLE]

Assuming infinitesimally small initial $\bar{u}$, we have

[TABLE]

Finally, since

[TABLE]

one has that $\bar{f}-\bar{u}^{2}/2$ is conserved; namely,

[TABLE]

is a constant. (Having $g^{2}<0$ is possible but does not lead to an instability; see below.) Hence, we obtain that $\bar{u}$ satisfies a nonlinear-oscillator equation

[TABLE]

The effective potential $\Theta(\bar{u})$ is plotted in Fig. 1(a). (A small nonzero initial value of $\bar{u}$ causes a slight alteration of $\Theta$, but the qualitative picture remains the same.) The steady-state solution $\bar{u}=0$ is unstable if $g^{2}>0$, which signifies the presence of a linear instability (namely, the MI) with the growth rate

[TABLE]

which is in agreement with Eq. (25). Beyond the linear regime, i.e., when $\bar{u}^{4}$ is no longer negligible, Eq. (54) can also be integrated exactly, yielding

[TABLE]

This solution corresponds to the initial condition $\bar{a}(\tau\to-\infty)=\bar{a}_{0}$ and $\bar{u}(\tau\to-\infty)=0$, and the origin on the time axis is chosen such that the ZF attains the maximum amplitude at $\tau=0$. (A comparison between this analytic solution and simulations is presented in Sec. IV.1) In our original variables, this maximum amplitude is given by

[TABLE]

or more explicitly,

[TABLE]

The 4MT dynamics is numerically illustrated in Fig. 1(b). Unlike the exact solution (56), a finite initial perturbation of $\bar{u}$ results in oscillations with a finite period [Fig. 1(b), blue curve with circles]. We also propose a toy model to illustrate the transition from oscillations to saturation of $\bar{u}$, which adds an ad hoc damping of the sidebands to mimic their coupling to higher harmonics. Specifically, we replace Eq. (49) with

[TABLE]

where $\nu$ is some positive constant. Then, as we increase $\nu$, a gradual transition from oscillations to saturation of $\bar{u}$ is observed as shown in Fig. 1(b). In particular, at very large $\nu$ such that $\nu\gg d/d\tau$, one has

[TABLE]

and Eq. (49) gives

[TABLE]

This equation also has an exact solution

[TABLE]

which describes monotonic saturation of the ZF [Fig. 1(b), green curve with squares]. This shows that, qualitatively, the ZF saturation can be explained as a result of phase mixing that occurs due to the primary-wave coupling to higher harmonics. For quantitative predictions, a more rigorous approach is proposed below based on exploring the drifton phase-space dynamics.

III.2 Control parameter $\boldsymbol{N}$ in the GO limit

Let us study the nonlinear stage of the MI using the WKE [Eq. (34)], which naturally accounts for all DW sidebands. As a starting point, we appeal to the findings from Ref. (Zhu18b, ), where three types of drifton phase-space trajectories were identified: trapped, passing, and runaway. (The last one corresponds to driftons leaving to infinity along the $k_{y}$ axis while retaining finite $y$.) It was also found in Ref. (Zhu18b, ) that the drifton phase-space topology for a sinusoidal ZF can vary depending on how the ZF amplitude relates to the following two critical values:

[TABLE]

The GO limit corresponds to $u_{c,1}\ll u_{c,2}$. At $u<u_{c,1}$ (regime 1), trajectories of all three types are possible; at $u_{c,1}\leq u<u_{c,2}$ (regime 2), passing trajectories vanish; and finally, at $u>u_{c,2}$ (regime 3), trapped trajectories also vanish.

Figures 2 and 3 show the corresponding evolution of $W_{\lambda}$ calculated numerically from Eq. (38) at fixed ZF velocity, namely, $U=u\cos qy$ with constant $u$. Under the GO assumption, the mHME and oHME exhibit similar dynamics, so for clarity, we consider only the mHME ($\alpha=0$). Since $u$ is assumed stationary, Eq. (38) is linear in $W$, and the initial condition is chosen as [Eq. (43)]

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta. Figure 2 corresponds to $u<u_{c,1}$ (regime 1). In this case, the distribution $W(y,\boldsymbol{k})$ in the phase space is confined to small $|k_{y}|$ regions enclosed by the contours of passing and trapped trajectories, and no driftons reside on runaway trajectories (transport along $k_{y}$ is suppressed). The corresponding $\bar{W}_{m,n}$ are vanishingly small at large $m$ and $n$; i.e., only a finite number of harmonics are coupled. In contrast, Fig. 3 corresponds to $u_{c,1}<u<u_{c,2}$ (regime 2). In this case, passing trajectories are replaced by runaway trajectories, so driftons can propagate to much larger $|k_{y}|$ (transport along $k_{y}$ is not suppressed). Then, the corresponding $\bar{W}_{m,n}$ has a much wider distribution; i.e., many harmonics are coupled simultaneously.

Now, let us consider the ZF evolution that would be driven by the above dynamics (which we now consider prescribed for simplicity). As seen from Eq. (29), the ZF is driven by the gradient of the (minus) DW momentum flux,

[TABLE]

which can be numerically calculated from the Wigner function (Ruiz16, ),

[TABLE]

The values of $R(y,t)$ are plotted in the lowest rows in Figs 2 and 3. Since $\kappa_{\pm\lambda}^{2}=p_{{\rm D}}^{2}+(k_{y}\pm\lambda/2)^{2}$, runaways with large $|k_{y}|$ contribute little to $R$, so the global dynamics is largely determined by passing driftons. Particularly, consider the slope of $R$ at the ZF peak ($y=0$). In regime 1, this slope oscillates [Figs. 2(g)-(i)], hence causing oscillations of the ZF amplitude. In contrast, in regime 2, this slope quickly flattens and stays zero indefinitely [Figs. 3(g)-(i)]. This is due to the fact that driftons largely accumulate on runaway and trapped trajectories near the ZF troughs ($y=\pm\pi/q$) and thus cannot influence the ZF peak anymore.

Hence, whether a ZF will oscillate or saturate monotonically depends on whether the “control parameter” $N\doteq u_{{\rm max}}/u_{c,1}$ is smaller or larger than unity. One can also make this estimate more quantitative as follows. As discussed above, the DW spectrum is confined to small $|k_{y}|$ when $N\lesssim 1$; hence, the 4MT can be considered as a reasonable model. By using Eq. (58) for $u_{{\rm max}}$, we obtain

[TABLE]

where we adopted $\delta,\delta^{\prime}\gg 1$ for the GO limit. Then, $N$ can be expressed as

[TABLE]

where $\omega_{{\rm DW}}$ can be recognized as the (absolute value of) the characteristic DW frequency. In summary, the ZF oscillates if $\gamma_{{\rm MI}}\ll\omega_{{\rm DW}}$ and monotonically saturates otherwise.

III.3 Modifications due to full-wave effects

At large enough $q$, the WKE (34) ceases to be valid and one must take into account the deviations from the GO approximation. These deviations are called full-wave effects and can be understood from the coupling between harmonics in Eq. (37). At large $q$, the coupling coefficient $Q$ deviates from unity and becomes inhomogeneous [Eq. (39)]. Recall that

[TABLE]

hence the minimum value of $Q_{a}$ is achieved at $k_{y}+a/2=0$, and

[TABLE]

(Note that $\delta_{\alpha}>1$ is a necessary condition for the MI to happen, and hence $\min Q>0$.) When $\min Q\ll 1$, it can be seen from Eq. (39) that the modes along the diagonals $k_{y}=\pm\lambda/2$ are decoupled from the rest and the initial perturbation $\bar{W}_{0,0}$ propagates mainly along the diagonals, as demonstrated in Fig. 4.

When full-wave effects are important, the critical ZF amplitude can deviate from $u_{c,1}$ given by Eq. (63). Here, we study the critical ZF amplitude by numerically integrating Eq. (38) with stationary $u$ and recording the ZF drive at $y=0$, namely,

[TABLE]

Specifically, we calculate

[TABLE]

where $\overline{\cdot\cdot\cdot}$ is the time-average. (To exclude the initial transient dynamics, only the interval $0.25T<t<T$ is used for the time-averaging, where $T$ is the total integration time.) The critical amplitude is defined as the value of $u$ beyond which $\mathcal{R}$ becomes negligible.

As shown in Figs. 5(c) and (d), the critical amplitude starts to deviate from $u_{c,1}$ as $q$ increases. The following empirical correction can be adopted to account for the finite $q$-dependence of the critical amplitude:

[TABLE]

as seen in Figs. 5(c) and (d). Then, the control parameter (68) becomes

[TABLE]

where $\delta,$ $\delta^{\prime}$, and $\delta_{\alpha}$ are given by Eqs. (26) and (70). The ZF oscillates at $N\ll 1$ and saturates at $N\gtrsim 1$.

IV Numerical simulations

IV.1 Parameter scan

In order to test the above theory of the ZF fate beyond the linear stage, we numerically integrated both the the QL system [Eqs. (28) and (29)] and the NL system [Eqs. (1) and (2)] for various parameters such that

[TABLE]

The second requirement means $\delta_{\alpha}$ shall not be too large, i.e., $q^{2}$ shall not be too small, because then ZF harmonics with wave numbers that are multiples of $q$ have higher growth rates and outpace the fundamental harmonic at the linear stage (see more discussions in Sec. V.1). Overall there are five parameters that determine the system dynamics at the nonlinear stage: $L_{{\rm D}}$, $p$, $\beta$, $q$, and $\Phi_{0}$. For the mHME ($\alpha=0$), we vary only $L_{{\rm D}}$, $\beta$, $q$, and $\Phi_{0}$, because $L_{{\rm D}}$ and $p$ mainly appear as a combination $L_{{\rm D}}^{-2}+p^{2}$, and the remaining $p$ in Eq. (38) only defines the time scale and can be absorbed by a variable transformation $t\to pt$. For the oHME ($\alpha=1$), we vary only $p$, $\beta$, $q$, and $\Phi_{0}$, because it is hard to satisfy the requirement (75) when $L_{{\rm D}}$ is varied with other parameters fixed. For both the oHME and the mHME, the initial conditions are chosen to be

[TABLE]

where $\Phi_{0}$ is real

The simulation results of the QL and NL systems are shown in Figs. 6 (mHME) and 7 (oHME), where we plot the ZF energy

[TABLE]

versus time $t$. As predicted by our theory, in QL simulations, ZFs oscillate at the nonlinear stage if $N\ll 1$ and largely saturate monotonically if $N\gtrsim 1$. For comparison, we also plot the corresponding 4MT solutions (56) for the smallest-$N$ cases in each figure. The fact that the 4MT solutions agree with predictions of the QL theory confirms the 4MT applicability at the nonlinear stage of the MI in the weak-ZF limit.

Note that the ZF energy has a larger maximum at larger $N$. This is explained as follows. At larger $N$, driftons propagate to larger wave numbers and hence end up with lower energy, which is given by (Ruiz16, )

[TABLE]

Due to the total energy conservation, this leaves more energy for the ZFs.

In Figs. 6 and 7, NL simulation results are also shown for the same parameters. The transition from ZF oscillations to saturation is also recovered from these simulations. However, in terms of the oscillation amplitudes and frequencies, good agreement between QL and NL simulations occurs only when the GO approximation is satisfied with reasonable accuracy, say, at

[TABLE]

At smaller $\delta_{\alpha}$, NL simulations show much smaller amplitudes of the ZF oscillations due to the additional DW–DW self-interactions. A brief explanation of the discrepancy between QL and NL simulations is given in Sec. IV.2.

In Fig. 8, we show a snapshot of the DW vorticity $\tilde{w}(x,y)$, the corresponding Wigner function $W(y,k_{y})$, and the ZF profile $U(y)$ from a QL mHME simulation. The snapshot is taken at $\gamma_{{\rm MI}}t=7$ of the $\beta=2.0$ case in Fig. 6(b). This corresponds to the time when the ZF energy reaches the maximum and is about to reverse. The DW vortex structure at $U<0$ is clearly seen and corresponds to trapped driftons in the phase space. In contrast, there is almost no DW activity at $U>0$, because driftons follow passing trajectories and have left this region. At this stage, the ZF velocity is no longer sinusoidal, and have a deep trough and a flat peak. This shape of the ZF is due to the larger ZF drive [Eq. (66)] induced by the trapped trajectories at the ZF trough, and hence cause the ZF to have a larger local amplitude. After this moment of time, passing driftons return to the ZF top since the system is periodic in $y$, and reduce the ZF amplitude; correspondingly, the ZF energy oscillates, as in Fig. 6(b).

IV.2 Difference between the QL and NL models

Here, we discuss why good agreement between QL and NL systems can be achieved when the GO approximation is well satisfied, and why discrepancies arise otherwise. Recall that the MI growth rate $\gamma_{{\rm MI}}$ (25) derived from the 4MT is exact for the QL model but not for the corresponding NL model. Therefore, it is expected that the discrepancy can be attributed, at least partly, to the difference between the QL and NL growth rates. We compare these growth rates by comparing the relative amplitude between the second DW sideband and the first DW sideband at the linear stage, the former being excluded from the 4MT. From Eq. (3), we have

[TABLE]

as a measure of the importance of NL effects. Also, $|\varphi_{\boldsymbol{p}}|=|\Phi_{0}|$ is assumed at the linear stage. From Eq. (5), the coupling coefficients are

[TABLE]

Then

[TABLE]

In the above expression for $\epsilon_{{\rm NL}}$, the first term $|\varphi_{\boldsymbol{p}+\boldsymbol{q}}/\varphi_{\boldsymbol{q}}|$ can be estimated from the 4MT equation (8), which gives

[TABLE]

Since $\partial_{t}\varphi_{\boldsymbol{p}+\boldsymbol{q}}/\partial_{t}\varphi_{\boldsymbol{q}}\sim\varphi_{\boldsymbol{p}+\boldsymbol{q}}/\varphi_{\boldsymbol{q}}$ and $|T(\boldsymbol{q},-\boldsymbol{p},\boldsymbol{p}+\boldsymbol{q})|=pq$, this gives

[TABLE]

and thus

[TABLE]

where $\delta$ and $\delta_{\alpha}$ are given by Eqs. (26) and (70), respectively. When $\delta,\delta_{\alpha}\gg 1$, i.e., $q^{2}\ll p^{2}+(1-\alpha)L_{{\rm D}}^{-2}$, one obtains $\epsilon_{{\rm NL}}\ll 1$. Then, NL effect is expected to be small, and hence QL and NL simulations produce similar results. In contrast, when $\delta_{\alpha}$ approaches unity, $\epsilon_{{\rm NL}}$ can become of order one. Then, the QL model ceases to be an adequate approximation to the NL model.

The above estimate also indicates that, at least for the mHME, when $p^{2}\gtrsim L_{{\rm D}}^{-2}$, one has

[TABLE]

where $\epsilon_{{\rm GO}}\sim\delta^{-1}\ll 1$ under the GO approximation. Therefore, $\epsilon_{{\rm NL}}\ll 1$ at $\epsilon_{{\rm GO}}\ll 1$ and $\epsilon_{{\rm NL}}\gtrsim 1$ when $\epsilon_{{\rm GO}}$ approaches unity (and the MI vanishes when $\epsilon_{{\rm GO}}>1$). Hence, the applicability domain of the QL approximation is roughly the same as that of the GO approximation.

V Comparison with previous studies

V.1 Control parameters: $\boldsymbol{N}$ versus $\boldsymbol{M}$

Let us also compare our results with those in Refs. (Connaughton10, ; Gallagher12, ) where related simulations were performed. Specifically, Ref. (Connaughton10, ) reports NL oHME simulations for $p^{2}\gg q^{2}$ and $L_{{\rm D}}^{-2}=0$. Also, Ref. (Gallagher12, ) reports NL mHME simulations; however, the parameter is chosen such that $p^{2}\gg q^{2},L_{{\rm D}}^{-2}$, and the resulting dynamics is almost identical to that in the oHME. In both cases, the GO assumption is well satisfied ($\delta,\delta^{\prime},\delta_{\alpha}\approx 100\gg 1$). Hence, the difference between NL and QL simulations is expected to be small (see Sec. IV.2), and the results of Refs. (Connaughton10, ; Gallagher12, ) can be compared with ours within the scope of the QL approximation. (We have indeed been able to reproduce all the related results from QL simulations using the same parameters as in Refs. (Connaughton10, ; Gallagher12, ), but we choose not to duplicate the figures here.) Due to the similar choice of the parameters in Ref. (Connaughton10, ) and Ref. (Gallagher12, ), we compare with Ref. (Connaughton10, ) only.

Within the GO limit, our control parameter $N$ [Eq. (74)] is

[TABLE]

A different parameter was proposed in Ref. (Connaughton10, ), namely,

[TABLE]

The authors argue that $M\ll 1/3$ corresponds to ZF oscillations and $M\gtrsim 1/3$ corresponds to monotonic ZF saturation. However, in contrast to our quantitative derivations, only a qualitative argument is provided in Ref. (Connaughton10, ) (also see Sec. V.2). As a result, the parameter $M$ cannot describe the sensitive dependence on $q$ or $L_{{\rm D}}$, as our $N$ does in Figs. 6(a), 6(c), and 7(b). Our $N$ is also better in terms of the predictive power. For example, for the case of $\beta=1.2$ in Fig. 7(c), one has $M=0.47>1/3$, so the ZF is supposed to saturate; however, the ZF actually oscillates, which can be predicted by our $N=0.36<1$.

We note that our theory can also predict the outcomes of three numerical examples used in Ref. (Connaughton10, ), where $M=0.1$, $1$, and $10$, respectively. For $M=0.1$ and $M=10$, our control parameter takes the values $N=0.028$ and $N=3.94$, respectively; hence, our theory predicts that the ZF oscillates in the former case and saturates monotonically in the latter case, which is indeed observed in Ref. (Connaughton10, ). However, Ref. (Connaughton10, ) reports ZF saturation at $M=1$, while our $N=0.39<1$. It may seem then that our theory predicts ZF oscillations instead of saturation, which would be incorrect. Actually, due to the small $q$ in this case, high harmonics of the ZF have much larger linear growth rates and outpace the fundamental harmonic (Fig. 9). This makes our analytic model inapplicable, for it assumes a quasi-sinusoidal ZF. That said, if one calculates $N$ by replacing $q$ with the wave number of the dominant harmonic, then $N$ becomes larger than unity and the ZF saturation is readily anticipated.

V.2 Relevance of the Rayleigh–Kuo threshold for the ZF saturation

The authors of Ref. (Connaughton10, ) proposed a brief explanation of the physical meaning of their parameter $M$, which is as follows. First, they assumed that DWs transfer an order-one fraction of their energy to ZFs, so the ZF maximum amplitude is

[TABLE]

Second, they speculated that the ZF saturates when it reaches the Rayleigh–Kuo (RK) threshold (Kuo49, )

[TABLE]

We believe that this explanation is problematic, namely, for two reasons. First, the estimate (88) contradicts the 4MT estimate (58) that we have confirmed numerically (Figs. 6 and 7). A more accurate estimate within the GO regime would be $u_{{\rm max}}\gtrsim 2q\Phi_{0}$. (Note that this estimate reinstates the dependence on $q$, which is absent in $M$.) Second, the RK criterion does not describe the ZF saturation but rather determines the threshold of the instability of the Kelvin–Helmholtz type that destroys the ZF. (It is also called the “tertiary instability” by some authors (Rogers00, ; Zhu18a, ; Zhu18b, ; Kuo49, ; Zhu18c, ; Numata07, ; Kim2002, ; Rath18, ), and in previous studies we showed that this instability does not exist in the GO limit (Zhu18a, ; Zhu18b, ; Zhu18c, ).) Since the RK threshold corresponds to $u>u_{c,2}$ (Sec. III.2), and $u_{c,2}\gg u_{c,1}$ in the GO regime, we claim that ZFs saturate before the RK threshold is reached. In summary, we believe that our parameter $N$ is more substantiated than the parameter $M$ introduced in Refs. (Connaughton10, ; Gallagher12, ). We also emphasize that our theory withstands the test of numerical simulations, as shown in Sec. IV.1.

VI Conclusions

In this paper we propose a semi-analytic theory that explains the transition between the oscillations and saturation of collisionless ZFs within the QL HME. By analyzing phase-space trajectories of driftons within the GO approximation, we argue that the parameter that controls this transition is $N\sim\gamma_{{\rm MI}}/\omega_{{\rm DW}}$, where $\gamma_{{\rm MI}}$ is the MI growth rate and $\omega_{{\rm DW}}$ is the linear DW frequency. We argue that at $N\ll 1$, ZFs oscillate due to the presence of so-called passing drifton trajectories, and we derive an approximate formula for the ZF amplitude as a function of time in this regime. In doing so, we also extend the applicability of the popular 4MT model, which is commonly used for the linear stage, to nonlinear ZF–DW interactions. We also show that at $N\gtrsim 1$, the passing trajectories vanish and ZFs saturate monotonically, which can be attributed to phase mixing of higher-order sidebands. A modification of $N$ that accounts for effects beyond the GO limit is also proposed. These analytic results are tested against both QL and NL simulations. They also explain the earlier numerical results by Refs. (Connaughton10, ; Gallagher12, ) and offer a revised perspective on what the control parameter is that determines the transition from oscillations to saturation of collisionless ZFs.

Acknowledgements.

This work was supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, and also by the U.S. DOE through Contract No. DE-AC02-09CH11466.