$\boldsymbol{v}$. Also, $\nabla^{2}\doteq\partial_{x}^ {2}+\partial_{y}^{2}$ is the Laplacian. Finally, we note that the parameter $\beta$ can be scaled out of equation (\[eq:Equations\_mTHE\]) by replacing $(\phi,t,\hat{D})$ with $(\phi/\beta,\beta t,\hat{D}/\beta)$. Therefore, varying $\beta$ is effectively similar to varying the strength of $\hat{D}$.

The mTHE is “modified” compared to the original Terry–Horton model [@Terry82; @Terry83] in that the following operator $\hat{\alpha}$ is used: $$\hat{\alpha}\varphi=\tilde{\varphi}\doteq\varphi-\langle\varphi\rangle,\label{eq:Equations_alpha}$$ where $\langle\dots\rangle$ is the zonal average given by $$\langle\varphi\rangle\doteq\frac{1}{L_{y}}\int_{0}^{L_{y}}\varphi\,{\rm d}y$$ and $L_{y}$ is the system length along $y$. Equation (\[eq:Equations\_alpha\]) states that electrons respond only to the fluctuation (or DW) part of the potential, $\tilde{\varphi}$, but do not respond to the zonal-averaged (or ZF) part, $\langle\varphi\rangle$ [@St-Onge17; @Hammett93]. The operator $\hat{\delta}$ describes the phase difference between $n$ and $\varphi$ and determines the primary DW instability [@Terry82; @Terry83]. Note that (\[eq:Equations\_mTHE\]) reduces to the modified Hasegawa–Mima equation at $\hat{\delta}=0$ [@Hasegawa77; @Dewar07frontiers], where the total energy is conserved. The DW and the ZF part of the energy (per unit area) are given by $$E_{\rm DW}\doteq\frac{1}{2L_xL_y}\int{\rm d}x\,{\rm d}y\left[(\nabla_{\perp}\tilde{\varphi})^2+\tilde{
    \varphi}^2\right],\quad E_{\rm ZF}\doteq \frac{1}{2L_x}\int{\rm d}x\,\left(\partial_x\langle\varphi\rangle\right)^2,\label{eq:Energies}$$ where $L_x$ is the system length along $x$. Various forms of $\hat{\delta}$ can be used to model different primary instabilities [@Terry82; @Tang78]. Here, we follow [@St-Onge17] and use the following 