2.1. The shallow-water MHD model

We study the dynamics of the shallow-water magnetohydrodynamic (SWMHD) model for ideal, perfectly conducting fluid proposed originally by Gilman (Reference Gilman2000). Our main objective is to study the effect of critical levels on stability, and as a starting point we use Cartesian coordinates, which are easier for analysis. Let $(x,y)$ be the horizontal coordinates, and $z$ the vertical direction. Under the assumption that the horizontal scale is much greater than the vertical scale, the leading-order dynamics is characterised by horizontal velocities ${\boldsymbol {u}}_*=[u_*(x,y,t), v_*(x,y,t)]$, which are independent of $z$, and the depth of the shallow water $h_*(x,y,t)$. We use a star subscript following the notation of Hayashi & Young (Reference Hayashi and Young1987), to denote the total quantity, which may include a basic state, a linear disturbance and a mean-flow response. The magnetic field is also dominated by the horizontal field ${\boldsymbol {B}}_*=[B_{1*}(x,y,t),B_{2*}(x,y,t)]$.

The dimensionless governing equations are

(2.1)\begin{gather} \frac{\partial h_*}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left(h_*{\boldsymbol{u}}_*\right)=0, \end{gather}

(2.2)\begin{gather} \frac{\partial {\boldsymbol{u}}_*}{\partial t}+{\boldsymbol{u}}_*\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}}_*={-} F^{{-}2}\, \,{\boldsymbol{\nabla} h_*}+ {\boldsymbol{B}}_*\boldsymbol{\cdot}\boldsymbol{\nabla}{{\boldsymbol{B}}}_*, \end{gather}

(2.3)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}(h_*{\boldsymbol{B}}_*)=0, \end{gather}

(2.4)\begin{gather} \frac{\partial {\boldsymbol{B}}_*}{\partial t}+{\boldsymbol{u}}_*\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{B}}_*={\boldsymbol{B}}_*\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}}_*, \end{gather}

which are the continuity equation, the momentum equation, the divergence-free condition in terms of the horizontal magnetic field, and the induction equation. The depth $h_*$ has been rescaled by the vertical length scale, which may be taken as the depth $H$ of the shallow water. The coordinates $(x,y)$ and ${\boldsymbol {u}}_*$ have been rescaled by the characteristic horizontal length scale $L$ and velocity scale $U_0$, respectively. The Froude number $F$ is defined by $F=U_0/\sqrt {gH}$, and the magnetic field ${\boldsymbol {B}}_*$ has been rescaled by $U_0\sqrt {\mu \rho }$.

Equation (2.3) indicates that the divergence-free condition involves the depth of the shallow water, and we may use this to define a magnetic flux ${\boldsymbol {A}}_*=A_*{\boldsymbol {e}}_z$, where ${\boldsymbol {e}}_z$ is the unit vector in the $z$-direction, such that

(2.5)\begin{equation} h_*{\boldsymbol{B}}_*=\boldsymbol{\nabla}\times {\boldsymbol{A}}_*. \end{equation}

From (2.1) and (2.4), one can show that $A_*$ is conserved following a fluid particle:

(2.6)\begin{equation} \frac{\partial A_*}{\partial t}+{\boldsymbol{u}}_*\boldsymbol{\cdot}\boldsymbol{\nabla} A_*=0, \end{equation}

which provides an alternative description for the induction equation (2.4). For the boundary conditions, we take the normal components of velocity and magnetic field to vanish on boundaries located at dimensionless values $y=\pm 1$:

(2.7)\begin{equation} v_*=B_{2*}=0,\quad \mathrm{at} \ y={\pm} 1. \end{equation}

We also take $h_*$, ${\boldsymbol {u}}_*$, ${\boldsymbol {B}}_*$ and $A_*$ to be periodic in the $x$-direction, with period $2{\rm \pi} /k$, where $k$ is the spatial wavenumber.

Dellar (Reference Dellar2002) has shown that the SWMHD system admits a number of conserved quantities. The most common ones are momentum $M$, energy $E$, and cross-helicity $W$: we have

(2.8a–c)\begin{align} \frac{\mathrm{d}M}{\mathrm{d}t} &= \frac{\mathrm{d}}{\mathrm{d}t}\iint h_*u_*\, \mathrm{d}\kern0.06em x\, \mathrm{d} y=0, \quad \frac{\mathrm{d}W}{\mathrm{d}t} =\frac{\mathrm{d}}{\mathrm{d}t}\iint h_*{\boldsymbol{u}}_*\boldsymbol{\cdot}{\boldsymbol{B}}_*\, \mathrm{d}\kern0.06em x\, \mathrm{d} y=0,\nonumber\\ \frac{\mathrm{d}E}{\mathrm{d}t} &=\frac{\mathrm{d}}{\mathrm{d}t} \iint\frac{1}{2}\left(h_*|{\boldsymbol{u}}_*|^{2}+\frac{|\boldsymbol{\nabla} A_*|^{2}}{h_*}+\frac{h_*^{2}}{F^{2}}\right)\, \mathrm{d}\kern0.06em x\, \mathrm{d} y=0 \end{align}

in one periodic domain. We will study mainly the conservation of momentum. The conservation of the other two quantities will be discussed briefly at the end of § 4.4.

2.2. The linear instability equations

We now consider linear instability for the SWMHD system outlined in the previous subsection. For the basic state, we take the shallow water to have uniform depth when its surface is flat: without loss of generality, we select $h_*= 1$. We take a steady parallel flow and magnetic field pointing in the $x$-direction, with a shear in the $y$-direction: ${\boldsymbol {u}}_*=[U(y),0]$, ${\boldsymbol {B}}_*=[B(y),0]$. According to (2.5), the basic state for $A_*=A(y)$ is determined by $A'(y)=B(y)$. These are all taken to be smooth functions of $y$; there are no internal discontinuities or interfaces present in the systems that we study.

Upon the basic state, we add linear disturbances $h,u,v,a,b_1,b_2$ with

(2.9a–c)\begin{gather} h_*=1+\varepsilon\,h(x,y,t), \quad u_*=U(y)+\varepsilon\,u(x,y,t), \quad v_*=\varepsilon\,v(x,y,t), \end{gather}

(2.9d–f)\begin{gather} A_*=A(y)+\varepsilon\,a(x,y,t),\quad B_{1*}=B(y)+\varepsilon\,b_1(x,y,t), \quad B_{2*}=\varepsilon\,b_2(x,y,t), \end{gather}

where $\varepsilon$ is a small number representing the order of the amplitude of the linear disturbances. We substitute (2.9) into the full SWMHD model (2.1)–(2.7), and the order $\varepsilon$ terms yield the linearised governing equations. The linearised version of (2.5) gives

(2.10a,b)\begin{equation} b_1=a_y-Bh,\quad b_2={-}a_x, \end{equation}

which express the field components in terms of the flux $a$, the subscripts representing partial derivatives. Using (2.10), the linearisation of (2.1), (2.2) and (2.6) yields

(2.11)\begin{gather} h_t+Uh_x+u_x+v_y=0, \end{gather}

(2.12)\begin{gather}u_t+Uu_x+U'v={-}\frac{1}{F^{2}}\,h_x+Ba_{xy}-B'a_x-B^{2}h_x, \end{gather}

(2.13)\begin{gather}v_t+Uv_x={-}\frac{1}{F^{2}}\, h_y-Ba_{xx}, \end{gather}

(2.14)\begin{gather}a_t+Ua_x+Bv=0. \end{gather}

The boundary conditions are

(2.15)\begin{equation} v=a_x=0, \quad \mathrm{at} \ y={\pm} 1. \end{equation}

We seek a normal mode instability

(2.16)\begin{equation} (u,v,h,a)=[\hat{u}(y),\hat{v}(y),\hat{h}(y),\hat{a}(y)]\,{\rm e}^{\mathrm{i}k(x-ct)}+\mathrm{c.c.}, \end{equation}

where $k$ is the wavenumber, $c$ is the complex phase velocity, and c.c. represents the complex conjugate. Substituting (2.16) into (2.11)–(2.15), we obtain

(2.17)\begin{gather} \mathrm{i}k(U-c)\hat{h}+\mathrm{i}k\hat{u}+\hat{v}'=0, \end{gather}

(2.18)\begin{gather}\mathrm{i}k(U-c)\hat{u}+U'\hat{v}={-}\frac{\mathrm{i}k}{F^{2}}\,\hat{h}+\mathrm{i}k(B\hat{a}'-B'\hat{a}-B^{2}\hat{h}), \end{gather}

(2.19)\begin{gather}\mathrm{i}k(U-c)\hat{v}={-}\frac{1}{F^{2}}\,\hat{h}'+k^{2}B\hat{a}, \end{gather}

(2.20)\begin{gather}\mathrm{i}k(U-c)\hat{a}+B\hat{v}=0, \end{gather}

(2.21)\begin{gather}\hat{v}=\mathrm{i}k\hat{a}=0,\quad \mathrm{at} \ y={\pm} 1. \end{gather}

After some algebra, we obtain the relations

(2.22a–c)\begin{equation} \hat{u}={-}\frac{\hat{v}'}{\mathrm{i}k}-(U-c)\hat{h}, \quad \hat{v}={-}\frac{U-c}{\mathrm{i}k[(U-c)^{2}-B^{2}]F^{2}}\, \hat{h}' ,\quad \hat{a}={-}\frac{B}{\mathrm{i}k(U-c)}\, \hat{v}, \end{equation}

and a second-order ODE for $\hat {h}$,

(2.23)\begin{equation} \hat{h}''-\frac{2[(U-c)U'-BB']}{(U-c)^{2}-B^{2}}\, \hat{h}'-k^{2}\left\{1-F^{2}[(U-c)^{2}-B^{2}]\right\} \hat{h}=0, \end{equation}

with boundary conditions

(2.24a,b)\begin{equation} \hat{h}'({-}1)=0,\quad \hat{h}'(1)=0. \end{equation}

Equations (2.23) and (2.24) constitute an eigenvalue problem for the phase velocity $c$, and will be the main problem that we are going to consider. Because all coefficients except $c$ are real, complex phase velocities for normal mode solutions always appear in complex conjugates, i.e. $c={c_{{r}}}+\mathrm {i}{c_{{i}}}$ and $c={c_{{r}}}-\mathrm {i}{c_{{i}}}$. Hence we will consider only normal modes with positive ${c_{{i}}}$, which represent unstable disturbances.

Equations (2.23) and (2.24) are equivalent to the eigenvalue problem of Mak et al. (Reference Mak, Griffiths and Hughes2016), expressed by the equation in terms of $\hat {v}$. They have shown that two semicircle theorems exist for any unstable mode, which we quote below:

(2.25a,b)\begin{align} {c_{{r}}}^{2}+{c_{{i}}}^{2}\leq(U^{2}-B^{2})_{max},\quad \left({c_{{r}}}-\frac{U_{max}+U_{min}}{2}\right)^{2}+{c_{{i}}}^{2}\leq \left(\frac{U_{max}-U_{min}}{2}\right)^{2}-B^{2}_{min}, \end{align}

where $max$ and $min$ indicate the maximum and minimum values among all locations of $y$. The two semicircle rules in (2.25) are the same as those of two-dimensional MHD flow without a free surface (Gilman Reference Gilman1967; Chandra Reference Chandra1973; Hughes & Tobias Reference Hughes and Tobias2001), but note that the semicircle rules in spherical geometry can be different (cf. Watson Reference Watson1980; Gilman & Fox Reference Gilman and Fox1997). Equation (2.25) indicates that for an arbitrary prescribed $U$, if $B$ is sufficiently strong everywhere, then no unstable mode can exist, so the magnetic field must be weak somewhere in the domain for any instability to occur.

The governing ODE (2.23) becomes singular when $c-U=\pm B$. Such locations of $y$ are critical levels, which we define as $y_{B\pm }$, so

(2.26a,b)\begin{equation} U(y_{B+})-c+B(y_{B+})=0,\quad U(y_{B-})-c-B(y_{B-})=0. \end{equation}

The Frobenius solution for $\hat {h}$ around each critical level is

(2.27)\begin{equation} \hat{h}=C_{s\pm}\left[1+\frac{k^{2}}{2}\,(y-y_{B\pm})^{2}\log(y-y_{B\pm})+\cdots\right]+C_{r\pm}\left[(y-y_{B\pm})^{2}+\cdots \right] , \end{equation}

where $C_{s\pm }$ and $C_{r\pm }$ are constants. Although $\hat {h}$ converges as $y\rightarrow y_{B\pm }$, other disturbance components, $\hat {u}$, $\hat {v}$ and $\hat {a}$, all diverge. When the flow is unstable, $c$ has an imaginary part ${c_{{i}}}$ and thus the critical levels $y_{B\pm }$ are also complex: for small $c_{i}$, the imaginary part of (2.26) yields

(2.28)\begin{equation} \operatorname{\mathrm{Im}} y_{B\pm} =\frac{{c_{{i}}}}{\left(U'\pm B'\right)|_{y=\operatorname{\mathrm{Re}} y_{B\pm}}} . \end{equation}

Hence the singularity is avoided, but we nonetheless have locally large amplitudes since ${c_{{i}}}$ is found to be small in our study. We will see that the critical levels play crucial roles in the eigenvalue problem.

If the two critical levels $y_{B\pm }$ coalesce at $y_B$, where $B=0$, then the Frobenius solution about $y_B$ is

(2.29) $$\begin{align} \hat{h}&=C_s\left[1-\frac{k^{2}}{2}\,(y-y_B)^{2}-\left.\frac{U'U''-B'B''}{3({U'^{2}}-{B'^{2}})}\right|_{y=y_B}k^{2}(y-y_B)^{3}\log(y-y_B)+\cdots \right]\nonumber\\ &\quad +C_r\left[(y-y_B)^{3}+\cdots\right]. \end{align}$$

Unlike in (2.27), now the logarithmic singularity of the coalesced critical level depends essentially on the local curvatures of the basic field profiles. Using the relations in (2.22), one can show that for $\hat {u}$, $\hat {v}$ and $\hat {a}$, the singular behaviour of the coalesced critical level is weaker than the separated critical levels. When the two critical levels are close to each other but do not coalesce exactly, $\hat {h}$ has the characteristics of both (2.27) and (2.29), and we will study this problem in detail in the later part of the paper.

Before discussing the solution to the instability problem, we note that in non-magnetic hydrodynamic flows, an important quantity is the PV:

(2.30)\begin{equation} q_*=\frac{v_{*,x}-u_{*,y}}{h_*} . \end{equation}

It is conserved following fluid particles. Linearising it in the same manner, i.e. $q_*=Q(y)+q(x,y,t)$, its value for the basic flow is $Q=-U'(y)$ and for a linear disturbance is

(2.31)\begin{equation} q=v_x-u_y+hU_y. \end{equation}

According to the linear governing equations (2.11)–(2.14), the evolution of $q$ follows

(2.32)\begin{equation} q_t+Uq_x+vQ_y={-}Ba_{xxx}-(Ba_{xy}-B_ya_x-B^{2}h_x)_y. \end{equation}

When the magnetic field is absent, the left-hand side of (2.32) is the material conservation of PV, and it puts a strong constraint on the shallow-water flow. For example, for the linear profile of $U=-y$, which we will study later on, $Q_y\equiv 0$ and hence $q\equiv 0$ for all $x,y,t$ when the normal mode solution (2.16) is applied. This would imply no hydrodynamic critical levels at all. The magnetic field essentially breaks the PV conservation, and as a result, magnetic critical levels with singular behaviour still exist for linear shear flows.

We will solve the eigenvalue problem represented by (2.23) and (2.24) both numerically and asymptotically. The numerical method is a shooting method based on ode15s of Matlab. With an initial guess for $c$, which can be provided by the asymptotic solution, we integrate (2.23) from $y=1$ with $\hat {h}'(1)=0$ to $y=-1$. The value of $\hat {h}'(-1)$ then serves as an error, which provides a correction to $c$ to be reduced by means of Newton iteration. Typical numerical results will be given in § 3. The asymptotic analysis provides approximate analytical solutions for eigenvalues and eigenfunctions at large wavenumbers $k$; details will be elaborated in § 4.

2.3. Mean-flow response and momentum conservation

We further explore the mean-flow response of the system to the instability, an important aspect of nonlinearity. Through quadratic terms, the instability modifies the basic flow and field profiles, which could potentially modify the instability. We will also study the momentum conservation through the mean-flow responses, which can provide a mechanism of the instability.

We extend (2.9) to the next order of $\varepsilon$ to include the mean-flow modifications denoted by $\Delta H$, $\Delta U$, $\Delta V$, $\Delta A$, $\Delta B$ and $\Delta B_2$:

(2.33a,b)\begin{gather} h_*=1+\varepsilon\,h(x,y,t)+\varepsilon^{2}\,\Delta H(y,t), \quad u_*=U(y)+\varepsilon\,u(x,y,t)+\varepsilon^{2}\,\Delta U(y,t), \end{gather}

(2.33c,d)\begin{gather} v_*=\varepsilon\,v(x,y,t)+\varepsilon^{2}\,\Delta V(y,t), \quad A_*=A(y)+\varepsilon\,a(x,y,t)+\varepsilon^{2}\,\Delta A(y,t), \end{gather}

(2.33e,f)\begin{gather} B_{1*}=B(y)+\varepsilon\,b_1(x,y,t)+\varepsilon^{2}\,\Delta B(y,t), \quad B_{2*}=\varepsilon\, b_2(x,y,t)+\varepsilon^{2}\,\Delta B_2(y,t). \end{gather}

We will limit our attention to weak nonlinearity, so that the mean-flow response is weak compared to the linear disturbances. The first harmonics of linear disturbances, the $\exp ({\pm 2\mathrm {i}k(x-ct)})$ waves are also present at order $\varepsilon ^{2}$, but are not of interest in our study.

We denote the zonal average as

(2.34)\begin{equation} \overline{(\cdots)}=\frac{k}{2{\rm \pi}}\int_0^{2{\rm \pi}/k}(\cdots)\, \mathrm{d}\kern0.06em x. \end{equation}

Spatial periodicity implies that the zonal average of linear disturbances and their harmonics, as well as quadratic terms such as $uu_x$, are all zero. Substituting (2.33) into (2.5), selecting the order $\varepsilon ^{2}$ terms and then taking the zonal average, we have

(2.35a,b)\begin{equation} \Delta B={-}B\,\Delta H-\overline{hb_1}+\frac{\partial \Delta A}{\partial y} , \quad \Delta B_2={-}\overline{hb_2}. \end{equation}

Implementing the same procedure for (2.1), (2.2), (2.6) and (2.7), we obtain the mean-flow equations

(2.36)\begin{gather} \frac{\partial\Delta H}{\partial t}+\frac{\partial \Delta V}{\partial y}+\overline{(hv)_y}=0, \end{gather}

(2.37)\begin{gather} \frac{\partial \Delta U}{\partial {t}}+\overline{vu_y}+\Delta V \, U'=\overline{a_x({-}a_{yy}+h_yB+2hB')}, \end{gather}

(2.38)\begin{gather} \frac{\partial \Delta V}{\partial t}+\overline{uv_x}+\overline{vv_y}={-}\frac{1}{F^{2}}\,\frac{\partial\Delta H}{\partial y}+\overline{a_{xx}(hB-2a_y)}, \end{gather}

(2.39)\begin{gather} \frac{\partial\Delta A}{\partial t}+\overline{ua_x}+\overline{va_y}+\Delta V\, B=0, \end{gather}

(2.40)\begin{gather} \Delta V=0, \quad \mathrm{at} \ y={\pm} 1. \end{gather}

Note that the boundary condition for $\Delta A$, i.e. $\partial _x \Delta A=0$, is satisfied automatically since $\Delta A$ is the zonal average independent of $x$.

For momentum conservation, substituting (2.33a,b) into (2.8a) and collecting the $O(\varepsilon ^{2})$ terms, we have

(2.41)\begin{equation} \frac{\mathrm{d}{M_{{w}}}}{\mathrm{d}t}+\frac{\mathrm{d}{M_{{m}}}}{\mathrm{d}t}=0, \end{equation}

where

(2.42a,b)\begin{equation} {M_{{w}}}=\int_{{-}1}^{1} \overline{hu} \,\mathrm{d} y,\quad {M_{{m}}}=\int_{{-}1}^{1}(\Delta U+U\,\Delta H) \, \mathrm{d} y. \end{equation}

Following Hayashi & Young (Reference Hayashi and Young1987), we refer to ${M_{{w}}}$ and ${M_{{m}}}$ as the ‘wave momentum’ and ‘mean momentum’, respectively, since the former is composed of linear disturbance fields, while the latter are mean-flow modifications. It is straightforward to verify that (2.11)–(2.15) and (2.36)–(2.40) guarantee (2.41). The conservation of energy and cross-helicity may also be represented by the balance between the wave and mean components in the same fashion. We will discuss these briefly in § 4.4.

Solving the mean-flow system (2.36)–(2.40) can be complicated in general, but we will see that the instability is weak for the examples that we study, i.e. the growth rate $\omega _{i}=kc_{i}$ is of the order of $0.01$ (cf. figure 2), and this allows us to make significant simplifications and derive relatively compact results. In particular, since the mean-flow responses are driven by terms that are quadratic in the linear disturbances, their time dependence is $\exp (2\omega _{{i}} t)$. Hence in (2.36) and (2.38), the time derivatives of $\Delta H$ and $\Delta V$,

(2.43a,b)\begin{equation} \frac{\partial\Delta H}{\partial t} =2\omega_{{i}}\,\Delta H,\quad \frac{\partial\Delta V}{\partial t}=2\omega_{{i}}\,\Delta V, \end{equation}

are small compared to terms in $\Delta H$ and $\Delta V$ without time derivatives. So we may neglect the time derivative terms and find

(2.44a,b)\begin{equation} \Delta V={-}\overline{hv},\quad \frac{\partial \Delta H}{\partial y}={-}F^{2}\left[\overline{a_{xx}({-}hB+2a_y)}+\overline{uv_x}+\overline{vv_y}\right]. \end{equation}

In (2.37) and (2.39), however, there are no terms in $\Delta U$ and $\Delta A$ without time derivatives, so the quadratic terms drive $\partial _t\Delta U$ and $\partial _t\Delta A$ directly. Substituting in (2.44a), we find

(2.45)\begin{gather} \frac{\partial \Delta U}{\partial t}=\overline{vq}+\overline{a_x({-}a_{yy}+Bh_y+2B'h)}, \end{gather}

(2.46)\begin{gather} \frac{\partial \Delta A}{\partial t}={-}\overline{ua_x}+\overline{v({-}a_y+Bh)},\end{gather}

where the PV $q$ is defined in (2.31). Combining (2.46) and (2.35a), we find

(2.47)\begin{equation} \frac{\partial \Delta B}{\partial t}=\overline{({-}ua_x-va_y+Bvh)_y}, \end{equation}

with $O(\omega _{{i}})$ terms again neglected. Mean-flow equations similar to (2.45) and (2.47) have been derived by Gilman & Fox (Reference Gilman and Fox1997) in spherical coordinates. When the field is switched off, (2.45) becomes $\partial _t\Delta U=\overline {vq}$, which is the classical result for the mean-flow response in hydrodynamic flows (cf. Bühler Reference Bühler2014). In that case, $q\equiv 0$ rendered from PV conservation in our linear shear flow would simply indicate no mean-flow response at all. The magnetic field, however, breaks this simple state of affairs fundamentally, as we will see subsequently.

From (2.43a), (2.44b) and (2.45), we can deduce that the time derivative of the mean surface displacement $\partial _t \Delta H$ is order $O(\omega _i)$ smaller than that of the mean velocity $\partial _t\Delta U$, hence we will neglect the former in the time derivative of the mean momentum and let

(2.48)\begin{equation} \frac{\mathrm{d}M_{m}}{\mathrm{d}t}=\int_{{-}1}^{1}\frac{\partial\Delta U}{\partial t}\,\mathrm{d} y. \end{equation}

We will present the numerical solution of (2.45) and (2.47) in § 3 to show the acceleration of mean velocity and field. We will also analyse the momentum conservation (2.41) in § 4 to give a mechanism for the instability.

4.1. WKB solutions

We rewrite (2.23) as

(4.1)\begin{equation} \hat{h}''-\frac{[(U-c)^{2}-B^{2}]'}{(U-c)^{2}-B^{2}}\, \hat{h}'+l^{2}\hat{h}=0, \end{equation}

where

(4.2)\begin{equation} l^{2}={-}\lambda^{2}={-}k^{2}\left\{1-F^{2}[(U-c)^{2}-B^{2}]\right\}. \end{equation}

In the short-wavelength limit $k\gg 1$, we have $l,\lambda \gg 1$, thus (4.1) has WKB solutions. Since ${c_{{i}}}$ is a small number, we may assume $c\approx {c_{{r}}}$ in (4.1) and (4.2), hence $l^{2}$ and $\lambda ^{2}$ are approximately real, as long as we are not close to the critical levels. The height field $\hat {h}$ is wavelike when $l^{2}>0$ and evanescent when $\lambda ^{2}>0$; $l$ and $\lambda$ represent the approximate wavenumber and the exponential decay rate, respectively.

We take the modes localised near the left boundary, i.e. L1 and L2 in figure 2, as an example for the asymptotic analysis. The distribution of $l^{2}$ for the eigenfunction of figure 3 is shown in figure 7. There is a turning point located at $y = {y_{{t}}}$ where $l^{2}=0$. Hence $\hat {h}$ is wavelike in $-1< y<{y_{{t}}}$ and evanescent in $y>y_{t}$, which can also be seen in figure 3. The two critical levels $y_{B\pm }$ are in the evanescent region. They render a thin critical layer where the WKB solution fails. For convenience, we define their midpoint as

(4.3)\begin{equation} y_{B}=\frac{y_{B+}+y_{B-}}{2} , \end{equation}

representing the centre of the critical layer.

Figure 7. Distribution of $l^{2}$ for the mode of figure 3.

For general basic flow profiles, it is also necessary for the instability that $l^{2}>0$ near the boundary, so that $\hat {h}$ is wavelike and the surface-gravity mode can exist. For modes localised near the left boundary, this means

(4.4a,b)\begin{equation} U-c>\sqrt{B^{2}+F^{{-}2}} \quad \mathrm{or}\quad U-c<{-}\sqrt{B^{2}+F^{{-}2}}, \quad \mathrm{at} \ y={-}1. \end{equation}

It is also necessary that there is no critical level other than $y_{B\pm }$ between $y_B$ and $y=-1$ (otherwise that critical level would be the dominant one to determine the instability property and hence the subject to study). The continuous functions $U(y)-c\pm B(y)$ have designated signs at $y=-1$ indicated by (4.4), but are nearly zero at $y=y_B$, hence to guarantee that they have no other zeros in between, the signs of their derivatives at $y=y_{B}$ are also fixed; that is,

(4.5a,b)\begin{equation} U'\pm B'<0 \quad \mathrm{or} \quad U'\pm B'>0,\quad \mathrm{at} \ y=\operatorname{\mathrm{Re}} y_B, \end{equation}

corresponding to (4.4a) and (4.4b), respectively. For either case, $|U'|>|B'|$ at $y=\operatorname {\mathrm {Re}} y_B$. Whether the critical layer has (4.5a) or (4.5b) holding will determine the sign of $\operatorname {\mathrm {Im}} y_{B\pm }$ (see (2.28)), and therefore results in a similar derivation with numerous sign changes. We will use the combination of (4.4a) and (4.5a) for our derivation, which is the case for figure 1. The other situation of (4.4b) and (4.5b) will be noted briefly, and in fact, the resulting condition for instability is the same. The flow field in $y>y_B$ is not important since the disturbance is weak there. There could be other turning points or critical levels in $y>y_B$, as long as they are not close to $y_B$.

In $y>y_{B}$, we consider the WKB solution of (4.1) that decays exponentially:

(4.6)\begin{equation} \hat{h}=\mathcal{A} \sqrt{\frac{(U-c)^{2}-B^{2}}{\lambda}} \exp\left(-\int_{y_{B}}^{y}\lambda(y')\, \mathrm{d} y'\right), \end{equation}

where $\mathcal {A}$ is an arbitrary constant. In ${y_{{t}}}< y< y_B$, because of the critical layer, both exponential solutions exist, and $\hat {h}$ is expressed by

(4.7)\begin{equation} \hat{h}=\mathcal{A}_-\hat{h}_-{+}\mathcal{A}_+\hat{h}_+, \end{equation}

where $\mathcal {A}_-$ and $\mathcal {A}_+$ are constants, and $\hat {h}_-$ and $\hat {h}_+$ are the exponentially decaying and growing solutions, respectively:

(4.8)\begin{equation} \hat{h}_\pm = \sqrt{\frac{(U-c)^{2}-B^{2}}{\lambda}} \exp\left({\pm}\int_{y_{B}}^{y}\lambda(y') \,\mathrm{d} y'\right). \end{equation}

In $-1< y<{y_{{t}}}$, (4.7) is still applicable but we need to find the corresponding wavelike solutions of $\hat {h}_-$ and $\hat {h}_+$. Following the standard procedure to match across the turning point ${y_{{t}}}$ via Airy functions (cf. Hinch Reference Hinch1991; Bender & Orszag Reference Bender and Orszag2013), we find

(4.9a)\begin{gather} \hat{h}_- = 2 \sqrt{\frac{(U-c)^{2}-B^{2}}{l}}\,\varPsi \cos\left[\int_{y_{t}}^{y}l(y')\, \mathrm{d} y'+\frac{\rm \pi}{4}\right],\end{gather}

(4.9b)\begin{gather} \hat{h}_+ = \sqrt{\frac{(U-c)^{2}-B^{2}}{l}}\,\frac{1}{\varPsi} \cos\left[\int_{y_{t}}^{y}l(y')\, \mathrm{d} y'-\frac{\rm \pi}{4}\right], \end{gather}

where

(4.10)\begin{equation} \varPsi=\exp\left(\int_{{y_{{t}}}}^{y_{B}}\lambda\, \mathrm{d} y\right) \end{equation}

represents the exponential gain (loss) of the amplitude of $\hat {h}_-$ ($\hat {h}_+$) from $y_{B}$ to ${y_{{t}}}$. We will need to determine the relation between $\mathcal {A}$, $\mathcal {A}_-$ and $\mathcal {A}_+$ through analysis of the critical layer that connects (4.6) and (4.7).

4.2. Local solution in the critical layer

In the critical layer, we introduce a stretched coordinate

(4.11a,b)\begin{equation} \eta=\frac{y-y_B}{\delta},\quad \delta=k^{{-}1}\ll 1, \end{equation}

based on the short-wave limit. To derive a local equation for (2.23), we Taylor expand $U-c\pm B$ around their zeros $y_{B\pm }$, then substitute in the local coordinate (4.11) and take the leading two orders of $\delta$. After some algebra, we arrive at the local equation

(4.12)\begin{equation} \hat{h}_{\eta\eta}+\left(\gamma-\frac{2\eta}{\eta^{2}-D^{2}}\right)\hat{h}_\eta-\hat{h}=0, \end{equation}

with the two parameters

(4.13a,b)\begin{equation} D=\frac{y_{B+}-y_{B-}}{2\delta},\quad \gamma=\left.- \delta\, \frac{B'B''-U'U''}{B'^{2}-U'^{2}}\right|_{y=y_{B}}. \end{equation}

The quantity $D$ represents a rescaled distance between the two critical levels that are located at $\eta =\pm D$ in the local equation, hence we refer to $D$ as the ‘separation parameter’. The parameter $\gamma$ is determined by the curvature of the profiles of the basic velocity and magnetic field, and is therefore referred to as the ‘curvature parameter’. When the magnetic field vanishes, the separation parameter is $D=0$ and we recover the hydrodynamic critical level at $\eta =0$, with the curvature parameter $\gamma$ determined by the vorticity gradient $Q' = -U''$. When a magnetic field is involved, the current gradient $J' = -B''$ appears on an equal footing. It should be noted that although the curvature parameter $\gamma = O(\delta )$ is algebraically small, it plays a crucial role in the singularity and instability. This can be understood easily for the hydrodynamic case $D=0$: $\gamma$ determines the strength of the singularity at $\eta =0$, and the singularity becomes removable if $\gamma$ is absent. Similar arguments have also been given by Riedinger & Gilbert (Reference Riedinger and Gilbert2014).

To provide the connection condition for the outer WKB solution, we consider the behaviour of $\hat {h}(\eta )$ in an intermediate regime $y-y_B=O(\delta ^{{1}/{2}})$, or $\eta = O(\delta ^{-{1}/{2}})$. Applying the method of dominant balance (Bender & Orszag Reference Bender and Orszag2013) to (4.12), we find

(4.14a,b)\begin{equation} \hat{h}\sim\left\{\begin{array}{ll} \alpha\eta\,{\rm e}^{-\eta}, & \eta = O(\delta ^{-{1}/{2}}),\ \eta>0,\\[2pt] \alpha_-(-\eta)\,{\rm e}^{-\eta}+\alpha_+(-\eta)\,{\rm e}^{\eta}, & \eta=O(\delta ^{-{1}/{2}}), \ \eta<0, \end{array}\right. \end{equation}

where $\alpha$, $\alpha _-$ and $\alpha _+$ are constants. For positive $\eta$, only the decaying solution is included in accordance with the WKB solution. We will see that the key quantity that controls the instability is

(4.15)\begin{equation} \beta=\operatorname{\mathrm{Im}} \left(\frac{\alpha_+}{\alpha_-}\right), \end{equation}

representing the phase difference between the growing and decaying amplitude for negative $\eta$, and we will study its properties in detail. The relation between $\alpha _-$ and $\alpha$, representing the amplitudes on two sides of the critical layer, is noted in Appendix A.

When $D=0$, i.e. two critical levels overlap exactly, $\hat {h}$ can be represented by confluent hypergeometric functions, and we can derive the analytical connection condition. We can show that under the condition (4.5a) and $c_{i}\rightarrow 0^{+}$, in the limit of small $\gamma = O(\delta )$,

(4.16)\begin{equation} \beta|_{D=0}=\frac{\rm \pi}{2} \gamma. \end{equation}

This result has been derived by Riedinger & Gilbert (Reference Riedinger and Gilbert2014) for hydrodynamic shallow-water critical layers, but continues to apply when we include the field curvature in $\gamma$. Note that (4.16) is the leading-order solution for small $\gamma$; a more precise solution for finite $\gamma$ is given in (A5c) in Appendix A. When $D\neq 0$, we have not been able to derive the connection formula analytically. Equation (4.12) may be converted to Heun's equation (cf. Arscott et al. Reference Arscott, Slavyanov, Schmidt, Wolf, Maroni and Duval1995), but still, we have not been able to find analytical connection formulae for Heun's functions in the literature. Therefore, we will resort to numerical solutions of (4.12).

There are some issues to which we should pay attention when solving the equation numerically. First, we need to make sure that we select the correct branch when passing the logarithmic branch points. The Frobenius solutions of (4.12) around $\eta =\mp D$ are

(4.17)\begin{equation} \hat{h}=C_{s\pm}\left[1+\tfrac{1}{2}(\eta\pm D)^{2}\log(\eta\pm D)+\cdots\right]+C_{r\pm}\left[\frac{(\eta\pm D)^{2}}{k^{2}}+\cdots\right], \end{equation}

which are equivalent to (2.27). By definition, $y_{B\pm }$, $\eta$, $D$ and $\gamma$ are all slightly complex due the small growth rate ${c_{{i}}}$. Our aim is to approach the limit ${c_{{i}}}\rightarrow 0^{+}$ if numerically possible, so as to draw parallel conclusions with the analytical relation (4.16). But we also need to make sure that in (4.17) we select the same branch as (2.27) for the logarithm function, and therefore we require that $\operatorname {\mathrm {Im}} (\eta \pm D)$ and $\operatorname {\mathrm {Im}} (y-y_{B\pm })$ have the same sign. According to (2.28), the quantities $\operatorname {\mathrm {Im}} y_{B\pm }$ are negative for an unstable mode if we prescribe the condition (4.5a). So for convenience, we assume that $D$ and $\gamma$ are real but add a very small positive imaginary part to $\eta$ when integrating (4.12) numerically. In practice, we use $\operatorname {\mathrm {Im}}\eta = 10^{-6}$–$10^{-8}$; using different values of $\operatorname {\mathrm {Im}} \eta$ in this range does not affect the solution for $\hat {h}$ or $\beta$.

There is also a technical issue regarding the large $|\eta |$ limits. The asymptotic behaviour (4.14) becomes precise when $|\eta |$ is very large, but numerically, $\alpha _+(-\eta )\,{\rm e}^{\eta }$ becomes too small compared to $\alpha _-(-\eta )\,{\rm e}^{-\eta }$ for very large $-\eta$, and the value of the former may not be stored precisely in $\hat {h}$ within the usual numerical precision. To tackle this difficulty, we develop a numerical method that separately computes the solution corresponding to the limit of $\alpha _+(-\eta )\,{\rm e}^{\eta }$, based on shooting from both sides of the domain. Details of this method are presented in Appendix A.

A sample solution for $\hat {h}$ is given in figure 8(a): $\hat {h}$ in general decays exponentially, but becomes flat at the two critical levels, as predicted by the Frobenius solution (4.17). There is a phase shift of the decaying amplitude across the critical layer, which is rendered by a complex $\alpha _-/\alpha$. The numerical solution for the key quantity $\beta$ is shown in figure 8(b) in solid curves, as a function of the separation parameter $D$ and the curvature parameter $\gamma$. It is seen clearly that $\beta$ increases with $\gamma$, but decreases with $D$. In our subsequent analysis, we are concerned about the situation where $\beta$ is positive, as this is the condition that instability may arise. When the two critical levels overlap, i.e. $D=0$, our analytical solution (4.16) indicates that the curvature parameter $\gamma$ has to be positive for $\beta$ to be positive, which is also seen in the figure. When $D\neq 0$, we found that the decrease of $\beta$ with $D$ can be well fitted by a quadratic function

(4.18)\begin{equation} \beta\approx \beta |_{D=0}-2D^{2}. \end{equation}

The results of (4.18) are plotted in figure 8(b) in dashed lines. Hence a positive $\beta$ requires $|D|$ to be relatively small. According to the value of $\beta |_{D=0}$ provided by (4.16), we need $|D|\lesssim \sqrt {{\rm \pi} \gamma }/2$.

Figure 8. (a) The local solution of $\hat {h}$ at $D=0.3$, $\gamma =0.3$, and (b) $\beta$ as a function of $D$ and $\gamma$ for the condition (4.5a). In (b), solid lines represent the numerical solution for $\beta$, and dashed lines show the results of the empirical fitting (4.18).

4.3. Matching and eigenvalues

We take the limit $y\rightarrow y_B$ for the WKB solutions (4.6) and (4.7), then match them to the inner solution (4.14). This provides the connection conditions

(4.19a–d)\begin{equation} \frac{\mathcal{A}}{\alpha}=\sqrt{\frac{k^{3}}{U_B'^{2}-B_B'^{2}}}\, , \quad \frac{\mathcal{A}_-}{\mathcal{A}}=\frac{\alpha_-}{\alpha}, \quad\frac{\mathcal{A}_+}{\mathcal{A}-} =\frac{\alpha_+}{\alpha_-} , \quad \operatorname{\mathrm{Im}} \left(\frac{\mathcal{A}_+}{\mathcal{A}_-}\right)=\beta. \end{equation}

Finally, we incorporate the boundary conditions to determine the eigenvalue $c$. For modes localised near the left boundary, the eigenfunction $\hat {h}$ has an exponential decay structure, so its amplitude is much larger near the left boundary $y=-1$ than near the right boundary $y=1$. Therefore, the eigenvalue is primarily determined by the boundary condition at $y=-1$, namely $\hat {h}'(-1)=0$, which is

(4.20)\begin{equation} \hat{h}_-'({-}1)+\frac{\mathcal{A}_+}{\mathcal{A}_-}\, \hat{h}_+'({-}1)=0. \end{equation}

Since $|\hat {h}_-|\gg |\hat {h}_+|$ given $\varPsi \gg 1$, the dominant balance of (4.20) is

(4.21)\begin{equation} \hat{h}'_-({-}1)=0. \end{equation}

Substituting in (4.9a) and taking the leading order in terms of $l\gg 1$, we can derive an integral dispersion relation

(4.22)\begin{equation} \int_{{-}1}^{y_{t}}l(y',{c_{{r}}})\, \mathrm{d} y'= \left(n{\rm \pi}-\tfrac{3}{4}{\rm \pi}\right),\quad n=1,2,\ldots. \end{equation}

Equation (4.22) implicitly determines $c={c_{{r}}}$, which is real. In the absence of a critical layer, it determines the phase velocity of a neutral surface-gravity mode in the large wavenumber limit. The integer $n$ represents the index of the mode; i.e. L1 and L2 in figure 2 correspond to $n=1$ and $n=2$.

The exponentially small $\hat {h}_+$ in (4.20), however, can make $c$ complex and render an unstable flow. Suppose that including the $\hat {h}_+$ term causes a small correction to the phase velocity: $c_{r}\rightarrow c_{r}+\Delta c$. Then from (4.20) and (4.21), we have

(4.23)\begin{equation} \left.\frac{\partial}{\partial c} \hat{h}_-'({-}1)\right|_{c={c_{{r}}}}\Delta c+\left.\frac{\mathcal{A}_+}{\mathcal{A}_-}\, \hat{h}'_+({-}1)\right|_{c={c_{{r}}}}=0. \end{equation}

It is now apparent that $\operatorname {\mathrm {Im}} (\mathcal {A}_+/\mathcal {A}_-)=\beta$ can generate an imaginary part ${c_{{i}}}$ of $\Delta c$, which demonstrates how the critical layer can make the flow unstable. Substituting (4.9a) into (4.23), again taking the leading order of $l$, after some algebra we find

(4.24a,b)\begin{equation} {c_{{i}}}={-}\frac{\beta}{2\varPsi^{2}\displaystyle\int_{{-}1}^{{y_{{t}}}}\left(\displaystyle\frac{\partial l}{\partial c}\right)_{c={c_{{r}}}}\mathrm{d} y},\quad \frac{\partial l}{\partial c}={-}\frac{kF^{2}(U-c)}{\sqrt{F^{2}[(U-c)^{2}-B^{2}]-1}} . \end{equation}

As prescribed in (4.4a), $U>{c_{{r}}}$ near the left boundary $y=-1$, so an unstable mode ${c_{{i}}}>0$ requires $\beta >0$. The WKB solution for $\hat {h}$ is compared with the numerical solution in figure 3(a). The asymptotic solutions of ${c_{{r}}}$ and $c_{i}$ for the dispersion relation of figure 4 are plotted in figure 9. For both L1 and L2 modes, the asymptotic solution becomes inaccurate for smaller $k$, because for negative ${c_{{r}}}$, the additional critical level in figure 1(b) moves close to $y_{B\pm }$, and all critical levels may disappear. These effects are not included in the asymptotic analysis, and make the error of $c_{i}$ very sensitive to the error of ${c_{{r}}}$. But overall, the asymptotic analysis for $k\gg 1$ gives qualitatively good predictions to the eigenfunction and the eigenvalue, and in fact, it still works well when $k$ is of the order of unity, as we see in figure 9.

Figure 9. Comparison between the numerical solution (solid lines) and asymptotic solution (dashed lines) for $c_{r}$ and $c_{i}$ for the dispersion relation in figure 4. Dotted lines for $c_{i}$ are the solutions predicted by momentum conservation in (4.37).

Hence the crucial condition that the critical layer destabilises the mode localised near the left boundary is $\beta >0$, for the situation under the conditions (4.4a) and (4.5a). As shown in figure 8, a positive $\beta$ requires two conditions: the curvature parameter $\gamma$ should be positive and the separation parameter $D$ should be relatively small. The first condition requires that the curvature of the profiles of the basic velocity and field have designated signs in the critical layer, that is

(4.25)\begin{equation} B'B''-U'U''>0 ,\quad \mathrm{at} \ y=y_B.\end{equation}

The second condition, on the other hand, requires that the two critical levels should be close to each other:

(4.26)\begin{equation} |y_{B+}-y_{B-}|\lesssim \frac{\sqrt{{\rm \pi}\gamma}}{k} \end{equation}

approximately from (4.18). It is straightforward to check that for the case of fields satisfying (4.4b) and (4.5b), the condition for the critical layer to destabilise the mode near the left boundary is the same. We need to pay attention only to the fact that in this case, $U-{c_{{r}}}<0$ in (4.24), but we replace $\beta$ by $-\beta$ in (4.16) and figure 8(b), a consequence of the sign of $\operatorname {\mathrm {Im}} y_{B\pm }$ reversing (see (2.28)). From (4.13b), (4.16), (4.18) and (4.24), we also see that the maximum value of the growth rate is proportional to $\gamma$ if the other quantities in (4.24) remain similar, thus the growth rate is positively correlated with the value of $B'B''-U'U''$ in the critical layer. The conditions and properties of the instability exhibited in § 3 are therefore derived analytically.

If ${c_{{i}}}$ predicted by (4.24) becomes negative, then it does not represent a decaying normal mode. This is because ${c_{{i}}}<0$ implies $\operatorname {\mathrm {Im}} \eta <0$, which is inconsistent with our choice of $\operatorname {\mathrm {Im}} \eta \rightarrow 0^{+}$ in the critical layer based on the prescription of ${c_{{i}}}>0$ (cf. (2.28), (4.11), (4.17) and related discussion). In such situations, the normal mode cannot exist: it is destroyed by the critical layer. In the setting of an initial-value problem, however, the solution of the Laplace-transformed variable involves integrals in $y$, and one may deform the integral contour in the complex plane (Briggs et al. Reference Briggs, Daugherty and Levy1970). If we choose $\operatorname {\mathrm {Im}} y>\operatorname {\mathrm {Im}} y_{B\pm }>0$ on the contour, then $\operatorname {\mathrm {Im}} \eta$ is still positive even when $c_{i}<0$. In this way, we may recover the solution to the eigenvalue problem (2.23) and (2.24) for complex $y$, which is the quasi-mode to which we referred earlier in § 3. For the numerical solution of the quasi-modes shown in figures 2 and 4, we have applied the shooting method on the complex contour:

(4.27)\begin{equation} y=w+\mathrm{i}\alpha(1-w^{2}),\quad -1< w<1, \alpha>0, \alpha,w \ \mathrm{real}, \end{equation}

which is an arc connecting $y=-1$ and $y=1$ in the complex plane with $\operatorname {\mathrm {Im}} y>0$. For the theoretical foundations of recovering quasi-modes from initial value problems, see Briggs et al. (Reference Briggs, Daugherty and Levy1970).

4.4. Implications of momentum conservation

Finally, we study the momentum conservation in (2.41) to provide a mechanism for the critical-layer instability in this system. In studies of hydrodynamic instabilities, the conservation law has been used to represent the signature of different types of instabilities. For example: Rayleigh's instability features the conservation of mean momentum ${M_{{m}}}$ (Bühler Reference Bühler2014); the resonance instability of shallow-water flows with linear shear velocities features the cancellation of wave momentum ${M_{{w}}}$ with opposite signs (Balmforth et al. Reference Balmforth, del Castillo-Negrete and Young1997); and the critical-layer instability of shallow-water flows features a balance between the wave momentum ${M_{{w}}}$ and the mean momentum ${M_{{m}}}$ (Balmforth et al. Reference Balmforth, del Castillo-Negrete and Young1997; Riedinger & Gilbert Reference Riedinger and Gilbert2014).

The wave momentum is determined by the global structure of the unstable mode:

(4.28a,b)\begin{equation} {M_{{w}}}=2 \hat{M}_{{w}}\, {\rm e}^{2\omega_{{i}} t},\quad\hat{M}_{{w}}=\int_{{-}1}^{1} \tfrac{1}{2}(\hat{h}\hat{u}^{*}+\hat{u}\hat{h}^{*})\, \mathrm{d} y. \end{equation}

Although $\hat {u}$ is locally strong in the critical layer, because the critical layer is too thin, $\hat {h}\hat {u}^{*}+\hat {u}\hat {h}^{*}$ is still not strong enough to give a significant contribution to the integral from the critical layer. Outside the critical layer, the exponential decay of the eigenfunctions implies that the main contribution to the integral comes from the region $y\in [-1, {y_{{t}}}]$. Hence in the large-$k$ limit, using (2.22) and (2.23), one can derive that

(4.29)\begin{equation} \hat{M}_{{w}} \simeq{-}\int_{{-}1}^{{y_{{t}}}}\frac{U-c_{r}}{F^{2}[(U-c_{r})^{2}-B^{2}]}\, |\hat{h}|^{2}\, \mathrm{d} y. \end{equation}

Given conditions (4.4a) and (4.5a), we have $\hat {M}_{{w}}<0$, and hence for an unstable mode,

(4.30)\begin{equation} \frac{\mathrm{d} {M_{{w}}}}{\mathrm{d}t}=4\omega_{{i}}\hat{M}_{{w}} \, {\rm e}^{2\omega_{{i}} t}<0. \end{equation}

For the mean-flow momentum, we have found that the value of $\partial _t \Delta U$ is at $O(k^{-1}\omega _{i}|\hat {v}|^{2})$ outside the critical layer, which is $O(k^{-1})$ smaller compared to the wave momentum in the large-$k$ limit. Details of the computation are given in Appendix B. Hence it is the critical-layer mean-flow acceleration that balances the wave momentum, which is the same as for hydrodynamic shallow-water instabilities (Balmforth et al. Reference Balmforth, del Castillo-Negrete and Young1997; Riedinger & Gilbert Reference Riedinger and Gilbert2014). Let $\varDelta$ be half of the critical-layer thickness; then we choose $\varDelta = O(\delta ^{{1}/{2}})$ as we did in (4.14). Substituting (2.45), (2.31) and (2.11) into (2.48), we derive

(4.31)\begin{align} \frac{\mathrm{d}{M_{{m}}}}{\mathrm{d}t}&\simeq\int_{y_{B}-\varDelta}^{y_B+\varDelta }\frac{\partial \Delta U}{\partial t}\, \mathrm{d} y\nonumber\\ &=\int_{y_B-\varDelta}^{y_B+\varDelta}\left[-\overline{(uv)_y}+\overline{U'vh-u(h_t+Uh_x)}-\overline{(a_xa_y)_y}+\overline{a_x(Bh_y+2B'h)}\right] \mathrm{d} y. \end{align}

From the Frobenius solutions (2.27) or (2.29), together with (2.22), one can show that the dominant terms of (4.31) are $-\overline {(vu)_y}$ and $-\overline {(a_xa_y)_y}$, and other terms that involve the surface displacement $h$ are much smaller. Therefore,

(4.32)\begin{equation} \left.\frac{\mathrm{d}M_{m}}{\mathrm{d}t}\simeq{-}\left.\overline{(uv+a_xa_y)}\right. \right|_{y_B-\varDelta}^{y_B+\varDelta}. \end{equation}

Equation (4.32) indicates that the mean-flow response is determined by the jump of Reynolds stress $\overline {uv}$ and Maxwell stress $\overline {a_xa_y}$ across the critical layer, and it extends the result of the mean-flow response of hydrodynamic critical layers determined by the jump of the Reynolds stress (for example, Killworth & McIntyre Reference Killworth and McIntyre1985; Booker & Bretherton Reference Booker and Bretherton1967). At the edge of the critical layer $y_B\pm \varDelta$, using (2.22) with $\delta \ll \varDelta \ll 1$, we can show that the following relation holds:

(4.33)\begin{equation} \overline{a_xa_y}\simeq{-}\left(\frac{B_B'}{U_B'}\right)^{2}\overline{uv}. \end{equation}

So the Maxwell stress has the opposite sign to the Reynolds stress, similar to what was found by Gilman & Fox (Reference Gilman and Fox1997), and given $|B'_B|<|U'_B|$ in our problem (see (4.5a)), the former is weaker. We note that we should not interpret the Maxwell stress as the overall effect of the magnetic field, because the Reynolds stress is also controlled by the value of $B'B''-U'U''$ in the critical layer.

To find the value of $\mathrm {d}_t{M_{{m}}}$, we first substitute the normal mode solution, (2.16) and (2.22), into (4.32) to represent it in terms of $\hat {h}$: in the limit of large $k$,

(4.34)\begin{equation} \frac{\mathrm{d}{M_{{m}}}}{\mathrm{d}t} \simeq {\rm e}^{2\omega_{{i}} t}\left.\left[\frac{2 \operatorname{\mathrm{Im}} (\hat{h}^{'*}\hat{h}'')}{k^{3}F^{4}\,|(U-c)^{2}-B^{2}|}\right]\right|_{y_B-\varDelta}^{y_B+\varDelta}. \end{equation}

Then we use the asymptotic solution (4.14) to compute the value of $\hat {h}$ at $y_B\pm \varDelta$. Using the property that $D\delta \ll \varDelta \ll 1$ and the relations in (4.19), after some algebra, we derive a compact result in terms of $\beta$:

(4.35)\begin{equation} \frac{\mathrm{d}{M_{{m}}}}{\mathrm{d}t}\simeq 4 {\rm e}^{2\omega_{{i}} t}\, \frac{ |\mathcal{A}_-|^{2}\beta}{k F^{4}} . \end{equation}

When the magnetic field is switched off, from (2.22b), (4.16), (4.19) and (A5b), together with $\hat {h}|_{\eta =0}=\alpha$ for the local solution, a property of the relevant confluent hypergeometric function, it can be shown that (4.35) reduces to

(4.36)\begin{equation} \frac{\mathrm{d}{M_{{m}}}}{\mathrm{d}t}\simeq {\rm e}^{2\omega_{{i}} t}\, \frac{2{\rm \pi}|\hat{v}_c|^{2}U_c''}{k|U_c'|}, \end{equation}

where now we use the subscript $c$ to represent the hydrodynamic critical level $y_c$. Equation (4.36) is the classical result for the mean-flow momentum forced in the critical layer (Miles Reference Miles1957; Vekstein Reference Vekstein1998; Balmforth Reference Balmforth1999; Riedinger & Gilbert Reference Riedinger and Gilbert2014), with the important implication that its value is proportional to the local vorticity gradient $-U_c''$. So again, we have generalised this result to MHD flows through use of the phase difference parameter $\beta$. To balance $\mathrm {d}_t{M_{{w}}}<0$, (4.35) should be positive, which again requires $\beta >0$. The balance between (4.30) and (4.35) yields another expression for $c_{i}$:

(4.37)\begin{equation} c_{i}={-}\frac{|\mathcal{A}_-|^{2}\beta}{ k^{2} F^{4}\hat{M}_{{w}}} . \end{equation}

Therefore, momentum conservation indicates that the exponential growth of the wave momentum is driven by the mean-flow acceleration in the critical layer, and this serves as a mechanism for the instability. The result of (4.37) is also plotted in figures 9(b,c) using dotted lines, again giving qualitatively good predictions. Note that here we are using the precise numerical solutions of ${c_{{r}}}$ for the computation of (4.37), and that is why it is much more precise than the results of (4.24) at smaller wavenumbers.

If we study the conservation law of the energy $E$ and cross-helicity $W$ shown in (2.8b,c), then under the assumption of small $\varDelta$ and $\omega _{i}$, we can derive that the contributions of the critical layer to their mean components are

(4.38)\begin{gather} \frac{\mathrm{d}E_{{m},B}}{\mathrm{d}t}=U_B\int_{y_B-\varDelta}^{y_B+\varDelta}\frac{\partial\Delta U}{\partial t}\, \mathrm{d} y+B_B\int_{y_B-\varDelta}^{y_B+\varDelta}\frac{\partial \Delta B}{\partial t}\, \mathrm{d} y, \end{gather}

(4.39)\begin{gather} \frac{\mathrm{d}W_{{m},B}}{\mathrm{d}t}=U_B\int_{y_B-\varDelta}^{y_B+\varDelta}\frac{\partial\Delta B}{\partial t}\, \mathrm{d} y+B_B\int_{y_B-\varDelta}^{y_B+\varDelta}\frac{\partial \Delta U}{\partial t}\, \mathrm{d} y. \end{gather}

For the integral of the mean-field time derivative $\partial _t\Delta B$, according to (2.46) and (2.47),

(4.40)\begin{equation} \int_{y_B-\varDelta}^{y_B+\varDelta}\frac{\partial \Delta B}{\partial t}\, \mathrm{d} y=\left.\frac{\partial \Delta A}{\partial t} \right|_{y_B-\varDelta}^{y_B+\varDelta}, \end{equation}

but we show in (B2) that the value of $\partial _t \Delta A=0$ is nearly zero outside the critical layer, so the value of (4.40) is negligible. In our problem, $B_B\approx 0$ in the critical layer (the condition for the two critical levels to be close), so the mean cross-helicity generated in the critical layer given by (4.39) is negligible. Similarly, our basic flow also has $U_{B}\approx 0$ for the L modes, so the mean energy in the critical layer is also negligible. Hence in conclusion, in the conservation laws of energy and cross-helicity, the critical layer does not provide a source of mean-flow components that drives the growth of the outer flow. Instead, the conservation is achieved by the cancellation of various wave and mean components outside the critical layer. We will not investigate these balances further in this paper.