Transient Metastability and Selective Decay for the Coherent Zonal Structures in Plasma Drift Wave Turbulence

Abstract

The emergence of persistent zonal structures is studied in freely decaying plasma flows. The plasma turbulence with drift waves can be described qualitatively by the modified Hasegawa–Mima (MHM) model, which is shown to create enhanced zonal jets and more physically relevant features compared with the original Charney–Hasegawa–Mima model. We analyze the generation and stability of the zonal state in the MHM model following the strategy of the selective decay principle. The selective decay and metastable states are defined as critical points of the enstrophy at constant energy. The critical points are first shown to be invariant solutions to the MHM equation with a special emphasis on the zonal modes, but the metastable states consist of a zonal state plus drift waves with a specific smaller wavenumber. Further, it is found with full mathematical rigor that any initial state will converge to some critical point solution at the long-time limit under proper dissipation forms, while the zonal states are the only stable ones. The selective decay process of the solutions can be characterized by the transient visits to several metastable states, then the final convergence to a purely zonal state. The selective decay and metastability properties are confirmed by numerical simulations with distinct initial structures. One highlight in both theory and numerics is the tendency of Landau damping to destabilize the selective decay process.

Appendices

Appendix A: Generalized Dissipation form with Selective Decay

In this appendix, we show the derivation for the general dissipation form that is in agreement with the selective decay principle. As from the main text for the proof of selective decay, the major task is to construct the proper damping operators ${\mathcal {D}}\left( \Delta \right)$ on the right-hand side of (20) so that the Dirichlet quotient $\Lambda \left( t\right)$ stays monotonically decreasing in time. From the first- and second-order selective decay damping forms in Sect. 5.2, it can be summarized that the agreeable dissipation operators for selective decay should follow the general structure

$$\begin{aligned} {\mathcal {D}}\left( \Delta \right) \varphi =-\sum _{j=2}^{L}D_{j}\left[ \left( -\Delta +1\right) ^{j}{\tilde{\varphi }}+\left( -\partial _{x}^{2}\right) ^{j}\overline{\varphi }\right] +D_{1}\left( \Delta \varphi -{\tilde{\varphi }}\right) , \end{aligned}$$

(A.1)

with damping coefficients $D_{j}\ge 0,j\ge 2$. We have shown in (26) that the first-order term above with any constant value $D_{1}$ will not change the value of $\Lambda \left( t\right)$ during its time evolution. The separated damping effects on the fluctuation ${\tilde{\varphi }}$ and zonal mean $\overline{\varphi }$ are reasonable considering the different treatment of the two parts in the MHM model. Next, we derive the dynamical equations for the Dirichlet quotient with a single-order damping j from (A.1).

First from Eq. (20), we have found the dynamical equations for $\Lambda \left( t\right)$ subject to the damping with a single order of the Laplace operator applied on either the full potential function or its fluctuation part ${\tilde{\varphi }}=\varphi -\overline{\varphi }$

$$\begin{aligned} \frac{\mathrm{d}\Lambda }{\mathrm{d}t}= & {} -\frac{d_{j}}{E}U_{j+1},\quad \mathrm {with}\quad {\mathcal {D}}_{j}\varphi =d_{j}\left( -\Delta \right) ^{j}\varphi ,\\ \frac{\mathrm{d}\Lambda }{\mathrm{d}t}= & {} -\frac{d_{j}}{E}{\tilde{U}}_{j+1},\quad \mathrm {with}\quad \mathcal {{\tilde{D}}}_{j}\varphi =d_{j}\left( -\Delta \right) ^{j}{\tilde{\varphi }}. \end{aligned}$$

Then for the generalized damping form in (A.1), we can consider the effects componentwisely through the polynomial expansion of the damping operator

$$\begin{aligned} \left( -\Delta +1\right) ^{j}{\tilde{\varphi }}=\sum _{l=0}^{j}\lambda _{l}\left( -\Delta \right) ^{l}{\tilde{\varphi }}, \end{aligned}$$

with the coefficients $\lambda _{l}=\begin{pmatrix}j\\ l \end{pmatrix}$. Accordingly for the general damping of a single order j,

$$\begin{aligned} -D_{j}\left[ \left( -\Delta +1\right) ^{j}{\tilde{\varphi }}+\left( -\partial _{x}^{2}\right) ^{j}\overline{\varphi }\right] =-D_{j}\left[ \left( -\Delta \right) ^{j}\varphi +\sum _{l=0}^{j-1}\lambda _{l}\left( -\Delta \right) ^{l}{\tilde{\varphi }}\right] , \end{aligned}$$

we get the dynamical equation for $\Lambda \left( t\right)$ in the expansion form by adding up all the component contributions as

$$\begin{aligned} \frac{\mathrm{d}\Lambda }{\mathrm{d}t}=-\frac{D_{j}}{E}\left( U_{j+1}+\sum _{l=0}^{j-1}\lambda _{l}{\tilde{U}}_{l+1}\right) , \end{aligned}$$

(A.2)

where we use the notation $U_{j}=\overline{U}_{j}+{\tilde{U}}_{j}$ from (21) for the contributions from the zonal mean and fluctuation components

$$\begin{aligned} \begin{aligned} \overline{U}_{j}&= \left\| \left( -\partial _{x}^{2}\right) ^{\frac{j}{2}}\overline{\varphi }\right\| _{0}^{2}-\Lambda \left( t\right) \left\| \left( -\partial _{x}^{2}\right) ^{\frac{j-1}{2}}\overline{\varphi }\right\| _{0}^{2},\\ {\tilde{U}}_{j}&= \left\| \left( -\Delta \right) ^{\frac{j}{2}}{\tilde{\varphi }}\right\| _{0}^{2}-\Gamma \left( t\right) \left\| \left( -\Delta \right) ^{\frac{j-1}{2}}{\tilde{\varphi }}\right\| _{0}^{2}. \end{aligned} \end{aligned}$$

Then the task is to reorganize the right-hand side of (A.2) into a summation of nonpositive quantities.

Next, we show the derivation of the recursive relations between the quantities defined in (21) and (23)

$$\begin{aligned} {\tilde{U}}_{j+1}={\tilde{S}}_{j}+\Gamma {\tilde{U}}_{j},\quad U_{j+1}=S_{j}+\Lambda \overline{U}_{j}+\Gamma {\tilde{U}}_{j},\quad S_{1}=-\Lambda U_{1}. \end{aligned}$$

(A.3)

The third relation is already derived in (22) directly from the definition of the Dirichlet quotient. The first two relations are the results from an integration by parts, that is, to get the fluctuation part

$$\begin{aligned} \begin{aligned} {\tilde{U}}_{j+1}=&\int \left| \left( -\Delta \right) ^{\frac{j+1}{2}}{\tilde{\varphi }}\right| ^{2}-\Gamma \left| \left( -\Delta \right) ^{\frac{j}{2}}{\tilde{\varphi }}\right| ^{2}\\ =&\int \left[ \left| \left( -\Delta \right) ^{\frac{j+1}{2}}{\tilde{\varphi }}+\Gamma \left( -\Delta \right) ^{\frac{j+1}{2}-1}{\tilde{\varphi }}\right| ^{2}-2\Gamma \left( \nabla \left( -\Delta \right) ^{\frac{j}{2}}{\tilde{\varphi }}\right) \left( \left( -\Delta \right) ^{\frac{j}{2}-\frac{1}{2}}{\tilde{\varphi }}\right) \right. \\&\left. -\Gamma ^{2}\left| \left( -\Delta \right) ^{\frac{j+1}{2}-1}{\tilde{\varphi }}\right| ^{2}-\Gamma \left| \left( -\Delta \right) ^{\frac{j}{2}}{\tilde{\varphi }}\right| ^{2}\right] \\ =&\int \left| \left( -\Delta \right) ^{\frac{j+1}{2}}{\tilde{\varphi }}+\Gamma \left( -\Delta \right) ^{\frac{j+1}{2}-1}{\tilde{\varphi }}\right| +\int \Gamma \left| \left( -\Delta \right) ^{\frac{j}{2}}{\tilde{\varphi }}\right| ^{2}-\Gamma ^{2}\left| \left( -\Delta \right) ^{\frac{j-1}{2}}{\tilde{\varphi }}\right| ^{2}\\ =&\int \left| \left( -\Delta \right) ^{\frac{j+1}{2}}{\tilde{\varphi }}+\Gamma \left( -\Delta \right) ^{\frac{j-1}{2}}{\tilde{\varphi }}\right| +\Gamma {\tilde{U}}_{j}. \end{aligned} \end{aligned}$$

Above in the second line, remind the notation $\left( -\Delta \right) ^{\frac{1}{2}}=\nabla$, thus integration by parts can be applied for the second term. In a similar fashion, we can find the relation in the zonal mean modes by applying the same trick. Therefore, by introducing the definition for the positive-definite components,

$$\begin{aligned} {\tilde{S}}_{j}=\left\| \left( -\Delta \right) ^{\frac{j+1}{2}}{\tilde{\varphi }}-\Gamma \left( -\Delta \right) ^{\frac{j-1}{2}}{\tilde{\varphi }}\right\| _{0}^{2},\quad {\bar{S}}_{j}=\left\| \left( -\partial _{x}^{2}\right) ^{\frac{j+1}{2}}\overline{\varphi }-\Lambda \left( -\partial _{x}^{2}\right) ^{\frac{j-1}{2}}\overline{\varphi }\right\| _{0}^{2}, \end{aligned}$$

the above two identities are reached. Notice that we have different coefficients $\Lambda \left( t\right) =\Gamma \left( t\right) +1$ in the zonal mean and fluctuation parts.

Now we can derive the final form of the dynamics of (A.2) by applying the identities (A.3) recursively from the original equation. The leading term $U_{j+1}$ can be expanded into all the lower order terms

$$\begin{aligned} \begin{aligned} U_{j+1}&= S_{j}+\sum _{l=1}^{j-1}\left( \Gamma ^{j-l}{\tilde{S}}_{l}+\Lambda ^{j-l}\overline{S}_{l}\right) +\Lambda ^{j}U_{1}+\Gamma ^{j}{\tilde{U}}_{1}-\Lambda ^{j}{\tilde{U}}_{1}\\&= S_{j}+\sum _{l=1}^{j-1}\left( \Gamma ^{j-l}{\tilde{S}}_{l}+\Lambda ^{j-l}\overline{S}_{l}\right) -\Lambda ^{j-1}S_{1}+\Gamma ^{j}{\tilde{U}}_{1}-\Lambda ^{j}{\tilde{U}}_{1}\\&= S_{j}+\sum _{l=2}^{j-1}\left( \Gamma ^{j-l}{\tilde{S}}_{l}+\Lambda ^{j-l}\overline{S}_{l}\right) +\left( \Gamma ^{j-1}-\Lambda ^{j-1}\right) {\tilde{S}}_{1}+\left( \Gamma ^{j}-\Lambda ^{j}\right) {\tilde{U}}_{1}. \end{aligned} \end{aligned}$$

We only need to attend to the last nondefinite term above. Again we can expand the coefficient in the polynomial form and notice $\lambda _{j}=1$

$$\begin{aligned} \left( \Gamma ^{j}-\Lambda ^{j}\right) {\tilde{U}}_{1}=\left[ \Gamma ^{j}-\left( 1+\Gamma \right) ^{j}\right] {\tilde{U}}_{1}=-\sum _{l=0}^{j-1}\lambda _{l}\Gamma ^{l}{\tilde{U}}_{1}. \end{aligned}$$

For each component of the above summation with index l, using the relation ${\tilde{U}}_{j+1}={\tilde{S}}_{j}+\Gamma {\tilde{U}}_{j}$ inversely, we find the further expansion

$$\begin{aligned} \begin{aligned} -\lambda _{l}\Gamma ^{l}{\tilde{U}}_{1}&= \lambda _{l}\Gamma ^{l-1}\left( {\tilde{S}}_{1}-{\tilde{U}}_{2}\right) \\&= \lambda _{l}\Gamma ^{l-1}{\tilde{S}}_{1}+\lambda _{l}\Gamma ^{l-2}\left( {\tilde{S}}_{2}-{\tilde{U}}_{3}\right) \\&= \lambda _{l}\sum _{i=1}^{l}\Gamma ^{l-i}{\tilde{S}}_{i}-\lambda _{l}{\tilde{U}}_{l+1}. \end{aligned} \end{aligned}$$

Therefore, by taking the summation of all the components we get

$$\begin{aligned} -\sum _{l=0}^{j-1}\lambda _{l}\Gamma ^{l}{\tilde{U}}_{1}=\sum _{l=1}^{j-1}\lambda _{l}\sum _{i=1}^{l}\Gamma ^{l-i}{\tilde{S}}_{i}-\sum _{l=0}^{j-1}\lambda _{l}{\tilde{U}}_{l+1}. \end{aligned}$$

Again the first part above is positive definite, and the second part then can be exactly canceled by the rest terms in the full dynamics (A.2). Combining all the above results, we finally reach the form for the total damping contributions from the j-th order dissipation operator

$$\begin{aligned} \begin{aligned} U_{j+1}+\sum _{l=0}^{j-1}\lambda _{l}{\tilde{U}}_{l+1}&= S_{j}+\sum _{l=2}^{j-1}\left( \Gamma ^{j-l}{\tilde{S}}_{l}+\Lambda ^{j-l}\overline{S}_{l}\right) \\&\quad +\left( \Gamma ^{j-1}-\Lambda ^{j-1}\right) {\tilde{S}}_{1}+\sum _{i=1}^{j-1}\Gamma ^{-i}{\tilde{S}}_{i}\sum _{l=i}^{j-1}\lambda _{l}\Gamma ^{l}\\&= S_{j}+\sum _{l=2}^{j-1}\left( \Gamma ^{j-l}{\tilde{S}}_{l} +\Lambda ^{j-l}\overline{S}_{l}\right) \\&\quad +\left( \Lambda ^{j-1}-1\right) \Gamma ^{-1}{\tilde{S}}_{1}+\sum _{i=2}^{j-1}\Gamma ^{-i}{\tilde{S}}_{i}\sum _{l=i}^{j-1}\lambda _{l}\Gamma ^{l}. \end{aligned} \end{aligned}$$

Above in the first line, we just change the order of summation for the last term, and notice that the first term in the summation with $i=1$ in the last summation can be combined with the second term with ${\tilde{S}}_{1}$, that is,

$$\begin{aligned} \Gamma ^{-1}{\tilde{S}}_{1}\sum _{l=1}^{j-1}\lambda _{l}\Gamma ^{l}= & {} \Gamma ^{-1}{\tilde{S}}_{1}\left( \sum _{l=0}^{j}\lambda _{l}\Gamma ^{l}-\Gamma ^{j} -1\right) \\= & {} \left( 1+\Gamma \right) ^{j}\Gamma ^{-1}{\tilde{S}}_{1}-\left( \Gamma ^{j-1}+\Gamma ^{-1}\right) {\tilde{S}}_{1}, \end{aligned}$$

and combining the coefficients

$$\begin{aligned} \left( \Gamma ^{j-1}-\Lambda ^{j-1}\right) +\Gamma ^{-1}\sum _{l=1}^{j-1}\lambda _{l}\Gamma ^{l}=\Lambda ^{j}\Gamma ^{-1}-\Lambda ^{j-1}-\Gamma ^{-1}=\Lambda ^{j-1}\Gamma ^{-1}-\Gamma ^{-1}. \end{aligned}$$

In summary, the final dynamical equation for the Dirichlet quotient $\Lambda \left( t\right)$ under the general j-th oder ($j>1$) damping operator in (A.1) can be found to satisfy the following form

$$\begin{aligned} \frac{\mathrm{d}\Lambda }{\mathrm{d}t}=-\frac{D_{j}}{E}\left[ S_{j}+\sum _{l=2}^{j-1}\left( \Lambda ^{j-l}\overline{S}_{l}+\sum _{i=l}^{j}\lambda _{i}\Gamma ^{i-l}{\tilde{S}}_{l}\right) +\left( \Lambda ^{j-1}-1\right) \Gamma ^{-1}{\tilde{S}}_{1}\right] , \end{aligned}$$

(A.4)

with $\Lambda =\Gamma +1$ and $\lambda _{l}=\begin{pmatrix}j\\ l \end{pmatrix}$ the coefficients before the $x^{l}$ term in the polynomial expansion of $\left( x+1\right) ^{j}$. The right-hand side of the above equation is always negative. Therefore, we conclude that $\Lambda \left( t\right)$ is a monotonically decreasing function in time with a lower bound. The same selective decay principle still applies in the general case.

Appendix B: A Counter-Example with Dissipation on Potential Vorticity Alone that Violates Selective Decay

We have shown in Sect. 5.2 of the main text that the damping form, $D\left( \Delta q-{\tilde{q}}\right)$, gives the convergence to the selective decay state. The second part in the damping form $-D{\tilde{q}}$ includes a pure effect on the fluctuations. Here as a counter-example, we show the second component is essential in maintaining the monotonicity of the Dirichlet quotient in the MHM model.

For the case with only damping on the potential vorticity

$$\begin{aligned} {\mathcal {D}}\varphi =D\Delta q=D\left( \Delta ^{2}\varphi -\Delta {\tilde{\varphi }}\right) , \end{aligned}$$

the dynamical equation for the Dirichlet quotient becomes

$$\begin{aligned} \begin{aligned} \frac{\mathrm{d}\Lambda }{\mathrm{d}t}&=-D\left( \left\| \nabla {\tilde{\zeta }}+\Gamma \nabla {\tilde{\varphi }}\right\| _{0}^{2}+\left\| \partial _{x}\overline{\zeta }+\Lambda \partial _{x}\overline{\varphi }\right\| _{0}^{2}\right) \\&\quad +D\Lambda \left( \left\| \nabla {\tilde{\varphi }}\right\| _{0}^{2}-\Gamma \left\| {\tilde{\varphi }}\right\| _{0}^{2}\right) . \end{aligned} \end{aligned}$$

(B.1)

Without the zonal state $\overline{\varphi }\equiv 0$, it can be seen from Poincaré inequality that the right-hand side of (B.1) is still nonpositive definite just as the CHM case. However, with the effect of a nonzero zonal flow, the term on the second line above is indefinite about its sign. The last indefinite term reflects the interactions between the fluctuation and zonal mean state through the entire Dirichlet quotient $\Lambda \left( t\right)$ that includes ratios of both mean and fluctuation parts. Without the detailed dynamics, it is hard to determine the energy transfer mechanism between the zonal mean and the fluctuation. To show this, consider a small nonzonal perturbation added on a zonal solution

$$\begin{aligned} \varphi _{0}=A\cos \sqrt{\Lambda _{l}+1}x+\epsilon \cos \left( \frac{2\pi }{L}\mathbf {k\cdot x}\right) , \end{aligned}$$

with $\Lambda _{l}=\left( \frac{2\pi }{L}l\right) ^{2}$, $\Lambda _{k}=\left( \frac{2\pi }{L}k\right) ^{2}$, $k>l$ and $\epsilon ^{2}<A^{2}$. Then we can calculate the Dirichlet quotient for this initial state as

$$\begin{aligned} \Lambda _{l}+1<\Lambda \left( 0\right) =\frac{\left( \Lambda _{k}+1\right) ^{2}\epsilon ^{2}+\left( \Lambda _{l}+1\right) ^{2}A^{2}}{\left( \Lambda _{k}+1\right) \epsilon ^{2}+\left( \Lambda _{l}+1\right) A^{2}}<\Lambda _{k}+1. \end{aligned}$$

Substituting the state into the right-hand side of the Eq. (B.1), we have the estimation for the initial transient state dynamics with the state $\varphi _{0}$

$$\begin{aligned} \begin{aligned} \frac{\mathrm{d}\Lambda }{\mathrm{d}t}\ge&-D\left[ \Lambda _{k}\left( \Lambda _{k}+1-\Lambda \left( 0\right) \right) ^{2}\epsilon ^{2}+\left( \Lambda _{l}+1\right) \left( \Lambda _{l}+1-\Lambda \left( 0\right) \right) ^{2}A^{2}\right] \\&+D\Lambda \left( 0\right) \frac{\left( \Lambda _{l}+1\right) \left( \Lambda _{k}-\Lambda _{l}\right) A^{2}}{\left( \Lambda _{k}+1\right) \epsilon ^{2}+\left( \Lambda _{l}+1\right) A^{2}}\epsilon ^{2}\\ \ge&-D\left[ \left( \Lambda _{k}-\Lambda _{l}\right) ^{2}\left( \Lambda _{k}\epsilon ^{2}+\left( \Lambda _{l}+1\right) A^{2}\right) +\frac{\left( \Lambda _{l}+1\right) ^{2}\left( \Lambda _{k}-\Lambda _{l}\right) }{\left( \Lambda _{k}+1\right) +\left( \Lambda _{l}+1\right) }\epsilon ^{2}\right] . \end{aligned} \end{aligned}$$

Therefore, the right-hand side of the equation is larger than zero if

$$\begin{aligned} \epsilon ^{2}>\frac{\left[ \left( \Lambda _{k}+1\right) ^{2}-\left( \Lambda _{l}+1\right) ^{2}\right] \left( \Lambda _{l}+1\right) }{\left[ \left( \Lambda _{l}+1\right) ^{2}-\Lambda _{k}\left( \Lambda _{k}+1\right) \right] \left( \Lambda _{k}+1\right) }A^{2}. \end{aligned}$$

Then by taking the wavenumbers satisfying $\Lambda _{k}\left( \Lambda _{k}+1\right)<\left( \Lambda _{l}+1\right) ^{2}<\left( \Lambda _{k}+1\right) ^{2}$, the Dirichlet quotient will increase in the initial state. Inversely. The larger value of $\Lambda \left( t\right)$ further implies the generation of more higher wavenumber fluctuation modes, thus to push the quotient to even larger values. As a result, this example with special initial state shows that the monotonic decrease of the Dirichlet quotient might be violated with the pure damping form $D\Delta q$. Then the selective decay principle is difficult to guarantee in this case.

Appendix C: Dynamical Convergence to the Zonal Mean Flow

For the convergence to a purely zonal state, we have proved in the main text using the convergence of the infinite integral in the dynamical equation of $\Lambda \left( t\right)$. Here as an alternative approach, we directly show the convergence to zero in the ratio of energy fluctuation from the dynamical equations for the mean and fluctuation parts.

We consider the convergence to a purely zonal state with the dissipation form $-D_{2}\left( -\Delta q+{\tilde{q}}\right)$. In this case, we consider the dynamical equations for the ratios of zonal energy and fluctuation energy

$$\begin{aligned} \frac{{\tilde{E}}\left( t\right) }{E\left( t\right) }+\frac{\overline{E}\left( t\right) }{E\left( t\right) }=1, \end{aligned}$$

with ${\tilde{E}}=\frac{1}{2}\left( \left\| \nabla {\tilde{\varphi }}\right\| _{0}^{2}+\left\| {\tilde{\varphi }}\right\| _{0}^{2}\right)$ the energy in the fluctuation and $\overline{E}=\frac{1}{2}\left\| \partial _{x}\overline{\varphi }\right\| _{0}^{2}$ the energy in the zonal state. First, we have the dynamics for the total energy E and the energy in the fluctuation ${\tilde{E}}$ for this damping form from (6)

$$\begin{aligned} \begin{aligned} \frac{\mathrm{d}E}{\mathrm{d}t}&=-D_{2}\left( \left\| \Delta \varphi \right\| _{0}^{2}+2\left\| \nabla {\tilde{\varphi }}\right\| _{0}^{2}+\left\| {\tilde{\varphi }}\right\| _{0}^{2}\right) ,\\ \frac{\mathrm{d}{\tilde{E}}}{\mathrm{d}t}-\left( \partial _{x}\overline{v},\overline{{\tilde{u}}{\tilde{v}}}\right) _{0}&=-D_{2}\left( \left\| \Delta {\tilde{\varphi }}\right\| _{0}^{2}+2\left\| \nabla {\tilde{\varphi }}\right\| _{0}^{2}+\left\| {\tilde{\varphi }}\right\| _{0}^{2}\right) . \end{aligned} \end{aligned}$$

Notice that there is the interaction term $\left( \partial _{x}\overline{v},\overline{{\tilde{u}}{\tilde{v}}}\right) _{0}$ between the mean and fluctuation due to the nonlinear interaction in the mean energy equation. Then we can find the dynamical equation for the ratio ${\tilde{E}}/E$ through the above two equations

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t}\left( \frac{{\tilde{E}}}{E}\right)= & {} \frac{1}{E^{2}}\left( \dot{{\tilde{E}}}E-{\tilde{E}}{\dot{E}}\right) \\= & {} \frac{1}{E}\left( \partial _{x}\overline{v},\overline{{\tilde{u}}{\tilde{v}}}\right) _{0}\\&-\frac{D_{2}}{2E^{2}}\left[ \left( \left\| \Delta {\tilde{\varphi }}\right\| _{0}^{2}+2\left\| \nabla {\tilde{\varphi }}\right\| _{0}^{2}+\left\| {\tilde{\varphi }}\right\| _{0}^{2}\right) \left( \left\| \nabla {\tilde{\varphi }}\right\| _{0}^{2}+\left\| {\tilde{\varphi }}\right\| _{0}^{2}+\left\| \partial _{x}\overline{\varphi }\right\| _{0}^{2}\right) \right. \\&-\left. \left( \left\| \Delta \varphi \right\| _{0}^{2}+2\left\| \nabla {\tilde{\varphi }}\right\| _{0}^{2}+\left\| {\tilde{\varphi }}\right\| _{0}^{2}\right) \left( \left\| \nabla {\tilde{\varphi }}\right\| _{0}^{2}+\left\| {\tilde{\varphi }}\right\| _{0}^{2}\right) \right] \\= & {} \frac{1}{E}\left( \partial _{x}\overline{v},\overline{{\tilde{u}}{\tilde{v}}}\right) _{0}-\frac{D_{2}}{2E^{2}}\left( W\left\| \partial _{x}\overline{\varphi }\right\| _{0}^{2}-E\left\| \partial _{x}^{2}\overline{\varphi }\right\| _{0}^{2}\right) \\= & {} \frac{1}{E}\left( \partial _{x}\overline{v},\overline{{\tilde{u}}{\tilde{v}}}\right) _{0}-\frac{D_{2}}{2E}\left( \Lambda \left\| \partial _{x}\overline{\varphi }\right\| _{0}^{2}-\left\| \partial _{x}^{2}\overline{\varphi }\right\| _{0}^{2}\right) \\= & {} \frac{1}{E}\left( \partial _{x}\overline{v},\overline{{\tilde{u}}{\tilde{v}}}\right) _{0}-\frac{D_{2}}{E}\left( \Lambda \left( E-{\tilde{E}}\right) -\left( W-{\tilde{W}}\right) \right) \\= & {} \frac{1}{E}\left( \partial _{x}\overline{v},\overline{{\tilde{u}}{\tilde{v}}}\right) _{0}+\frac{D_{2}}{E}\left( \Lambda {\tilde{E}}-{\tilde{W}}\right) \\\le & {} \frac{1}{E}\left( \partial _{x}\overline{v},\overline{{\tilde{u}}{\tilde{v}}}\right) _{0}-D_{2}\left( 1+\Lambda _{1}-\Lambda \left( t\right) \right) \frac{{\tilde{E}}}{E}. \end{aligned}$$

Above we use the relations $\frac{W}{E}=\Lambda \left( t\right)$ and ${\tilde{W}}\ge \left( 1+\Lambda _{1}\right) {\tilde{E}}$. On the other hand, we have the estimation for the nonlinear interaction term

$$\begin{aligned} \frac{1}{E}\left| \left( \partial _{x}^{2}\overline{\varphi },\overline{{\tilde{u}}{\tilde{v}}}\right) _{0}\right| \le \frac{1}{2E}\int \left| \partial _{x}^{2}\overline{\varphi }\right| \left( {\tilde{u}}^{2}+{\tilde{v}}^{2}\right) \le \frac{1}{2E}\left\| \partial _{x}^{2}\overline{\varphi }\right\| _{\infty }\left\| \nabla {\tilde{\varphi }}\right\| _{0}^{2}. \end{aligned}$$

With the selective decay principle satisfied with the eigenvalue $\Lambda _{*}$, we can find that the upper bounds for the total energy and enstrophy decay to zero in the exponential rates

$$\begin{aligned} \begin{aligned} \left\| \nabla \varphi \right\| _{0}&\le \left\| \nabla \varphi \left( 0\right) \right\| _{0}\mathrm{e}^{-D\Lambda _{*}t},\\ \left\| \zeta \right\| _{0}&\le \left\| \zeta \left( 0\right) \right\| _{0}\mathrm{e}^{-D\Lambda _{*}t}. \end{aligned} \end{aligned}$$

Assuming the solution $\overline{\varphi }$ is smooth on a bounded domain, then it implies that the maximum value of zonal vorticity is bounded by any small value, $\left\| \partial _{x}^{2}\overline{\varphi }\right\| _{\infty }\le c$, as time goes on. Therefore for any small value $\epsilon >0$, after large enough time $t>T$, the nonlinear interaction term can always be controlled

$$\begin{aligned} \frac{1}{E}\left| \left( \partial _{x}^{2}\overline{\varphi },\overline{{\tilde{u}}{\tilde{v}}}\right) _{0}\right| \le \frac{c}{2E}\left\| \nabla {\tilde{\varphi }}\right\| _{0}^{2}=\epsilon \frac{{\tilde{E}}}{E}. \end{aligned}$$

The second term in the dynamics of $\Lambda \left( t\right)$ then becomes negative when $1+\Lambda _{1}>\Lambda \left( t\right)$ at some point of the time, and is guaranteed in later times due to the monotonicity of $\Lambda \left( t\right)$. Thus the ratio ${\tilde{E}}/E$ is always decreasing in time after the quotient $\Lambda \left( t\right)$ reaches the value below $\Lambda _{1}+1$.

Notice that we achieve the above result based on the special damping form $-D_{2}\left( -\Delta q+{\tilde{q}}\right)$, thus it is less general than the argument in the main text that can include an additional anti-damping operator as a forcing effect. Still it offers a rigorous proof for the decay of the fluctuation mode and the final convergence to the zonal structure shown in the numerical results.