3.1 Primary instability

It is straightforward to show that $\left \{ \varphi , w \right \} = 0$ for Fourier eigenmodes of the form

(3.1)\begin{equation} \varphi=\varphi_{\boldsymbol{k}}\exp\left({\textrm{i}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}-\textrm{i}\varOmega_{\boldsymbol{k}}t}\right)+\textrm{c.c.}, \end{equation}

where $\boldsymbol {k}=(k_{x},k_{y})$. Therefore, a Fourier eigenmode is an exact solution of the system provided that $\varOmega _{\boldsymbol {k}}$ satisfies the following relation:

(3.2a,b)\begin{align} \omega_{\boldsymbol{k}}\doteq\textrm{Re}\,\,\varOmega_{\boldsymbol{k}}=\frac{\beta k_{y}(1+k^{2})}{(1+k^{2})^{2}+\delta_{0}^{2}k_{y}^{2}},\quad \gamma_{\boldsymbol{k}}\doteq\textrm{Im}\,\,\varOmega_{\boldsymbol{k}}=\frac{\beta\delta_{0}k_{y}^{2}}{(1+k^{2})^{2}+\delta_{0}^{2}k_{y}^{2}}-\alpha_{\boldsymbol{k}}D_{\boldsymbol{k}}. \end{align}

Here, $k^{2}\doteq k_{x}^{2}+k_{y}^{2}$, $D_{\boldsymbol {k}}=1+0.01k^{2}$ and we have used (2.7). Also, $\alpha _{\boldsymbol {k}}=1$ for $k_{y}\neq 0$ and $\alpha _{\boldsymbol {k}}=0$ for $k_{y}=0$, and hence a ZF ($k_{y}=0$) corresponds to $\varOmega _{\boldsymbol {k}}=0$, i.e. to a stationary state. From (3.2a,b), it is seen that when $D_{\boldsymbol {k}}=0$, $\gamma _{\boldsymbol {k}}$ is maximized at $(k_{x},k_{y})=(0,1)$. A non-zero $D_{\boldsymbol {k}}$ can modify the value of $\boldsymbol {k}$ that maximizes $\gamma _{\boldsymbol {k}}$, but for the chosen form of $\hat {D}$, (2.8), this modification is very small. Therefore, if one numerically simulates (2.1) with small random noise as the initial conditions, then nonlinear interactions can be neglected at first and coherent DW structures will grow exponentially with typical wavenumber $\boldsymbol {k}\approx (0,1)$, as seen in figure 1.

3.2 Secondary instability

When many Fourier modes are present and have grown to a finite amplitude, the nonlinear term in (2.1) becomes important. This can be seen from the Fourier representation, $\varphi =\sum _{\boldsymbol {k}}\varphi _{\boldsymbol {k}}(t)\exp \left (\textrm {i}\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}\right )$, where (2.1) is written as

(3.3)\begin{equation} \frac{\textrm{d}\varphi_{\boldsymbol{k}}}{\textrm{d}t}=-\textrm{i}\varOmega_{\boldsymbol{k}}\varphi_{\boldsymbol{k}}+\frac{1}{2}\sum_{\boldsymbol{k}_{1},\boldsymbol{k}_{2}}T(\boldsymbol{k},\boldsymbol{k}_{1},\boldsymbol{k}_{2})\delta_{\boldsymbol{k},\boldsymbol{k}_{1}+\boldsymbol{k}_{2}}\varphi_{\boldsymbol{k}_{1}}\varphi_{\boldsymbol{k}_{2}} \end{equation}

and $\delta _{\boldsymbol {k}_1, \boldsymbol {k}_2}$ is the Kronecker symbol. Also,

(3.4)\begin{equation} T(\boldsymbol{k},\boldsymbol{k}_{1},\boldsymbol{k}_{2})\doteq-\frac{\bar{k}_{1}^{2}-\bar{k}_{2}^{2}}{\bar{k}^{2}}(\boldsymbol{k}_{1}\times\boldsymbol{k}_{2})\boldsymbol{\cdot}\hat{\boldsymbol{z}} \end{equation}

are the coefficients that govern the nonlinear mode coupling, $\bar {k}^{2}$ is defined as

(3.5)\begin{equation} \bar{k}^{2}\doteq\alpha_{\boldsymbol{k}}+k^{2}-\textrm{i}\delta_{0}k_{y}, \end{equation}

and similarly for $\bar {k}_{1}^{2}$ and $\bar {k}_{2}^{2}$.

Due to nonlinear interactions, ZFs can be generated from DWs, which process is known as the secondary instability. Here, we use the 4MT model to analyse this instability, namely, by considering a primary DW with $\boldsymbol {k}=(0,k_{y})$, a ZF with $\boldsymbol {q}=(q_{x},0)$, and two DW sidebands with $\boldsymbol {k}_{\pm }=(\pm q_{x},k_{y})$. Assume that the ZF is small, so the exponential growth of the primary DW is unaffected; i.e. $\varphi _{\boldsymbol {k}}=\varphi _{0}\exp \left (-\textrm {i}\varOmega _{\boldsymbol {k}}t\right )$, with $\varphi _{0}$ being a constant. Then, from (3.3), the equations that describe the ZF and the sidebands are as follows (St-Onge Reference St-Onge2017):

(3.6)\begin{gather} \textrm{d}_{t}\varphi_{\boldsymbol{q}}=\frac{k_{y}\textrm{e}^{\gamma_{\boldsymbol{k}} t}}{q_{x}}\left[\left(q_{x}^{2}-\textrm{i}\delta_{+}\right)\varphi_{\boldsymbol{k}_{+}}\varphi_{0}^{*}\textrm{e}^{\textrm{i}\omega_{\boldsymbol{k}}t}-\left(q_{x}^{2}+\textrm{i}\delta_{+}\right)\varphi_{\boldsymbol{k}_{-}}^{*}\varphi_{0}\textrm{e}^{-\textrm{i}\omega_{\boldsymbol{k}}t}\right],\end{gather}

(3.7)\begin{gather} \textrm{d}_{t}\varphi_{\boldsymbol{k}_{+}}=-\textrm{i}\varOmega_{\boldsymbol{k}_{+}}\varphi_{\boldsymbol{k}_{+}}+T(\boldsymbol{k}_{+},\boldsymbol{k},\boldsymbol{q})\varphi_{0}\varphi_{\boldsymbol{q}}\textrm{e}^{-\textrm{i}\varOmega_{\boldsymbol{k}}t}, \end{gather}

(3.8)\begin{gather} \textrm{d}_{t}\varphi_{\boldsymbol{k}_{-}}=-\textrm{i}\varOmega_{\boldsymbol{k}_{-}}\varphi_{\boldsymbol{k}_{-}}+T(\boldsymbol{k}_{-},\boldsymbol{k},-\boldsymbol{q})\varphi_{0}\varphi_{\boldsymbol{q}}^{*}\textrm{e}^{-\textrm{i}\varOmega_{\boldsymbol{k}}t}, \end{gather}

where $\delta _{+}\doteq \delta _{\boldsymbol {k}}+\delta _{\boldsymbol {k}_{+}}=2\delta _{0}k_{y}$. We have also used $\varOmega _{\boldsymbol {k}}=\omega _{\boldsymbol {k}}+\textrm {i}\gamma _{\boldsymbol {k}}$. These equations can be combined to yield a single time-evolution equation for the ZF amplitude $\varphi _{\boldsymbol {q}}$

(3.9)\begin{equation} \frac{\textrm{d}^{3}\varphi_{\boldsymbol{q}}}{\textrm{d}t^{3}}-A\,\frac{\textrm{d}^{2}\varphi_{\boldsymbol{q}}}{\textrm{d}t^{2}}+(B-C)\,\frac{\textrm{d}\varphi_{\boldsymbol{q}}}{\textrm{d}t}-D\varphi_{\boldsymbol{q}}=0. \end{equation}

Here, $A=2\gamma _{+}$, $B=\omega _{-}^{2}+\gamma _{+}^{2}$, $C,D\propto |\varphi {}_{0}\textrm {e}^{\gamma _{\boldsymbol {k}}t}|^{2}$, $\gamma _{+}\doteq \gamma _{\boldsymbol {k}}+\gamma _{\boldsymbol {k}_{+}}$ and $\omega _{-}\doteq \omega _{\boldsymbol {k}}-\omega _{\boldsymbol {k}_{+}}$. The derivation of (3.9) can be found in St-Onge (Reference St-Onge2017). Expressions for $C$ and $D$ can also be found there but will not be important for our discussion; however, note that, compared to St-Onge (Reference St-Onge2017), we have absorbed the coefficient $\textrm {e}^{\gamma _{\boldsymbol {k}}t}$ into the definitions of $C$ and $D$.

When $C$ and $D$ are much larger than $A$ and $B$, $\varphi _{\boldsymbol {q}}$ can grow ‘super-exponentially’ (Rogers et al. Reference Rogers, Dorland and Kotschenreuther2000; St-Onge Reference St-Onge2017), i.e. as an exponential of an exponential. This is also known as the secondary KH instability (Rogers et al. Reference Rogers, Dorland and Kotschenreuther2000). In the opposite case, when $A$ and $B$ dominate over $C$ and $D$, the non-constant solution of (3.9) is approximately

(3.10)\begin{equation} \varphi_{\boldsymbol{q}}\propto\exp\left({(\gamma_{+}\pm\textrm{i}\omega_{-})t}\right). \end{equation}

Since $\gamma _{+}$ decreases as $|q_{x}|$ increases (see (3.2a,b)), the growth rate is maximized at the lowest ZF wavenumber $|q_{x}|=2{\rm \pi} /L_{x}$. In other words, the box-scale ZF grows fastest, with the growth rate given by $\gamma _{+}\approx 2\gamma _{\boldsymbol {k}}$, i.e. twice the growth rate of the primary DW instability.

In the following, we show that exponential growth of the ZF at the box scale is more common than the super-exponential growth, provided that the characteristic amplitude $\varphi _0$ of the initial random noise is small enough. At first, both the primary DW and the sidebands grow exponentially,

(3.11)\begin{equation} |\varphi_{\boldsymbol{k}}|\sim|\varphi_{\boldsymbol{k}_{\pm}}|\sim|\varphi_{0}|\textrm{e}^{\gamma_{\boldsymbol{k}}t}, \end{equation}

while the ZF amplitude remains at the noise level. Then, DWs grow for some time $t_{p}$ before they begin to affect ZFs. Assume that at $t=t_{p}$, the box-scale ZF with the amplitude $\varphi _{\boldsymbol {q}}\sim \varphi _{0}$ starts to grow with the growth rate $\gamma _{+}\approx 2\gamma _{\boldsymbol {k}}$; then, $\delta _{+}=2\delta _{0}k_{y}\gg q_{x}^{2}$, and we have from (3.6) that

(3.12)\begin{equation} |\partial_{t}\varphi_{\boldsymbol{q}}|\sim 2\gamma_{\boldsymbol{k}}|\varphi_{{0}}|\approx\frac{2|\varphi_{\boldsymbol{k}}||\varphi_{\boldsymbol{k}_{+}}|k_{y}\delta_{+}}{q_{x}}\approx\frac{2k_{y}\delta_{+}|\varphi_{{0}}\textrm{e}^{\gamma_{\boldsymbol{k}}t_{p}}|^{2}}{q_{x}}.\end{equation}

This leads to

(3.13)\begin{equation} C,D\propto\frac{q_{x}\gamma_{\boldsymbol{k}}|\varphi_{0}|}{2\delta_{0}k_{y}^{2}}. \end{equation}

Therefore, $C$ and $D$ are small when the initial noise level $|\varphi _{0}|$ is small enough; hence, the assumptions made above are self-consistent, namely, $A$ and $B$ are indeed much larger than $C$ and $D$, and the box-scale ZF with wavenumber $q_{x}=2{\rm \pi} /L_{x}$ grows fastest with the growth rate $2\gamma _{\boldsymbol {k}}$.

The secondary instability will persist for some time $t_{s}$ until ZFs grows up to a finite amplitude that is enough to significantly distort the DW structure. Using the result from Zhu, Zhou & Dodin (Reference Zhu, Zhou and Dodin2018b), this amplitude can be estimated as follows (also see (B 11)):

(3.14)\begin{equation} \varphi_{c}=\frac{\beta/q_{x}}{2(1+k_{y}^{2})-q_{x}^{2}}. \end{equation}

At $\varphi _{\boldsymbol {q}}\ll \varphi _{c}$, DWs do not ‘see’ the ZF and hence keep growing exponentially, while at $\varphi _{\boldsymbol {q}}\gtrsim \varphi _{c}$ the system enters the fully nonlinear regime. Therefore, $t_{s}$ is the time when the ZF amplitude grows from $\varphi _{0}$ to $\varphi _{c}$, and it can be estimated as follows:

(3.15)\begin{equation} t_{s}=\frac{1}{2\gamma_{\boldsymbol{k}}}\ln\frac{\varphi_{c}}{\varphi_{0}}. \end{equation}

Note that (3.14) is obtained from the modified Hasegawa–Mima system, so it is based on the assumption that $\delta _0 = 0$. For non-zero $\delta _0$, it is modified accordingly (see (B 11)), but the above estimate is sufficient for our qualitative description.

By the time when the system enters the fully nonlinear regime, the DW amplitude becomes $|\varphi _{\boldsymbol {k}}|\sim \varphi _0\exp {\gamma _{\boldsymbol {k}}(t_{s}+t_{p})}$, which can be estimated from (3.12) and (3.15) as

(3.16)\begin{equation} |\varphi_{\boldsymbol{k}}|\sim\sqrt{\frac{q_{x}\gamma_{\boldsymbol{k}}\varphi_{c}}{2\delta_{0}k_{y}^{2}}}. \end{equation}

From (2.6a,b), the corresponding DW and ZF energies are as follows:

(3.17a,b)\begin{equation} E_{\textrm{ZF}}\sim\frac{\beta^{2}}{8(1+k_{y}^{2})^{2}},\quad E_{\textrm{DW}}\sim\frac{\beta\gamma_{\boldsymbol{k}}}{8\delta_{0}k_{y}^{2}}, \end{equation}

where we assumed $q_x^2\ll 1+k_y^2$. Using (3.2a,b) for $\gamma _{\boldsymbol {k}}$ and assuming $D_{\boldsymbol {k}}=0$ for simplicity, we obtain

(3.18)\begin{equation} \frac{E_{\textrm{ZF}}}{E_{\textrm{DW}}}\sim 1 + \frac{\delta_0^2k_y^2}{(1+k_y^2)^2}. \end{equation}

This shows that the ZF energy and the DW energy are roughly equal when the system enters the fully nonlinear regime, since $\delta _0$ and $k_y$ are of order unity. This conclusion will be used to estimate the ZF curvature in § 5.

These predictions are in agreement with numerical simulations (figure 2). This indicates that the 4MT captures the basic dynamics of the primary and the secondary instabilities. However, as shown below, the 4MT does not capture essential features of the TI, and thus more accurate models are needed to describe the TI and the Dimits shift.

3.3 Tertiary instability

In the fully nonlinear regime, DW turbulence becomes inhomogeneous and localized at the extrema of the ZF velocity $U$ (figure 1). To understand the DW dynamics in this case, let us linearize (2.1) to obtain

(3.19)\begin{equation} \partial_{t}\tilde{w}+U\partial_{y}\tilde{w}-(\beta+U'')\partial_{y}\tilde{\varphi}+\hat{D}\tilde{w}=0, \end{equation}

where

(3.20a,b)\begin{equation} \tilde{w}=(\nabla^{2}-1+\textrm{i}\hat{\delta})\tilde{\varphi},\quad U''\doteq\textrm{d}^{2}U(x)/\textrm{d}x^{2}. \end{equation}

For given boundary conditions in $x$, eigenmodes of (3.19) can be searched for in the form

(3.21a,b)\begin{equation} \tilde{w}=w(x)\exp\left({\textrm{i}(k_{y}y-\omega t)}\right),\quad\tilde{\varphi}=\left(\frac{\textrm{d}^{2}}{\textrm{d}x^{2}}-k_{y}^{2}-1+\textrm{i}\delta_{0}k_{y}\right)^{-1}\tilde{w}, \end{equation}

which leads to the following equation for $w(x)$:

(3.22a,b)\begin{equation} \omega w=\hat{H}w,\quad \hat{H}(\hat{x},\hat{k}_x)\doteq k_{y}\hat{U}+k_{y}(\beta+\hat{U}'')\hat{\bar{k}}^{-2}-\textrm{i}\hat{D}, \end{equation}

where

(3.23a–c)\begin{equation} \hat{U}=U(\hat{x}),\quad \hat{k}_x=-\textrm{i}\,\textrm{d}/\textrm{d}x,\quad\hat{\bar{k}}^{2}=1+k_{y}^{2}+\hat{k}_x^2-\textrm{i}\delta_{0}k_{y}. \end{equation}

If an eigenvalue $\omega$ exists and

(3.24)\begin{equation} \gamma_{\textrm{TI}}\doteq\textrm{Im}\,\,\omega>0, \end{equation}

then the perturbation grows exponentially. This is the TI.

Equation (3.19) does not have an analytic solution for an arbitrary profile $U$, but a general understanding can be developed by considering special cases. In Zhu et al. (Reference Zhu, Zhou, Ruiz and Dodin2018c), we considered the ZF velocity profile

(3.25)\begin{equation} U(x)=u\cos q_{x}x, \end{equation}

with $\hat {\delta } = \hat {D} = 0$. In this case, the system exhibits an instability of the KH type provided that $q_{x}^{2}>1$ and $q_x^2u>\beta$. In Zhu et al. (Reference Zhu, Zhou, Ruiz and Dodin2018c), we also discussed a generalization to periodic non-sinusoidal profiles. However, generalizing those results to non-zero $\hat {\delta }$ and $\hat {D}$ is challenging. The common approach is to adopt the 4MT again, i.e. to assume a DW perturbation with $\boldsymbol {k}=(0,k_{y})$ and two sidebands with $\boldsymbol {k}_{\pm }=({\pm }q_{x},k_{y})$ as small perturbations (Kim & Diamond Reference Kim and Diamond2002; St-Onge Reference St-Onge2017; Rath et al. Reference Rath, Peeters, Buchholz, Grosshauser, Seiferling and Weikl2018; Zhu, Zhou & Dodin Reference Zhu, Zhou and Dodin2018a). In particular, St-Onge (Reference St-Onge2017) derived $\gamma _{\textrm {TI}}$ within the 4MT and estimated the Dimits shift by finding a sufficient condition for $\gamma _{\textrm {TI}}=0$. However, the 4MT-based approach is not entirely satisfactory, because the ZF is typically far from sinusoidal, as seen in simulations. Even more importantly, the 4MT approach ignores the fact that there are multiple TI modes with different growth rates. As we show below, understanding the variety of these modes is essential for understanding the Dimits shift.

Let us assume the same sinusoidal ZF profile (3.25) as in St-Onge (Reference St-Onge2017) for now, and let us calculate the corresponding eigenmodes (3.22a,b) numerically, assuming periodic boundary conditions $x$. In this case, we can search for solutions in the form

(3.26)\begin{equation} w(x)=\sum_{n=-N}^{N}w_{n}\textrm{e}^{\textrm{i}nq_{x}x}, \end{equation}

where $N$ is some large enough integer. In other words, we truncate the Fourier series by keeping only the first $2N+1$ Fourier modes. This turns (3.22a,b) into a vector equation for $\{w_{-N}, \ldots w_0, \ldots w_N\}$, where $\hat {H}$ becomes a $(2N+1)\times (2N+1)$ matrix. Then, one finds $2N+1$ eigenmodes with complex eigenfrequencies. Typical numerical eigenmodes are illustrated in figure 3. It is seen that the TI-mode structure is localized at the maximum ($x=0$) or minimum ($x=-{\rm \pi} /q_{x}$) of the ZF velocity and has either even or odd parity because of the symmetry of $U$. Within the figure, the eigenmodes localized at the ZF minimum can be labelled by the integer $m=0,1,2,\ldots$, which also indicates the parity of $w(x)$. Eigenmodes localized near the ZF maximum can be labelled similarly. Note that in order for a mode to be localized, the ZF must be large-scale, namely, $q_x^2\ll 1+k_y^2$, which is consistent with numerical simulations.

Figure 3. The first four tertiary eigenmodes found numerically using the ZF velocity profile (3.25). The ordering is such that $\gamma _{\textrm {TI}}$ decreases from left to right. The first two eigenmodes are runaway and trapped modes, respectively. The parameters are $\beta =6$, $\delta _{0}=1.5$, $q_{x}=0.4$ and $u=10$. The first row shows the eigenmode structures $\tilde {w}(x,y)=\textrm {Re}\,[w(x)\textrm {e}^{\textrm {i}k_{y}y}]$ ((3.21a,b), colour), the ZF velocity $U$ (green curve) and the analytic mode structure $\tilde {w}=\textrm {Re}\,[\mathsf {H}_m(x)\exp (S+\textrm {i}k_y y)]$ ((4.6a,b), dashed contour), where $m=0$ for (a) and (b), $m=1$ for (c) and $m=2$ for (d). The second row shows the corresponding Wigner function $W(x,k_x)$ ((A 5), colour) and the isosurfaces of the drifton Hamiltonian $\mathcal {H}$ ((B 3), dashed contour). The striped structure of $W$ away from the actual location of DW quanta is a signature of a quantum-like ‘cat state’ (Weinbub & Ferry Reference Weinbub and Ferry2018).

Apart from the eigenmode structures, we also show in figure 3 their corresponding Wigner functions $W(x,k_x)$ (A 5) and contour plots of the drifton Hamiltonian $\mathcal {H}$ (B 3). The Wigner function can be understood as the distribution function of ‘driftons’ (DW quanta) in the $(x,k_{x})$ phase space (Smolyakov & Diamond Reference Smolyakov and Diamond1999; Ruiz et al. Reference Ruiz, Parker, Shi and Dodin2016; Zhu et al. Reference Zhu, Zhou, Ruiz and Dodin2018c), and its shape is expected to align with the contours of $\mathcal {H}$ . Then, eigenmodes are naturally centred at phase-space equilibria of $\mathcal {H}$, namely,

(3.27)\begin{equation} \partial_{x}\mathcal{H}=\partial_{k_{x}}\mathcal{H}=0\quad\Rightarrow\quad U' = k_x=0.\end{equation}

This explains eigenmode localization near extrema of $U$. [Strictly speaking, (3.27) stems from our approximation of sinusoidal flow (3.25), which ensures that $U'$ and $U'''$ become zero at same locations. Nevertheless, (3.27) remains a good approximation as long as ZFs are large-scale, i.e. $|U'''/\bar {k}^2|\ll U'$.] Maxima of $U$ (even $n$) correspond to phase-space islands encircled by ‘trapped’ trajectories, and minima of $U$ (odd $n$) correspond to saddle points passed by the ‘runaway’ trajectories (Zhu et al. Reference Zhu, Zhou and Dodin2018a,Reference Zhu, Zhou and Dodinb,Reference Zhu, Zhou, Ruiz and Dodinc). Hence, we call the modes localized near maxima and minima of $U$ trapped and runaway modes, respectively. (See appendix B for more discussions on drifton phase-space trajectories.) In § 4, we provide analytic calculation of the TI growth rates based on the above observations.

4.1 Analogy with a quantum harmonic oscillator

As seen in figure 3, tertiary modes are centred at the phase-space equilibria. Based on this, let us expand the Hamiltonian up to the second order both in $x$ and in $\hat {k}_x$. Specifically, we approximate the ZF velocity with a parabola

(4.1)\begin{equation} U\approx U_0+\frac{1}{2}\,\mathcal{C}x^2, \end{equation}

where $U_0$ is the local ZF velocity and $\mathcal {C}\doteq U''(0)$ is the local ZF curvature. For the sinusoidal velocity (3.25), this corresponds to $U_0={\pm }u$ and $\mathcal {C}={\mp }q_x^2 u$. We also make the approximation that $\hat {D}\approx D_0\doteq D_{\boldsymbol {k}=(0,k_y)}$ and

(4.2a,b)\begin{equation} \hat{k}^{-2}\approx k_0^{-2}+k_0^{-4}\frac{\textrm{d}^2}{\textrm{d}x^2},\quad k_0^2\doteq 1+k_y^2-\textrm{i}\delta_0 k_y. \end{equation}

Then, the Hamiltonian operator $\hat {H}$ (3.22a,b) is approximated as

(4.3)\begin{equation} \hat{H}\approx k_y U_0 +\frac{1}{2}k_y\mathcal{C}\hat{x}^2+k_y(\beta+\mathcal{C})\left(k_0^{-2}+k_0^{-4}\frac{\textrm{d}^2}{\textrm{d}x^2}\right)-\textrm{i}D_0, \end{equation}

and the corresponding eigenmode equation (3.22a,b) becomes

(4.4)\begin{equation} \left(-\tau^2\frac{\textrm{d}^2}{\textrm{d}x^2}+x^2\right)w=\varepsilon w. \end{equation}

It is the same equation that describes a quantum harmonic oscillator, except that here the coefficients are complex; specifically,

(4.5a,b)\begin{equation} \tau^2\doteq-\frac{2}{k_0^4}\left(1+\frac{\beta}{\mathcal{C}}\right),\quad\varepsilon\doteq\frac{2[\omega-k_y U_0+\textrm{i}D_0-k_y(\beta+\mathcal{C})/k_0^2]}{k_y\mathcal{C}}. \end{equation}

Note that the coefficients are different at minima and maxima of $U$, as they depend on the sign of $\mathcal {C}$. Also note that for runaway modes, we have shifted the coordinate as $x\to x+{\rm \pi} /q_{x}$ to recentre the ZF minimum at $x=0$.

Following the standard procedure known from quantum mechanics (Sakurai Reference Sakurai and Tuan1994), one can show that the asymptotic behaviour of the solution at large $|x|$ is

(4.6a,b)\begin{equation} w(x)\sim \textrm{e}^{S(x)},\quad S(x)=-\frac{x^2}{2\tau}=-\left[\frac{\textrm{i}(1+k_y^2)+\delta_0k_y}{2\sqrt{2(1+\beta/\mathcal{C})}}\right]x^2. \end{equation}

To ensure that $w\to 0$ at large $|x|$, we require $\textrm {Im}\sqrt {1+\beta /\mathcal {C}}>0$ if $1+\beta /\mathcal {C}<0$. We also assumed that $\delta _0,k_y>0$. Then, letting $w=\phi (x)\exp S(x)$, we obtain

(4.7)\begin{equation} \phi''-\frac{2x}{\tau}\phi'+\frac{\varepsilon-\tau}{\tau^2}\phi=0. \end{equation}

Solutions are $\phi =\mathsf {H}_m(x/\sqrt {\tau })$, where $\mathsf {H}_m$ are Hermite polynomials, $m=0,1,2,\ldots$, and

(4.8)\begin{equation} \varepsilon=(2m+1)\tau. \end{equation}

Therefore, for each sign of $\mathcal {C}$, eigenmodes are labelled by $m$. In figure 3, these approximate solutions are compared with numerical solutions of (3.22a,b). In the following, we shall focus on the two modes with $m = 0$, since they are the most unstable. In this case, $\phi =\mathsf {H}_0$ is constant and $\varepsilon =\tau$. This corresponds to $\tilde {w}=\textrm {Re}\,[\exp (S+\textrm {i}k_y y)]$, and the eigenfrequencies are found from (4.5a,b) to be

(4.9a,b)\begin{equation} \omega=\bar{\varOmega}+k_yU_0-\frac{\textrm{i}k_y\mathcal{C}\sqrt{(1+\beta/\mathcal{C})/2}}{1+k_y^2-\textrm{i}\delta_0 k_y},\quad \bar{\varOmega}=\frac{k_y(\beta+\mathcal{C})}{k_0^2}-\textrm{i}D_0. \end{equation}

Here, $\bar {\varOmega }$ is the primary-mode eigenfrequency $\varOmega _{\boldsymbol {k}}$ (3.2a,b) modified by $\mathcal {C}$, $k_yU_0$ is the local Doppler shift, and the remaining term in $\omega$ vanishes at zero $\mathcal {C}$. Note that at $\mathcal {C}=0$, $\omega$ reduces to the primary-mode frequency $\varOmega _{\boldsymbol {k}}$ at $\boldsymbol {k}=(0,k_y)$. Hence, TI modes found here can be interpreted as standing primary modes modified by ZFs. Accordingly, the TI growth rate $\gamma _{\textrm {TI}}$ approaches the primary-instability growth rate in the limit $\mathcal {C} \to 0$.

Let us examine the validity of our approximation in (4.3). First, the parabolic approximation of $U$ is valid if the mode spatial width in $x$, which is determined by $\mathcal {C}$, is much smaller than $q_x^{-1}$, which is the characteristic scale of ZFs. Specifically, for the sinusoidal ZF (3.25), we have $\mathcal {C}=q_x^2 u$, so the parabolic approximation (4.1) is valid at small enough $q_x$ and large enough $u=\mathcal {C}/q_x^2$. Second, the expansion of $\hat {k}^{-2}$ in (4.2a,b) is valid at $|k_x^2|\ll |k_0^2|$, where $k_x$ is the characteristic mode wavenumber in $x$. From (4.4), $k_x$ can be estimated as $k_x=1/\sqrt {\tau }$. Then, the requirement $|k_x^2|\ll |k_0^2|$ leads to $|\sqrt {2(1+\beta /\mathcal {C})}|\gg 1$, or equivalently, $|\mathcal {C}|\ll \beta$, and one expects that the approximation in (4.2a,b) becomes invalid as $|\mathcal {C}|$ approaches $\beta$. Therefore, in the following, we restrict our consideration to the parameter regime $|\mathcal {C}|<\beta$, which is also the regime relevant to our numerical simulations.

The TI growth rate $\gamma _{\textrm {TI}}$ is obtained by taking the imaginary part of $\omega$. Within the regime $|\mathcal {C}|<\beta$, let us introduce the notation

(4.10)\begin{equation} \bar{\gamma}\doteq\textrm{Im}\,\bar{\varOmega}=\frac{\delta_0k_y^2(\beta+\mathcal{C})}{(1+k_y^2)^2+\delta_0^2k_y^2}-D_0, \end{equation}

which is the primary-instability growth rate $\gamma _{\boldsymbol {k}}$ (3.2a,b) modified by $\mathcal {C}$. Then, for the runaway mode (labelled with superscript ‘$R$’), which corresponds to $\mathcal {C}>0$, one has

(4.11)\begin{equation} \gamma_{\textrm{TI}}^{R}=\bar{\gamma}-\sqrt{\frac{1}{2}\left(1+\frac{\beta}{\mathcal{C}}\right)}\,\frac{(1+k_y^2)k_y\mathcal{C}}{(1+k_y^2)^2+\delta_0^2k_y^2}. \end{equation}

For the trapped mode (labelled with superscript ‘$T$’), which corresponds to $\mathcal {C}<0$, one has

(4.12)\begin{equation} \gamma_{\textrm{TI}}^{T}=\bar{\gamma}+\sqrt{\frac{1}{2}\left(\left|\frac{\beta}{\mathcal{C}}\right|-1\right)}\,\frac{\delta_0k_y^2\mathcal{C}}{(1+k_y^2)^2+\delta_0^2k_y^2}.\end{equation}

Figure 5(c) shows that these formulas are in good agreement with our numerical calculations of the eigenvalues. Also, we have verified (not shown) that the results in figure 5(c) are insensitive to $q_x$ as long as $q_x$ is small, more specifically, $q_x^2\ll 1+k_y^2$.

Figure 4. The structure of the tertiary modes with $m=0$ in the zonal-velocity profile (3.25). The parameters are $\delta _{0}=1.6$, $\beta =5$, $q_{x}=0.2$, $u=50$ (hence, $\mathcal {C}={\pm }2$) and $\hat {D}$ given by (2.8). These parameters result in $\gamma _{\textrm {TI}}^{R}=-0.276$ and $\gamma _{\textrm {TI}}^{T}=-0.587$. $(a)$ The Wigner function $W(x,k_x)$ of the runaway mode (colour), the local $U$ (magenta curve) and the runaway trajectory (dashed curve; see (4.21)). Note that we have shifted the coordinates as $x\to x+{\rm \pi} /q_x$ to recentre the ZF minimum at $x=0$. $(b)$ The structure of each term in (4.17) calculated from $W$ of the runaway mode in $(a)$. $(c)$ The Wigner function $W(x,k_x)$ of the trapped mode (colour), the local $U$ (magenta curve) and isosurfaces of $\mathcal {H}$ (dashed contours; see (4.14a,b)). In this figure, $\Delta x$ and $\Delta k_x$ denote the characteristic widths of the mode in the $x$ and $k_x$ directions, respectively. $(d)$ Same as $(b)$ but for the trapped mode.

Figure 5. $(a)$ The empirical factor $\eta _{R}$ (4.24) as a function for $\delta _0$: numerical values of a sinusoidal ZF (3.25) (blue circles) versus the fitting formula (black curve). The parameters are $\beta =6$, $q_x=0.4$, $u=10$ and $\hat {D}$ given by (2.8). It is found that $\eta _{R}$ is not sensitive to $u$. $(b)$ Same as $(a)$ except for $\eta _{T}$ (4.28a,b). It is found that $\eta _{T}$ is not sensitive to $u$ at $u< \beta /q_x^2$. $(c)$ The TI growth rates versus $|\mathcal {C}|=q_x^{2}u$ at $\delta _{0}=1.5$, $\beta =6$, $q_{x}=0.4$ and varying $u$. Black curves: numerical solutions of (3.22a,b) indicated by the superscript ‘$N$’. Multiple branches are shown, with the two most unstable branches being the runaway mode and the trapped mode. Blue dashed curve and red dash-dotted curve: analytic formulas (4.11) and (4.12). The superscript ‘1’ corresponds to predictions made using the approach described in § 4.1. Blue circles and red squares: analytic formulas (4.20) with $\eta _{{R}}=0.595$ and (4.27) with $\eta _{{T}}=1.2$. The superscript ‘2’ corresponds to predictions made using the approach described in § 4.2.

Notably, while the trapped-mode growth rate always decreases with $|\mathcal {C}|$, the runaway-mode growth rate can increase at large $\mathcal {C}$ if $\delta _0$ is large. In fact, at $\mathcal {C}\gg \beta$, (4.11) becomes

(4.13)\begin{equation} \gamma_{\textrm{TI}}^{R}\approx \left[\frac{\delta_0k_y^2}{(1+k_y^2)^2+\delta_0^2k_y^2}-\sqrt{\frac{1}{2}}\,\frac{(1+k_y^2)k_y}{(1+k_y^2)^2+\delta_0^2k_y^2}\right]\mathcal{C}, \end{equation}

which predicts that $\gamma _{\textrm {TI}}^{R}$ increases with $\mathcal {C}$ if $\delta _0>(k_y + k_y^{-1})/\sqrt {2}$. Therefore, it is possible that the TI can develop in strong ZFs, but the physical mechanism is very different from the KH mode, as will be discussed in § 4.3. (Strictly speaking, (4.11) becomes invalid at $\mathcal {C}\gg \beta$. Nevertheless, we have verified from numerical calculations (not shown) that at $k_y=1$, $\gamma _{\textrm {TI}}^{R}$ indeed increases with $\mathcal {C}$ at large $\mathcal {C}$, if $\delta _0\gtrsim 1.7$.)

4.2 Alternative approach

An alternative formula for $\gamma _{\textrm {TI}}$ can be obtained using the WME for the Wigner function $W$ of the fluctuations $\tilde {w}$ (appendix A). This approach is somewhat more accurate because the Hamiltonian is expanded only in $x$ but not in $k_x$. As in § 4.1, let us assume $U=U_0+\mathcal {C}x^2/2$. Then, $U''=\mathcal {C}$ is constant, $U'''$ vanishes and the drifton Hamiltonian is simplified down to (appendix A)

(4.14a,b)\begin{equation} \mathcal{H}=k_{y}U_{0}+\frac{1}{2}k_{y}\mathcal{C}x^{2}+k_{y}(\beta+\mathcal{C})\textrm{Re}\,\left(\frac{1}{\bar{k}^{2}}\right),\quad\varGamma=k_{y}(\beta+\mathcal{C})\textrm{Im}\left(\frac{1}{\bar{k}^{2}}\right)-D_{\boldsymbol{k}}, \end{equation}

where $\bar {k}^2=1+k_x^2+k_y^2-\textrm {i}\delta _0k_y$. Then, the WME (A 4) acquires the form

(4.15)\begin{equation} \frac{\partial W}{\partial t} = k_{y}\mathcal{C}x\,\frac{\partial W}{\partial k_{x}}-V_{{g}}\,\frac{\partial W}{\partial x}+2\varGamma W+\frac{\partial Q}{\partial x},\end{equation}

where

(4.16)\begin{equation} V_{{g}}(k_{x})\doteq\frac{\partial\mathcal{H}}{\partial k_{x}}=k_{y}(\beta+\mathcal{C})\frac{\partial\textrm{Re}\,(1/\bar{k}^{2})}{\partial k_{x}} \end{equation}

is the drifton group velocity. (Details of the drifton dynamics are discussed in appendix B.) The value of $Q$ is given by (A 11), but it is not important for our calculations, because we are interested only in the spatial integral of (4.15). Since $V_{{g}}$ and $\varGamma$ are independent of $x$, integrating (4.15) over $x$ leads to

(4.17)\begin{equation} 2\gamma_{\textrm{TI}}f_{1}=k_{y}\mathcal{C}\,\frac{\partial f_{2}}{\partial k_{x}}+2\varGamma f_{1}, \end{equation}

where we have replaced $\partial _{t}$ with $2\gamma _{\textrm {TI}}$ and introduced

(4.18a,b)\begin{equation} f_{1}(k_{x})\doteq\int W\,\textrm{d}x,\quad f_{2}(k_{x})\doteq\int xW\,\textrm{d}x. \end{equation}

The functions are shown in figure 4 for the runaway mode and for the trapped mode, respectively. Note that from comparing (4.17) with (4.9a,b), it is seen that $f_1$ is associated with the modified frequency $\bar {\varOmega }$, namely, $\varGamma =\bar {\gamma }=\textrm {Im}\,\bar {\varOmega }$; meanwhile, $\partial f_2/\partial k_x$ is associated with the additional term in (4.9a,b) that vanishes at $\mathcal {C}=0$.

To obtain $\gamma _{\textrm {TI}}$ from (4.17), one needs to find the relation between $f_1$ and $f_2$. Let us first consider the runaway mode. As shown in figure 4(a), the Wigner function of this mode peaks along $x=x_{{R}}(k_{x})$, which is the runaway trajectory that passes through the saddle point of $\mathcal {H}$ at $x=k_x=0$, and is given by (4.21) below. Therefore, let us adopt $f_{2}\approx x_{{R}}f_{1}$; then,

(4.19)\begin{equation} \frac{\partial f_{2}}{\partial k_{x}}\approx\frac{\partial x_{{R}}}{\partial k_{x}}f_{1}+x_{{R}}\frac{\partial f_{1}}{\partial k_{x}}.\end{equation}

With this assumption, let us evaluate (4.17) at $k_{x}=0$, where $\partial f_{1}/\partial k_{x}=0$ because $f_1$ is even in $k_x$ due to the symmetry of $\hat {H}$ (see figure 4); then, we find

(4.20)\begin{equation} \gamma_{\textrm{TI}}^{R}= \left. \left(\varGamma+\frac{k_{y}\mathcal{C}}{2\eta_{{R}}}\frac{\partial x_{{R}}}{\partial k_{x}}\right)\right|_{k_{x}=0}. \end{equation}

Here, the first term $\varGamma$ is given by (4.14a,b). The second term is negative because ${\partial x_{{R}}/\partial k_{x}<0}$ (see (4.23) below). The coefficient $\eta _{{R}}>0$ is an empirical factor that compensates for the inaccuracy of (4.19). We proceed to determine $x_{{R}}(k_{x})$ and $\eta _{{R}}$. The runaway trajectory $x_{{R}}$ is determined from (4.14a,b) by equating $\mathcal {H}$ to its value at the origin $(x,k_{x})=(0,0)$ and solving $x$ as a function of $k_{x}$. This gives

(4.21)\begin{equation} x_{{R}}(k_{x})=\pm\sqrt{2\left(1+\frac{\beta}{\mathcal{C}}\right)} \sqrt{\frac{1+k_{y}^{2}}{(1+k_{y}^{2})^{2}+\delta_{0}^{2}k_{y}^{2}}-\frac{1+k_{y}^{2}+k_{x}^{2}}{(1+k_{y}^{2}+k_{x}^{2})^{2}+\delta_{0}^{2}k_{y}^{2}}}, \end{equation}

where the plus sign is for $k_x<0$ and the minus sign is for $k_x>0$. Figure 4(a) demonstrates that this solution indeed correlates well with the actual runaway-mode structure. Also note that $x_{{R}}$ is finite, namely,

(4.22)\begin{equation} x_{{R}}(k_{x}=\infty)=-\sqrt{2\left(1+\frac{\beta}{\mathcal{C}}\right)\frac{1+k_{y}^{2}}{(1+k_{y}^{2})^{2}+\delta_{0}^{2}k_{y}^{2}}}. \end{equation}

From (4.21), we obtain

(4.23)\begin{equation} \left. \frac{\partial x_{{R}}}{\partial k_{x}}\right|_{k_{x}=0}=-\frac{\sqrt{2(1+\beta/\mathcal{C})\left[(1+k_{y}^{2})^{2}-\delta_{0}^{2}k_{y}^{2}\right]}}{(1+k_{y}^{2})^{2}+\delta_{0}^{2}k_{y}^{2}}. \end{equation}

Notably, $\partial x_{{R}}/\partial k_{x}$ becomes zero at $\delta _{0}=|k_{y}+k_{y}^{-1}|$, which corresponds to the transition from runaway to trapped trajectory at the ZF minimum, as shown in figure 8.

Now, let us consider the correction factor $\eta _{R}$, which can be formally defined as

(4.24)\begin{equation} \eta_{{R}}\doteq\left. \left(\frac{\partial x_{{R}}}{\partial k_{x}}\,\frac{f_{1}}{\partial f_{2}/\partial k_{x}}\right)\right|_{k_{x}=0}. \end{equation}

We determine $\eta _{{R}}$ numerically from the eigenmode structures obtained in § 3.3. It can be shown that if $\hat {D}=0$, then rescaling $t\to {k_y\mathcal {C}t}/{(1+k_y^2)}$, $x\to x\sqrt {1+k_y^2}$ and $k_x\to {k_x}/{\sqrt {1+k_y^2}}$ leaves only two parameters in the WME (4.15), namely, $\mathcal {C}/\beta$ and $\delta _0k_y/(1+k_y^2)$; hence, $\eta _{{R}}$ mainly depends on these two parameters. Numerically, we see that $\eta _{{R}}$ changes little as $\mathcal {C}/\beta$ varies from zero to unity. Meanwhile, the dependence of $\eta _{{R}}$ on $\delta _{0}k_{y}/(1+k_{y}^{2})$ is shown in figure 5(a), which suggests the following approximation:

(4.25)\begin{equation} \eta_{{R}}\approx\sqrt{1-\left(\frac{\delta_{0}k_{y}}{1+k_{y}^{2}}\right)^{2}}. \end{equation}

Then, (4.20) is simplified as

(4.26)\begin{equation} \gamma_{\textrm{TI}}^{R}\approx\varGamma|_{k_{x}=0}-\sqrt{\frac{1}{2}\left(1+\frac{\beta}{\mathcal{C}}\right)}\frac{(1+k_{y}^{2})k_{y}\mathcal{C}}{(1+k_{y}^{2})^{2}+\delta_{0}^{2}k_{y}^{2}}. \end{equation}

Remarkably, this formula is identical to (4.11) that was obtained in § 4.1 by drawing the analogy with a quantum harmonic oscillator.

The above approach can also be applied to the trapped mode. Similarly to (4.20), the trapped-mode growth rate can be expressed as follows:

(4.27)\begin{equation} \gamma_{\textrm{TI}}^{T}=\left. \left(\varGamma+\frac{k_{y}\mathcal{C}}{2\eta_{T}}\,\frac{\Delta x}{\Delta k_{x}}\right)\right|_{k_{x}=0}, \end{equation}

where

(4.28a,b)\begin{equation} \frac{\Delta x}{\Delta k_{x}}\doteq\frac{\sqrt{2(|\beta/\mathcal{C}|-1)\left[(1+k_{y}^{2})^{2}-\delta_{0}^{2}k_{y}^{2}\right]}}{(1+k_{y}^{2})^{2}+\delta_{0}^{2}k_{y}^{2}},\quad\eta_{{T}}\doteq\left.\left(\frac{\Delta x}{\Delta k_{x}}\,\frac{f_{1}}{\partial f_{2}/\partial k_{x}}\right)\right|_{k_{x}=0}. \end{equation}

Here, $\mathcal {C}<0$, and we consider the regime $\beta /\mathcal {C}<-1$. Also, $\Delta x/\Delta k_{x}$ is not the slope of the runaway trajectory but the ratio of the $x$-axis radius and the $k_{x}$-axis radius of the elliptic trapped trajectory near $(x,k_{x})=(0,0)$ in figure 4(c). (The value of $\Delta x/\Delta k_{x}$ becomes zero at $\delta _0=|k_y+k_y^{-1}|$, which corresponds to the transition from a single island to two islands, as shown in figure 8.) The coefficient $\eta _{T}$ is determined numerically. As shown in figure 5(b), $\eta _{T}$ can be approximated as

(4.29)\begin{equation} \eta_{T}\approx\sqrt{\left(\frac{1+k_y^2}{\delta_0k_y}\right)^2-1}\end{equation}

at $\beta /\mathcal {C}<-1$, when the mode is well localized in phase space. In this case, (4.20) becomes identical to (4.12).

These results show that the alternative approach adopted here is in agreement with the one we used in § 4.1 if we use the fitting formula (4.25) for $\eta _{R}$ and (4.29) for $\eta _{T}$. If these factors are calculated numerically instead, then the alternative approach is slightly more accurate, as seen in figure 5(c).

4.3 Connection with the Kelvin–Helmholtz instability

The above analysis shows that the TI can be considered as a primary instability modified by ZFs. As seen from figure 5, the growth rate $\gamma _{\textrm {TI}}$ decreases with $|\mathcal {C}|$ in general. Therefore, the TI is very different from the KHI, which develops only in strong ZFs. To study the relation between the TI and the KHI, we numerically solve (3.22a,b) for various $q_{x}$ and $\delta _{0}$ and explore how the mode structure changes with these parameters. The results are shown in figure 6.

Figure 6. Numerical solutions of (3.22a,b) illustrating the relation between the runaway mode and the KH mode at $\beta = 6$, $u = 10$ and various $q_x$ and $\delta _0$. $(a)$ At $\delta _0 = 0$ and $q_x = 1.6$, the unstable mode is the KH mode, which has a global structure as discussed in Zhu et al. (Reference Zhu, Zhou and Dodin2018a,Reference Zhu, Zhou, Ruiz and Dodinc). $(b)$ The KH mode transitions to an ‘intermediate’ mode as $\delta _0$ is increased from $\delta _0=0$ to $\delta _0 = 1.5$ while keeping $q_x = 1.6$ fixed. $(c)$ The corresponding evolution of $\gamma$ with $\delta _0$ at constant $q_x = 1.6$. Blue curves show multiple branches of eigenmodes, but only one branch (the KH mode) is unstable. $(d)$ $\gamma$ as a function of $q_x$ at constant $\delta _0 = 1.5$. As $q_x$ decreases, the intermediate mode analytically continues into the runaway TI mode in figure 3. See the main text for details.

First, consider figure 6(a), which shows a global (not localized) KH mode that corresponds to $q_x = 1.6$ and $\delta _0 = 0$. This KH mode has been discussed in Zhu et al. (Reference Zhu, Zhou and Dodin2018a); it is global because the ZF is small-scale, specifically, $q_x^2 > 1$. Next, let us increase $\delta _0$ from zero up to $\delta _0=1.5$ while keeping $q_x=1.6$ fixed. Then, the original KH mode transforms into an ‘intermediate’ mode shown in figure 6(b). It is not a pure KHI, because dissipation (i.e. non-zero $\delta _0$) is now important, but it is not quite the TI either, because $q_x^2$ is large and the mode localization is less pronounced. Our theory does not apply to such modes, but we have calculated the growth rate numerically as a function of $\delta _0$, as shown in figure 6(c). Finally, with $\delta _0 = 1.5$ fixed, let us reduce $q_x$. The mode localization improves and the instability rates goes down at first, as seen in figure 6(d). But eventually, when $q_x$ has become small enough ($q_x \sim 0.6$), the mode transforms into the runaway mode that we introduced earlier (figure 4) and our theory becomes applicable.

This shows that in principle, the KH mode can be continuously transformed into the runaway mode. However, the KHI and TI are fundamentally different in physical mechanisms, because the TI is due to dissipation and $\gamma _{\textrm {TI}}^{R}$ is determined by $\delta _0$, while the KHI requires a strongly sheared flow and has $\gamma _{\textrm {KHI}} \sim k_y u$. Since typical large-scale ZFs seen in simulations have $q_x^2 \ll 1$, the TI is more relevant to them than the KHI.