Critical-layer instability of shallow-water magnetohydrodynamic shear flows

Abstract

In this paper, the instability of shallow-water shear flow with a sheared parallel magnetic field is studied. Waves propagating in such magnetic shear flows encounter critical levels where the phase velocity relative to the basic flow, $c-U(y)$, matches the Alfvén wave velocities $\pm B(y)/\sqrt {\mu \rho }$, based on the local magnetic field $B(y)$, the magnetic permeability $\mu$, and the mass density of the fluid $\rho$. It is shown that when the two critical levels are close to each other, the critical layer can generate an instability. The instability problem is solved, combining asymptotic solutions at large wavenumbers and numerical solutions, and the mechanism of instability explained using the conservation of momentum. For the shallow-water magnetohydrodynamic system, the paper gives the general form of the local differential equation governing such coalescing critical layers for any generic field and flow profiles, and determines precisely how the magnetic field modifies the purely hydrodynamic stability criterion based on the potential vorticity gradient in the critical layer. The curvature of the magnetic field profile, or equivalently the electric current gradient $J' = - B''/\mu$ in the critical layer, is found to play a complementary role in the instability.

1. Introduction

The shallow-water equations are an essential idealised model for the study of large-scale dynamics in geophysical fluid systems. Their relatively simple form has provided significant insights in our understanding of waves, instabilities and turbulence in the oceans and atmosphere. More recently, Gilman (2000) extended the application of the shallow-water model to the dynamics of the solar tachocline by incorporating a magnetic field. In the present study, we will examine the stability of such shallow-water magnetohydrodynamic (MHD) systems.

The stability properties of hydrodynamic shallow-water flows have been studied extensively, and we now know of a number of instabilities with distinctive features. They include Rayleigh's instability, which is related to an inflectional point of the shear flow profile (Blumen, Drazin & Billings 1975), resonant instability generated by interaction between two neutral modes (Satomura 1981; Hayashi & Young 1987), critical-layer instability induced by singularities of neutral modes (Balmforth 1999; Riedinger & Gilbert 2014), and radiative instability caused by waves radiating outwards in an unbounded domain (Ford 1994; Riedinger & Gilbert 2014).

In the context of astrophysical flows, such as in the solar tachocline, magnetic fields are present and will generally modify the stability properties of the flow. Although the elasticity of field lines suggests a stabilising effect, in reality this additional coupling can lead to new modes of instability, as occurs for example in the magnetorotational instability (Balbus & Hawley 1991). The instability of two-dimensional shear flow with a parallel magnetic field has been investigated by a number of researchers. It is well known that a strong magnetic field has a stabilising effect. Gilman (1967), Chandra (1973) and Hughes & Tobias (2001) have shown that modified versions of Howard's (1961) semicircle rule exist when the field is included. The magnetic field reduces the possible domain in which the complex phase velocity may reside, and instability will be suppressed if the magnetic field is sufficiently strong everywhere.

The role of a weaker magnetic field, however, is more subtle, and researchers have found situations where it may have a destabilising effect. Kent (1968), Stern (1963) and Chen & Morrison (1991) have studied the instability problem analytically in the zero wavenumber limit, in which case the dispersion relation reduces to an equation involving a simple integral. They have demonstrated various examples where the magnetic field may destabilise an otherwise stable flow; e.g. a parabolic profile of the field can destabilise a plane Couette flow (Chen & Morrison 1991). Guided by these theoretical studies, Tatsuno & Dorland (2006), Lecoanet et al. (2010) and Heifetz et al. (2015) have computed unstable modes numerically at finite small wavenumbers.

The instability of shallow-water MHD systems has been studied by Mak, Griffiths & Hughes (2016). They consider basic velocity profiles of unstable shear layers and jets, and examine the effect of a uniform magnetic field on these classical instabilities. Their results demonstrate that the field plays mainly a stabilising role. The same semicircle rule of Hughes & Tobias (2001) exists, so a sufficiently strong magnetic field can suppress any instability. Increasing the field always reduces the maximum unstable growth rate, but in some situations may increase the (albeit small) growth rates for long wavelength modes.

The instability analyses cited so far are all based on flows in Cartesian geometry, but for the solar tachocline, spherical geometry is a better representation. Gilman & Fox (1997, 1999) considered two-dimensional MHD flows in a thin spherical shell with the basic differential rotation profile of the Sun and various magnetic fields. They found a ‘joint instability’: either the shear flow or the magnetic field is stable by itself, but the system is unstable when they are present together. Gilman & Dikpati (2002) and Dikpati, Gilman & Rempel (2003) studied the effect of a free surface on these joint instabilities by employing the shallow-water MHD model of Gilman (2000). They show that the free surface has a weak effect on the instability as long as the effective gravity is not too small, but as this parameter is decreased, the instability is eventually completely suppressed as the shell thickness of the reference state tends to zero at certain latitudes. Márquez-Artavia, Jones & Tobias (2017) revealed yet another type of instability for shallow-water MHD flow on a sphere: it is induced purely by the free surface and the magnetic field, and exists even when the basic flow is quiescent or is solid-body rotation. A review of MHD instabilities in spherical shells has been given by Gilman & Cally (2007). The consequences of these instabilities may include transition to turbulence, magnetic reconnection (Cally 2001; Cally, Dikpati & Gilman 2003), and the generation of Rossby waves in the solar tachocline (Dikpati & McIntosh 2020).

In the present study, we investigate the effect of critical levels on the instability of shallow-water MHD systems. Critical levels appear as singularities of steady waves propagating in shear flows if the fluid system has no dissipation. In hydrodynamic problems with flow profile $U(y)$, say, a critical level is a location $y= {y_{{c}}}$ at which the phase velocity of waves matches the basic flow velocity, i.e. $c=U(y)$. When a parallel magnetic field is added, the critical levels become locations $y=y_{B\pm }$ where $c-U(y)$ matches the Alfvén wave velocity $\pm B(y)/\sqrt {\rho \mu }$, with $B(y)$ the magnetic field profile, $\rho$ the mass density of the fluid, and $\mu$ the magnetic permeability. These magnetic critical layers have been found to play crucial roles in a wide range of phenomena and applications, including sunspots (Sakurai, Goossens & Hollweg 1991), solar wind (Chen & Hasegawa 1974), hot Jupiters (Hindle, Bushby & Rogers 2021), and tokamak reactors (Mok & Einaudi 1985).

In hydrodynamic stability theory, critical layers play the key role in driving instabilities in a wide variety of flows. They include the shallow-water flows that we have mentioned, baroclinic flows (Bretherton 1966), stratified flows with horizontal shear (Wang & Balmforth 2018), and, perhaps most famously, wind flowing over water generating surface water waves (Miles 1957). All these instabilities share a common mechanism: the critical layer generates a finite amount of mean-flow momentum (or potential vorticity, equivalently). Thus, given that the total momentum is conserved, the mean-flow momentum in the critical layer must be balanced by unsteady motion in the outer flow, and this balance can be thought of as driving the instability.

Therefore, these instabilities are sometimes referred to as ‘critical-layer instabilities’ (Bretherton 1966; Riedinger & Gilbert 2014). However, the magnetic-field-induced instabilities found in previous studies cannot be understood in terms of this mechanism. For example, in Tatsuno & Dorland (2006), Lecoanet et al. (2010) and Heifetz et al. (2015), it can be inferred that the mean-flow modification is anti-symmetric due to the parity property of the unstable modes. As a result, the mean-flow momentum generated in two critical layers cancels out completely, and so cannot be seen to drive the growth of the outer flow. Instead, the instability may be interpreted as a result of adding more magnetic free energy to the system (Lecoanet et al. 2010) or through the interaction between vorticity waves (Heifetz et al. 2015).

In this paper, we will reveal a new kind of magnetic critical-layer instability in shallow-water MHD systems. It shares a similar mechanism with the hydrodynamic critical-layer instabilities, but magnetic critical levels have distinctive properties. Unlike hydrodynamic instabilities where there is usually a single critical level, here we have two critical levels close to each other, hence their interaction plays a crucial role in the instability. Also, in the hydrodynamic shallow-water system, the potential vorticity (PV) is conserved, and it largely controls the dynamics of the flow. In particular, the sign of the background PV gradient $Q' = -(U'/H)'$ (where $H$ is the depth of shallow water) determines whether the critical layer has a stabilising or destabilising effect (Balmforth 1999; Riedinger & Gilbert 2014). A magnetic field, however, breaks the PV conservation and changes the dynamics significantly. An example of the impact of such a loss of PV conservation has been presented by Dritschel, Diamond & Tobias (2018) for the fundamental problem of the evolution of two-dimensional vortices. We will show that the loss of PV conservation has a dramatic impact on the MHD shallow-water system as well: this is evident in the study of the basic flow of linear shear, i.e. vanishing background vorticity gradient $-U''(y)=0$, and a constant shallow-water depth $H$. In the absence of a magnetic field, PV conservation would imply that there is no critical level at all. The presence of a magnetic field, however, brings back the singular behaviour at the critical levels, hence the possibility of critical-layer instability. The general instability criterion that we obtain combines the PV gradient $Q'$ with the analogous quantity of the electric current gradient $J' = - B''/\mu$ in a key ‘curvature parameter’ $\gamma$ that appears in the local ODE for the critical layer.

Finally, we investigate the mean-flow response of the instability, and show that it is localised strongly in the critical layer. We explain the instability mechanism via the conservation of momentum following the paradigm of Hayashi & Young (1987), which is a balance between the ‘mean momentum’ and the ‘wave momentum’. The mean momentum is just the momentum of mean-flow response, while the wave momentum is the coupling between linear waves of velocity and surface displacements. We show that the mean momentum generated in the critical layer must be balanced by the exponential growth of the wave momentum, and this can be understood as a mechanism for the instability. This mechanism is similar to that of hydrodynamic critical-layer instabilities, but we will show that the magnetic field also controls the mean momentum via the Maxwell stress and the electric current gradient $J'=-B''/\mu$.

The layout of the paper is as follows. In § 2, we present the equations for the problem. The shallow-water MHD system of Gilman (2000) is given in § 2.1, the eigenvalue problem for the linear instability is derived in § 2.2, and the equations for the mean-flow responses and momentum conservation are obtained in § 2.3. In § 3, we present numerical solutions to the instability problem and the mean-flow response for typical basic-flow profiles, and summarise the instability criterion and the instability mechanism. In § 4, we derive the asymptotic solution for the instability problem at large wavenumbers for general basic-flow profiles, and explain the instability mechanism via momentum conservation. We conclude in § 5 and compare the new instability to those discussed in the literature.

2. The governing equations

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 (2000). 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 (1987), 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 (2002) 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. (2016), 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 1967; Chandra 1973; Hughes & Tobias 2001), but note that the semicircle rules in spherical geometry can be different (cf. Watson 1980; Gilman & Fox 1997). 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 (1987), 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 (1997) 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 2014). 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.

3. General results

In this section, we present typical numerical solutions of the eigenvalue problem, and give general conclusions regarding the conditions for the instability. We use two basic flow profiles to present concrete numerical results:

(3.1a,b)\begin{equation} U={-}y,\quad B=\tfrac{1}{2}y(y+1), \end{equation}

and

(3.2a,b)\begin{equation} U={-}y+\tfrac{1}{3}y^{2},\quad B=\tfrac{1}{2}y(y+1). \end{equation}

It is known that the basic-flow vorticity gradient $-U''$ is responsible for hydrodynamic critical-layer instabilities. In order to exclude these instabilities and demonstrate the impact of breaking PV conservation, in the first example we use a profile that has $U''=0$ everywhere. We will show that the magnetic field itself can induce a new kind of instability. In the second example, we demonstrate how a non-zero $U''$ affects the instability. Given that the flow already has a critical-layer instability without the magnetic field (cf. Balmforth 1999), we study how the field modifies it. A sketch of the two profiles with the critical levels identified is shown in figure 1. The field is relatively weak between $y=-1$ and $y=0$, hence the radii in the semicircle rule (2.25) remain positive, which retains the possibility for instability.

Figure 1. A sketch of the basic flow of (a) $U=-y$, $B=y(y+1)/2$, and (b) $U=-y+y^{2}/3$, $B=y(y+1)/2$, with two critical levels $y_{B\pm }$ for ${c_{{r}}} = 0.2$. There is another critical level in (b) that is not labelled.

The numerical solutions of the dispersion relation for the basic state (3.1), i.e. linear $U$ and a parabolic $B$, are shown in figure 2. We plot the real part of the phase velocity ${c_{{r}}}$ and the unstable growth rate $\omega _{{i}}=k{c_{{i}}}$ versus the wavenumber $k$. The solid lines represent normal mode solutions, while the dotted lines represent ‘quasi-modes’, which we will explain later in more detail. We have plotted four modes: L1 and R1 represent the first surface-gravity mode (cf. Balmforth 1999) localised near the left and right boundaries, respectively (see figure 3 for the eigenfunctions of L1). Similarly, L2 and R2 represent the second such modes. In figure 2(a), the dash-dot lines ${c_{{r}}}=1$ and ${c_{{r}}}=-2$ are the conditions that the critical level $y_{B-}$ is on the boundaries $y=-1$ and $y=1$, respectively (cf. figure 1a). In ${c_{{r}}}>1$ or ${c_{{r}}}<-2$, modes have no critical level and they are neutral. In the central region $-2< c_{r}<1$, at least one critical level is inside the domain $-1< y<1$. It is seen that the critical levels destroy most of the normal modes, turning them into quasi-modes. On the segments of solid lines where the normal modes survive, they become unstable, as indicated by the positive growth rates in figure 2(b). For L1 and L2, unstable modes appear at ${c_{{r}}}\approx 0$, whereas for R1 and R2, they appear at ${c_{{r}}}\approx 1$. This is related to the fact that for the profile of (3.1), when ${c_{{r}}}=0$ or ${c_{{r}}}=1$, the two critical levels coalesce at $y=0$ or $y=-1$ where $B=0$ (cf. figure 1a). These instabilities are induced essentially by the critical layers. Since $U''\equiv 0$, they are distinct from the hydrodynamic critical-layer instabilities; the local magnetic field plays the crucial role in the destabilisation, as we will elaborate subsequently. In figure 2(b), we also see two narrow peaks of unstable growth rates, L1-R1 and L1-R2. They are the resonant instabilities induced by two modes with nearly the same phase velocity, i.e. they correspond to the intersections of curves in figure 2(a).

Figure 2. Numerical solution giving (a) ${c_{{r}}}$ and (b) the unstable growth rate $\omega _{{i}}=k {c_{{i}}}$, for the basic flow profile $U=-y$, $B=y(y+1)/2$ at $F=2$. Solid lines represent normal modes governed by (2.23) and (2.24), and dotted lines represent quasi-modes. ‘L1’ represents the ‘first’ mode localised near the ‘left’ boundary, and similar definitions apply for the labelling of other modes. In (a), the dash-dot lines at ${c_{{r}}}=1$ and ${c_{{r}}}=-2$ represent the conditions that critical level $y_{B-}$ is on the boundaries $y=-1$ and $y=1$, respectively.

Figure 3. Eigenfunctions (a) $\hat {h}$, (b) $\hat {v}$, (c) $\hat {a}$ and (d) $\hat {u}$, for the unstable normal mode L1 of figure 2 at $k=3$. In (a) we have also plotted the Wentzel–Kramers–Brillouin (WKB) solution studied in § 4, which uses the asymptotic eigenvalue and has been normalised by $\mathcal {A}=1$ in (4.6). The amplitude of the numerical solution is chosen by fitting to the asymptotic solution in (a).

The eigenfunctions of $\hat {h}$, $\hat {v}$, $\hat {a}$ and $\hat {u}$ for an L1 unstable mode at $k=3$ are shown in figure 3. As stated above, the wavelike structure is localised near the left boundary $y=-1$, representing the surface-gravity mode there, and the two critical levels are close to each other near $y=0$. For $\hat {v}$, $\hat {a}$ and $\hat {u}$, there are very strong amplitude gradients in the critical layer, and because it contains two adjacent critical levels interacting with each other, the critical-layer flow is more distorted than those of hydrodynamic critical layers (see e.g. Drazin & Reid 1982).

The dotted lines in figure 2 represent ‘quasi-modes’, being dotted to indicate that these ‘modes’ are not actual solutions to the eigenvalue problem, but arise only if we deform the path of $y$ into a contour in the complex plane between $y=-1$ and $y=1$. By this means, we obtain non-trivial solutions of (2.23) and (2.24), which are referred to as quasi-modes. Such computations usually appear when we solve an initial-value problem that involves integrals in $y$, the paths of which can be deformed in the complex plane. For large times, a quasi-mode behaves like a decaying normal mode (the decay rates are shown in figure 2b), but also involves the continuous spectrum. In the early stage, however, it can contribute to transient algebraic growth under certain initial conditions (Balmforth, del Castillo-Negrete & Young 1997). For detailed properties and behaviours of quasi-modes, see Briggs, Daugherty & Levy (1970), Balmforth, Llewellyn Smith & Young (2001) and Turner & Gilbert (2007). We will explain briefly the formation of quasi-modes in our problem, and our method to compute them, in § 4.3.

The dispersion relation for the basic state (3.2), i.e. parabolic profiles for both $U$ and $B$, is shown in figure 4. The general features are very similar to those in figure 2. The L1 and L2 unstable modes are again located where ${c_{{r}}}$ is close to zero. However, we notice that there is no longer any instability of the R1 and R2 modes. In addition, the growth rates of the L1 and L2 modes have been enhanced significantly. These are essentially the effects of $U''$ in the critical layer, which we will elaborate later on. Since the basic velocity profile is unstable itself, in figure 5 we plot the instability both with and without the magnetic field. The purely hydrodynamic instability (red dashed curves) has a broader unstable waveband, since there is no additional restriction for the critical-layer instability other than the sign of $U''$ (Balmforth 1999). Thus in the full system (blue curves), the magnetic field has the effect of narrowing the unstable waveband and also seems to inhibit the resonant instability significantly. However, the magnetic field enhances the largest growth rate of the L1 mode.

Figure 4. Numerical solution giving (a) ${c_{{r}}}$ and (b) the unstable growth rate $\omega _{{i}}=k {c_{{i}}}$, for the basic flow profile $U=-y+y^{2}/3$, $B=y(y+1)/2$ at $F=2$. Line styles and notations have the same meanings as in figure 2.

Figure 5. Growth rates of unstable L1 modes for $U=-y+y^{2}/3$ with and without the magnetic field, for $F=2$.

In summary, we see that when the magnetic field is present, the critical-layer instability may arise only when the two critical levels are close to each other in one critical layer where $B\approx 0$. We will show this analytically in § 4.3 using the asymptotic analysis at large $k$, and we find that the closeness is described by

(3.3)\begin{equation} |y_{B+}-y_{B-}|\lesssim \sqrt{\frac{\rm \pi}{k^{3}}\left|\frac{B'B''-U'U''}{U'^{2}-B'^{2}}\right|_{y=y_B}} . \end{equation}

If the two critical levels are well separated, then the magnetic field is found to be stabilising and even hydrodynamic instabilities are suppressed. Note that only the two critical levels closest to the boundary where the surface-gravity mode is localised count. For example, for figure 1(b), if we study the instability of the L1 and L2 modes, then we consider only $y_{B+}$ and $y_{B-}$, and we do not count the additional (unlabelled) critical level further away, since disturbances are much weaker there (similarly to figure 3). Also note that the two critical levels must come from each of $U-c-B=0$ and $U-c+B=0$; if both of them belong to $U-c-B=0$ (or both to $U-c+B=0$), then the critical levels still have a stabilising effect even if they are close to each other.

Our asymptotic analysis also indicates that once the two critical levels are close, the key quantity that determines the instability is $B'B''-U'U''$ in the critical layer. In particular, for modes localised near the left boundary (i.e. L1, L2, etc.), the condition for instability is that in the critical layer,

(3.4)\begin{equation} B'B''-U'U''>0. \end{equation}

By performing a rotation of the domain, the condition for the instability of the modes localised near the right boundary is that in the critical layer,

(3.5)\begin{equation} B'B''-U'U''<0. \end{equation}

These conditions are generalisations of hydrodynamic critical-layer instabilities based on the vorticity gradient $-U''$ (Balmforth 1999; Riedinger & Gilbert 2014), and they can well explain the numerical results that we just presented. For the profile of (3.1), $B'B''-U'U''={1}/{2}$ at $y=0$, so (3.4) is satisfied, and L1 and L2 are destabilised when the critical levels are near $y=0$. Similarly, $B'B''-U'U''=-{1}/{2}$ at $y=-1$, so (3.5) is satisfied, and R1 and R2 are destabilised when the critical levels are near $y=-1$. Since $U''=0$ for this profile, it is the magnetic field that plays the key role in the destabilisation via the current gradient $J' = -B''$. When we include curvature of $U$, in the profile (3.2), $B'B''-U'U''={7}/{6}$ at $y=0$, so unstable modes L1 and L2 again exist, but $B'B''-U'U''={11}/{18}$ at $y=-1$, which violates (3.5), so instability of the modes R1 and R2 no longer occurs, as shown in figure 4. Once (3.4) or (3.5) is satisfied, the largest growth rate increases with the value of $|B'B''-U'U''|$ in the critical layer if the unstable wavenumbers remain similar, as we see in figures 2, 4 and 5 for the L1 modes. We note that our asymptotic analysis is based on large $k$, but we find that for the conditions that we study, it still gives qualitatively good results even when $k$ is of the order of unity.

The numerical results for $\partial _t\Delta U$ and $\partial _t\Delta B$ from (2.45) and (2.47) for the unstable mode of figure 3 are plotted in figure 6, normalised by the exponential growth $\exp (2\omega _{{i}} t)$. The mean-flow responses are strongly localised around the two critical levels. The flow response $\Delta U$ generally exhibits two jets forced in opposite directions. The profile of $\Delta B$ is a little different: it is extended in both directions at each critical level, which may be understood as a result of stretching caused by the mean-flow jet, following Alfvén's theorem.

Figure 6. Mean-flow acceleration of (a) streamwise velocity and (b) magnetic field corresponding to the unstable mode of figure 3.

Another prominent feature of figure 6 is that there is almost no mean-flow or mean-field response outside the critical layer: we find that the amplitudes of both are of magnitude $0.01$ or smaller. In non-magnetic hydrodynamic mean-flow theory, $\partial _t\Delta U\approx 0$ outside the critical layer may be inferred from the ‘non-acceleration rule’ (which also applies when the PV is not zero). This rule states that the mean-flow velocity is not accelerated if the waves are steady and there is no dissipation (see Bühler 2014). Apparently, adding a streamwise magnetic field does not change this in our problem. In Appendix B, we give a mathematical proof that $\partial _t\Delta U$ and $\partial _t\Delta B$ are both zero outside the critical layer in the limit of neutral stability, hence the use of analytical continuation implies $\partial _t\Delta U$ and $\partial _t\Delta B$ are of order $\omega _{{i}}$, which is small. We note that in hydrodynamic wave–mean-flow interaction, the derivation of the non-acceleration rule depends strongly on PV conservation, so it is a little surprising that it still holds when the magnetic field breaks this conservation. Our mathematical derivation in Appendix B shows that in each of the quadratic terms of (2.45)–(2.47), if the wave is steady, i.e. $c_{i}=0$, then the two components of linear waves have phase difference ${\rm \pi} /2$, hence their product is still a wave with zero mean value. We do not have a deeper physical explanation at present, nor can we extend this conclusion to more general flows.

The momentum conservation represented by (2.41) can provide an explanation for the mechanism of the instability. As we see in figure 6, the mean-velocity acceleration $\partial _t\Delta U$ is very strong in the critical layer. We can show that its integral in $y$ over the critical layer has a non-trivial value, which represents a source of mean momentum $M_{m}$. Thus it drives the exponential growth of the outer flow following the conservation of momentum. We will demonstrate details of this mechanism in § 4.4, taking advantage of the large-wavenumber asymptotics.

4. The asymptotic analysis

In order to better understand the instability and obtain conclusions for general smooth profiles, we perform an asymptotic analysis at large wavenumbers. This allows us to derive the instability criteria exhibited in § 3 analytically. We will combine Wentzel–Kramers–Brillouin (WKB) solutions through the bulk of the flow and a local analysis near the critical levels, highlighting the effects of the singularities, and then derive an asymptotic solution for the eigenvalue $c$. The methodology is similar to Riedinger & Gilbert's (2014) analysis of shallow-water instability and Wang & Balmforth's (2018) analysis of strato-rotational instability, but here we have a more complicated critical layer since there are two critical levels inside. We will also study the conservation law of momentum in detail to provide a mechanism for the instability. We will study only the instability induced by critical layers as the principal goal of this paper, though we note that critical layers may also affect the resonant instability, a topic that we leave for further research.

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 1991; Bender & Orszag 2013), 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 (2014).

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 2013) 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 (2014) 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. 1995), 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. 1970). 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. (1970).

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 2014); 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. 1997); 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. 1997; Riedinger & Gilbert 2014).

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. 1997; Riedinger & Gilbert 2014). 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 1985; Booker & Bretherton 1967). 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 (1997), 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 1957; Vekstein 1998; Balmforth 1999; Riedinger & Gilbert 2014), 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.

5. Conclusions and remarks

In this paper, we have studied the linear instability of shallow-water flow with a magnetic field parallel to the basic-flow velocity. We have combined an asymptotic analysis in the short-wavelength limit and a numerical shooting method to solve the instability problem. We paid special attention to the magnetic critical levels, which are located where the Doppler-shifted velocity matches the Alfvén wave velocity, i.e. where $c-U(y)=\pm B(y)$ in our dimensionless system. The critical levels appear as singularities for neutral modes, and generate pronounced wave amplitudes and mean-flow responses in their vicinity, namely in the critical layers. We have shown that when two critical levels are close to each other, they may induce an instability. If the two critical levels are separated, then the magnetic field has a strong stabilising effect.

The centrepiece of our analysis is a local equation for the critical layer, which has two parameters: a ‘separation parameter’ $D$, which represents a rescaled distance between two critical levels, and a ‘curvature parameter’ $\gamma$, representing a combination of the curvature of the field and velocity profile. We have shown that the critical-layer instability may be generated if $D$ is sufficiently small and $\gamma$ has a designated sign. These conditions may be used to study the instability of generalised profiles of velocity $U(y)$ and field $B(y)$. In order for the instability to happen, $B(y)$ needs to be very weak somewhere, the simplest case being where $B(y)$ passes through zero, and then the two critical levels can both reside there, close to each other. As for the profile curvature, the requirement is that the value of $B'B''-U'U''$ in the critical layer should be positive (negative) if the critical layer is to destabilise a surface-gravity mode localised near the left boundary at $y=-1$ (the right boundary at $y=1$). This result generalises that for hydrodynamic instabilities based on the vorticity gradient $-U''$, bringing in the electric current gradient $-B''$ on an equal footing. If these conditions are satisfied, then provided that the surface-gravity mode localised on the boundary exists and there are no other critical levels closer to the boundary, the critical-layer instability will arise.

We have explained the mechanism of the instability via the conservation of momentum, following the framework of Hayashi & Young (1987). There is a balance between the ‘mean momentum’, which consists of the mean-flow modifications, and the ‘wave momentum’, which combines surface displacements and velocities of linear waves. We demonstrate that the critical layer produces a finite amount of mean momentum, which drives the exponential growth of the wave momentum and makes the flow unstable. This mechanism is similar to the critical-layer instability of hydrodynamic shallow-water shear flows (Balmforth 1999; Riedinger & Gilbert 2014), but we note the importance of the magnetic field: the Maxwell stress forces the mean-flow response on an equal footing to the Reynolds stress, and the mean-flow momentum in the critical layer is again controlled by the local value of $B'B''-U'U''$ in the critical layer.

The magnetic field can play a fundamental role in destabilising the flow, so the critical-layer instability reported here is different from the shallow-water MHD instability studied by Mak et al. (2016), which is driven primarily by the hydrodynamic shear. The instability here is also quite different from the field-induced instabilities in two-dimensional shear flows reported in previous studies (Stern 1963; Kent 1968; Chen & Morrison 1991; Tatsuno & Dorland 2006; Lecoanet et al. 2010; Heifetz et al. 2015) owing to the free surface in our problem. The analytical studies of these instabilities (Stern 1963; Kent 1968; Chen & Morrison 1991) were undertaken in the limit of zero wavenumber, and as the numerical solutions of Tatsuno & Dorland (2006), Lecoanet et al. (2010) and Heifetz et al. (2015) confirm, the instability exists when the wavenumbers are small. Our instability, on the other hand, exists for large wavenumbers, which is the feature of the surface-gravity mode. Also, for the numerical studies in these papers, the symmetry of the unstable modes makes the mean-flow modifications anti-symmetric; hence the mean-flow modifications in the two critical layers cancel each other and cannot drive an instability. Momentum conservation in their problems involves only the mean-flow momentum, unlike the balance between the wave momentum and mean momentum in our case.

The magnetic critical layers have also been shown to play important roles in the instability of fluids in spherical geometry: Gilman & Fox (1997, 1999) and Dikpati & Gilman (1999) have studied the instability of fluid in a thin spherical shell with toroidal magnetic field, and they found an energy reservoir for the instability, concentrated around the critical layers. These instabilities have quite a different nature from the instability that we study here: they can arise in the absence of a free surface, so the ‘wave momentum’ that serves as a fundamental element in our instability mechanism does not exist there. In a forthcoming paper, we will show that the instability is strongly related to the spherical geometry of the flow. In particular, the global structure of the unstable mode features a pattern of tilted basic toroidal field (cf. Cally 2001; Cally et al. 2003), and we have found that its interaction with the critical layer drives the instability.

The strong amplitudes in the critical layers suggest that the effects of nonlinearity and diffusion may also be important. Shukhman (1998a,b) has constructed weakly nonlinear theories for the evolution of magnetic critical layers, which we may adopt to extend the study in the present paper.