An analytical study of the MHD clamshell instability on a sphere

Abstract

This paper studies the instability of two-dimensional magnetohydrodynamic systems on a sphere using analytical methods. The underlying flow consists of a zonal differential rotation and a toroidal magnetic field is present. Semicircle rules that prescribe the possible domain of the wave velocity in the complex plane for general flow and field profiles are derived. The paper then sets out an analytical study of the ‘clamshell instability’, which features field lines on the two hemispheres tilting in opposite directions (Cally, Sol. Phys., vol. 199, 2001, pp. 231–249). An asymptotic solution for the instability problem is derived for the limit of weak shear of the zonal flow, via the method of matched asymptotic expansions. It is shown that when the zonal flow is solid body rotation, there exists a neutral mode that tilts the magnetic field lines, referred to as the ‘tilting mode’. A weak shear of the zonal flow excites the critical layer of the tilting mode, which reverses the tilting direction to form the clamshell pattern and induces the instability. The asymptotic solution provides insights into properties of the instability for a range of flow and field profiles. A remarkable feature is that the magnetic field affects the instability only through its local behaviour in the critical layer.

1. Introduction

Magnetohydrodynamic (MHD) instability is of significant importance to astrophysical flows. Magnetic fields are ubiquitous in stars and planets, and although the magnetic field can act as a restoring force, it can also destabilise the fluid, resulting in turbulence and flow rearrangement.

There are numerous MHD instabilities, differing in geometry and parameter regime. In this study, our particular focus is on two-dimensional MHD instability on a sphere. This has important applications to the solar tachocline, which is a thin transition layer between the Sun's radiative interior and the outer convection zone. The tachocline couples these two regions that have distinct properties, and it plays a pivotal role in solar physics. In particular, this thin shear layer is believed to be the seat of the solar dynamo (see, for example, Charbonneau 2014; Brun & Browning 2017). A strong azimuthal magnetic field generated in the tachocline rises up to the photosphere via magnetic buoyancy, and generates a variety of surface phenomena including sunspots, coronal loops and flares. The MHD instabilities in the solar tachocline can significantly modify the magnetic field and, thus, have a strong impact on the subsequent surface phenomena.

In this context, there have been numerous studies of MHD instabilities of zonal flow in a thin spherical shell coupled with a toroidal magnetic field. The differential rotation profile of the Sun may be modelled as the angular velocity (Newton & Nunn 1951)

(1.1)\begin{equation} \varOmega=r-s\mu^2, \end{equation}

where $r$ is the angular velocity at the equator, $s$ is the shear rate and $\mu =\cos \theta$, $\theta$ being the colatitude in spherical polar coordinates. The parameters $r$ and $s$ are both positive for solar differential rotation, as the Sun rotates faster at the equator. For hydrodynamic flows without magnetic fields, Watson (1981) found that the zonal flow (1.1) becomes unstable when $s/r>0.29$. For the Sun, $s/r=0.12\sim 0.17$ (see, for example, Gough 2007) and the flow is hydrodynamically stable. When a magnetic field is added, however, Gilman & Fox (1997) have shown that instabilities may be present for $s/r$ smaller than 0.29: these instabilities require weaker shear, implying that the magnetic field has a destabilising effect. Gilman & Fox (1997) named such instabilities ‘joint instabilities’ since they only arise when both the hydrodynamic shear and the magnetic field, which are stable separately, exist together. Joint instabilities exist for a wide range of toroidal magnetic fields, including broad field profiles with single or multiple nodes (Gilman & Fox 1997, 1999), and magnetic field bands localised at various latitudes (Dikpati & Gilman 1999; Gilman & Dikpati 2000).

The ‘singular points’ of the unstable modes play an important role in this joint instability. These points are located where the phase velocity of the mode relative to the basic zonal flow matches the characteristic velocity of the Alfvén waves. Here the wave equation becomes singular for ideal fluids; weak effects of viscosity or magnetic resistivity, unsteadiness or nonlinearity remove the singularity, but the disturbances still exhibit strong amplitudes locally. In hydrodynamic stability theory the singular points and their vicinity are named ‘critical levels’ and ‘critical layers’ (Drazin & Reid 1982), terminology we will adopt here. Gilman & Fox (1999) and Dikpati & Gilman (1999) have found that the disturbances change dramatically across the critical layers, and that the critical layers largely determine the spatial structure of the unstable modes. Moreover, various stresses that contribute to the energy of the instability are concentrated in the critical layers, suggesting that they are responsible for driving the instability.

Cally (2001) and Cally, Dikpati & Gilman (2003) computed the nonlinear evolution of the joint instability for two-dimensional flow on a sphere. For strong and broad magnetic field profiles, the magnetic field lines feature a ‘clamshell’ pattern in the early stage: the field lines are tilted in opposite directions on the two hemispheres, and they named this the ‘clamshell instability’. In the later nonlinear evolution, the field lines tilt over $90^\circ$ and reconnect at the equator. These behaviours remain similar when the more realistic physics of density stratification and vertical shear are added, though they produce more subtle vertical structures (Miesch, Gilman & Dikpati 2007). If an external force that maintains the zonal flow and the poloidal field is added (Miesch 2007), then the field lines do not tilt over. Instead, the clamshell pattern is maintained together with a strong mean toroidal field. The mean field and the unstable mode dominate in turn in a quasi-periodic manner.

In the present study we undertake analytical investigations to understand two- dimensional MHD instability on a sphere. We derive the semicircle rules that prescribe the domain of the complex phase velocities for general profiles of a zonal flow and a toroidal magnetic field. A semicircle rule was first derived by Howard (1961) for hydrodynamic instability in Cartesian geometry. It states that, for any unstable mode in an inviscid parallel shear flow $U(y)$, the phase velocity $c=c_{r}+\mathrm {i}c_{i}$ must lie in a semicircle with centre $(c_{r},c_{i})=( (U_{max}+U_{min})/2,0)$ and radius $(U_{max}-U_{min})/2$ above the real axis. It has become a celebrated result of hydrodynamic instability for its generality and simple form. Watson (1981), Thuburn & Haynes (1996) and Sasaki, Takehiro & Yamada (2012) extended the theory to hydrodynamic instability in spherical geometry, showing that this introduces additional terms in the radius of the semicircle. On the other hand, Howard & Gupta (1962), Gilman (1967), Chandra (1973), Cally (2000), Hughes & Tobias (2001) and Deguchi (2021) derived semicircle rules for MHD instability in Cartesian geometry. Their results indicate that the magnetic field can reduce the radii of the semicircles. Consequently, if the magnetic field is strong enough everywhere, the semicircles will disappear and stability is guaranteed. In this paper we derive semicircle rules for MHD instability in spherical geometry. An interesting phenomenon is that in this geometry, the magnetic field can increase the radii of the semicircles, something that never happens in the Cartesian case.

We also develop an asymptotic analysis for the clamshell instability. The limit that we consider is when the shear of the zonal flow is relatively weak, a typical situation in solar physics: in (1.1), $s/r$ is a small number for solar differential rotation. For the magnetic field, we consider the broad and strong field profiles that are typical for the clamshell instability (Cally 2001). The location where the magnetic field profile passes through zero, the node of the profile, is where the magnetic critical level sits and this plays a fundamental role in the instability (see, e.g. Gilman & Fox 1999). Here, in our first study of the problem, we limit our attention to the situation where there is only one such node of the field profile and we assume that its gradient there is non-zero. We derive the solution to the instability problem for general profiles from this family.

Our asymptotic analysis can help us better understand the clamshell instability in several ways. First, it can provide an analytical explanation for the mechanism of the instability: we show that the instability is caused by the interaction between the global tilting motion of the magnetic field and the critical level of the mode located at the node of the field profile. Second, it can give insights for general profiles of magnetic field and shear flows: which types of profiles are unstable and which are not. Finally, it can easily tackle the neutral stability limit where numerical solutions often suffer from resolution difficulties caused by singularities emerging in the corresponding eigenmodes.

The organisation of the paper is as follows. In § 2 we present the governing equations of two-dimensional MHD flow on a sphere, and the corresponding equations for the linear instability problem. The equations for the mean-flow response and angular momentum conservation are also given. In § 3 we derive the semicircle rules for the complex phase velocity and we discuss their applications. In § 4 we undertake an asymptotic analysis for the clamshell instability. In § 4.1 we show that the tilting mode exists for solid body rotation, and that a weak shear can excite its critical levels. Then we solve the eigenvalue problem in § 4.2 by matching the tilting mode and critical layer. The results of the eigenvalue problem, which yield a number of general conclusions, are discussed in §§ 5.1–5.4, and the conservation of angular momentum is shown to provide a mechanism for the instability in § 5.5. Concluding remarks are given in § 6.

2. Governing equations

In this section we present the equations of two-dimensional MHD on a sphere and derive the equations governing the linear instability problem. We consider the motion of an incompressible, inviscid, perfectly electrically conducting fluid with density $\rho$ and magnetic permeability $\mu$. The flow is on a sphere of radius $R$ with a characteristic velocity $U_0$. The MHD equations for the dimensionless velocity ${\boldsymbol {u}}$, magnetic field ${\boldsymbol {B}}$ and pressure $p$ are

(2.1)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}} =0, \end{gather}

(2.2)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{B}}=0, \end{gather}

(2.3)\begin{gather} \frac{\partial{\boldsymbol{u}}}{\partial t}+({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}) {\boldsymbol{u}}={-}\boldsymbol{\nabla} \left(p+\frac{1}{2} {|B|^2}\right)+({\boldsymbol{B}}\boldsymbol{\cdot} \boldsymbol{\nabla} ){\boldsymbol{B}}, \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}

where the length, time, velocity, magnetic field and pressure have been non- dimensionalised by $R$, $R/U_0$, $U_0$, $U_0\sqrt {\mu \rho }$ and $\rho U_0^2$, respectively. We use spherical polar coordinates $(r,\theta,\phi )$, where $r$ is the radius, $\theta$ is the co-latitude and $\phi$ is the longitude, and the corresponding unit vectors are $({\boldsymbol {e}}_r, {\boldsymbol {e}}_\theta, {\boldsymbol {e}}_\phi )$. We take the system to be two dimensional, as for a thin spherical shell, so that the velocity and magnetic field have no radial component (that is, in the ${\boldsymbol {e}}_r$ direction), and are independent of the radial coordinate $r$.

On writing ${\boldsymbol {u}}= v{\boldsymbol {e}}_\theta + u{\boldsymbol {e}}_\phi$ and ${\boldsymbol {B}}= b{\boldsymbol {e}}_\theta + a{\boldsymbol {e}}_\phi$, where $u$, $v$, $a$ and $b$ are functions of $\theta$, $\phi$ and $t$, (2.1)–(2.4) become

(2.5)\begin{gather} \frac{\partial u}{\partial \phi}+\frac{\partial}{\partial \theta}(v\sin \theta)=0, \end{gather}

(2.6)\begin{gather} \frac{\partial a}{\partial \phi}+\frac{\partial}{\partial \theta}(b\sin \theta)=0, \end{gather}

(2.7)\begin{gather} \frac{\partial u}{\partial t}+\frac{u}{\sin \theta}\frac{\partial u}{\partial \phi}+v\, \frac{\partial u}{\partial \theta}+\frac{uv\cos \theta}{\sin \theta}={-}\frac{1}{\sin \theta}\frac{\partial p}{\partial\phi}-\frac{b}{\sin \theta}\frac{\partial b}{\partial \phi}+b\,\frac{\partial a}{\partial \theta}+\frac{ab \cos \theta}{\sin \theta}, \end{gather}

(2.8)\begin{gather} \frac{\partial v}{\partial t}+\frac{u}{\sin \theta}\frac{\partial v}{\partial \phi}+v\,\frac{\partial v}{\partial \theta}-\frac{u^2\cos \theta}{\sin \theta}={-}\frac{\partial p}{\partial \theta}-a\,\frac{\partial a}{\partial \theta}+\frac{a}{\sin \theta}\frac{\partial b}{\partial \phi}-\frac{a^2\cos\theta}{\sin \theta}, \end{gather}

(2.9)\begin{gather} \frac{\partial a}{\partial t}-\frac{\partial}{\partial \theta}(ub-va)=0, \end{gather}

(2.10)\begin{gather} \frac{\partial b}{\partial t}+\frac{1}{\sin \theta}\frac{\partial}{\partial \phi}(ub-va)=0. \end{gather}

We consider an axisymmetric basic state consisting of a zonal flow and a toroidal magnetic field: $u=U$, $v=0$, $a=A$, $b=0$, where $U$ and $A$ vary with co-latitude $\theta$ but are independent of longitude $\phi$. The basic state pressure $P$ is therefore governed by the balance

(2.11)\begin{equation} \frac{\partial}{\partial \theta}\left(P+\frac{1}{2}{A^2}\right)= (U^2-A^2)\cot \theta. \end{equation}

We then study the instability of this state to small disturbances, i.e.

(2.12a–e)\begin{equation} u=U+u_\ell,\quad v=v_\ell,\quad a=A+a_\ell,\quad b=b_\ell, \quad p=P+p_\ell, \end{equation}

where the $\ell$ subscript denotes disturbances of linear instability. Substituting (2.12a–e) into (2.5)–(2.10) and linearising gives

(2.13)\begin{gather} \frac{\partial u_\ell}{\partial \phi}+\frac{\partial}{\partial \theta}(v_\ell\sin \theta)=0, \end{gather}

(2.14)\begin{gather} \frac{\partial a_\ell}{\partial \phi}+\frac{\partial}{\partial \theta}(b_\ell\sin \theta)=0, \end{gather}

(2.15)\begin{gather} \frac{\partial u_\ell}{\partial t}+\frac{U}{\sin \theta}\frac{\partial u_\ell}{\partial \phi}+\frac{\mathrm{d}U}{\mathrm{d}\theta}\, v_\ell+\frac{U\cos\theta}{\sin \theta}\, v_\ell={-}\frac{1}{\sin\theta}\frac{\partial p_\ell}{\partial \phi}+\frac{\mathrm{d}A}{\mathrm{d}\theta}\, b_\ell+\frac{A\cos\theta}{\sin \theta}\, b_\ell, \end{gather}

(2.16)\begin{gather} \frac{\partial v_\ell}{\partial t}+\frac{U}{\sin\theta}\frac{\partial v_\ell}{\partial \phi}-\frac{2U\cos \theta}{\sin \theta}\, u_\ell={-}\frac{\partial}{\partial\theta}\left(p_\ell+Aa_\ell\right)+\frac{A}{\sin\theta}\frac{\partial b_\ell}{\partial\phi}-\frac{2A\cos\theta}{\sin\theta}\, a_\ell, \end{gather}

(2.17)\begin{gather} \frac{\partial a_\ell}{\partial t}-\frac{\partial}{\partial \theta}\left(Ub_\ell-Av_\ell\right)=0, \end{gather}

(2.18)\begin{gather} \frac{\partial b_\ell}{\partial t}+\frac{1}{\sin\theta}\frac{\partial}{\partial \phi}(Ub_\ell-Av_\ell)=0. \end{gather}

For mathematical convenience, we then introduce the notation

(2.19a,b)\begin{equation} U(\theta)=\varOmega(\theta) \sin\theta,\quad A(\theta) =\beta(\theta) \sin\theta. \end{equation}

With the radius of the sphere as unity in our non-dimensional system, $\varOmega$ is the angular velocity and $\beta$ is the magnetic analogue as noted by Gilman & Fox (1997). For a flow and magnetic field that are smooth at the poles $\theta =0$, ${\rm \pi}$, the quantities $\varOmega$ and $\beta$ will tend to constants there. In the absence of a better term, we will refer to $\beta$, inaccurately, as the magnetic field from now on. In view of the divergence-free conditions (2.13) and (2.14), we introduce the streamfunction $\psi$ and the flux function $\chi$, such that

(2.20a–d)\begin{equation} u_\ell={-}\frac{\partial \psi}{\partial \theta} ,\quad v_\ell=\frac{1}{\sin\theta}\frac{\partial \psi}{\partial \phi} ,\quad a_\ell={-}\frac{\partial\chi}{\partial \theta} ,\quad b_\ell=\frac{1}{\sin\theta}\frac{\partial\chi}{\partial\phi}. \end{equation}

We combine (2.15) and (2.16) to eliminate the pressure $p$, and then we apply (2.20a–d). After some algebra, following Watson (1981), we derive the vorticity equation

(2.21)$$\begin{gather} \left(\frac{\partial}{\partial t}+\varOmega\frac{\partial}{\partial \phi}\right){\nabla}^2\psi-\frac{1}{\sin \theta}\frac{\mathrm{d}}{\mathrm{d}\theta}\left[\frac{1}{\sin\theta}\frac{\mathrm{d}}{\mathrm{d}\theta}(\varOmega\sin^2\theta)\right]\frac{\partial \psi}{\partial \phi} \nonumber\\ -\beta\frac{\partial}{\partial \phi}{\nabla}^2\chi+\frac{1}{\sin \theta}\frac{\mathrm{d}}{\mathrm{d}\theta}\left[\frac{1}{\sin\theta}\frac{\mathrm{d}}{\mathrm{d}\theta}(\beta\sin^2\theta)\right]\frac{\partial \chi}{\partial \phi}=0, \end{gather}$$

where the Laplacian on the spherical surface is

(2.22)\begin{equation} {\nabla}^2 f \equiv\frac{1}{\sin\theta}\frac{\partial }{\partial\theta}\left(\sin\theta \frac{\partial f}{\partial\theta}\right)+\frac{1}{\sin^2\theta}\frac{\partial^2 f}{\partial\phi^2}. \end{equation}

Similarly, using (2.20a–d) in the induction equation, (2.17) and (2.18), we obtain

(2.23)\begin{equation} \left(\frac{\partial}{\partial t}+\varOmega \frac{\partial}{\partial \phi}\right)\chi-\beta \frac{\partial \psi}{\partial \phi}=0. \end{equation}

Now we consider normal mode disturbances, replacing

(2.24)\begin{equation} (\psi, \chi)\rightarrow(\psi, \chi) \exp[{\mathrm{i}m(\phi-ct)}], \end{equation}

where $m$ is an integer representing the wavenumber in the longitudinal direction and $c$ is a complex constant representing the phase velocity. Using the substitution $\mu =\cos \theta$, (2.21) and (2.23) become ordinary differential equations in $\mu$,

(2.25)\begin{gather} (\varOmega-c)L\psi-\psi \frac{\mathrm{d}^2}{\mathrm{d}\mu^2}[\varOmega(1-\mu^2)]-\beta L\chi+\chi \frac{\mathrm{d}^2}{\mathrm{d}\mu^2}[\beta(1-\mu^2)]=0,\end{gather}

(2.26)\begin{gather} (\varOmega-c)\chi-\beta\psi=0, \end{gather}

where $L$ is the Legendre operator

(2.27)\begin{equation} Lf \equiv \frac{\mathrm{d}}{\mathrm{d}\mu}\left[(1-\mu^2) \frac{\mathrm{d}f}{\mathrm{d}\mu}\right]-\frac{m^2 f}{1-\mu^2} . \end{equation}

We note that (2.25) and (2.26) do not hold for axisymmetric disturbances, having $m=0$ (and indeed are written down having been divided throughout by $\mathrm {i}m$). The continuity equation (2.13) does not allow axisymmetric disturbances that remain finite at the poles and such disturbances are only possible when a free surface is present; see Gilman & Dikpati (2002). Without loss of generality, we take $m$ to be a positive integer.

Finally, for the subsequent analysis, it is convenient to introduce the variable $H$ defined by

(2.28)\begin{equation} H=\frac{\psi}{\varOmega-c}=\frac{\chi}{\beta} , \end{equation}

motivated by (2.26). Its governing equation is

(2.29)\begin{equation} (SH')'+\left[2(\varOmega-c)(\mu\varOmega)'-2\beta(\mu\beta)'-\frac{m^2S}{(1-\mu^2)^2}\right]H=0, \end{equation}

with

(2.30)\begin{equation} S=[(\varOmega-c)^2-\beta^2](1-\mu^2), \end{equation}

where the prime denotes a derivative with respect to $\mu$. The physical meaning of $H$ is that it is the Lagrangian displacement in the $\theta$ direction scaled by $\sin \theta$.

The boundary conditions for (2.25), (2.26) and (2.29) are that $\psi$, $\chi$ and $H$ must vanish at the poles $\mu =\pm 1$, these locations being singularities of the equations. This poses an eigenvalue problem for the phase velocity $c$. Our particular interest is the situation where we have a complex eigenvalue $c=c_{r}+\mathrm {i}c_{i}$ with $c_{i}>0$, indicating the presence of instability. In any case since (2.29) is real, solutions for $c$ always appear in complex conjugate pairs.

Equation (2.29) has singularities at locations $\mu _\star$ where $S=0$ or

(2.31)\begin{equation} \varOmega-c={\pm} \beta. \end{equation}

In general, $H$ diverges at such points if the eigenvalue $c$ is real, and $\mu _\star$ is referred to as a critical level. In the context of instability, i.e. when $c=c_{r}+\mathrm {i}c_{i}$, $c_{i}>0$, the solution remains analytic (for real $\mu$) but there are strong gradients of the eigenfunctions in the critical layer where $\varOmega -c_{r}\approx \pm \beta$, since $c_{i}$ is usually small. We will show that the critical layer plays a fundamental role in the instability problem.

A main contribution of this paper is an asymptotic analysis of the eigenvalue problem for the ‘clamshell instability’, presented in § 4 in detail. In the limit of weak shear of the zonal flow $\varOmega '$, we derive an asymptotic solution for the eigenvalue and eigenfunction using the method of matched asymptotic expansions, combining solutions in the bulk of the flow and the critical layer. We will also solve the problem numerically by adopting two methods from previous studies: one expands $\psi$ and $\chi$ using Legendre polynomials (Gilman & Fox 1999), and the other is a shooting method (Dikpati & Gilman 1999). The method of Legendre polynomial expansion computes all of the eigenvalues but it can be expensive. Similarly to Gilman & Fox (1999), we use this method when $\varOmega (\mu )$ and $\beta (\mu )$ are expressed by polynomials, so that their expansions merely involve several terms, resulting in a sparse matrix for the eigenvalue problem. The shooting method is fast and provides more precise solutions, but it needs a good guess for the eigenvalue. Such a good guess will either come from the method of Legendre polynomial expansion or the asymptotic solution derived in the limit of weak shear. When these two approaches are not available, some trial and error for the initial guess has to be performed.

In what follows, we will also be interested in the mean-flow response of the linear instability. We therefore set

(2.32a–d)\begin{equation} u=U+u_\ell+\Delta U,\quad v=v_\ell+\Delta V,\quad a=A+a_\ell+\Delta A,\quad b=b_\ell+\Delta B, \end{equation}

where the mean-flow modifications $\Delta U$, $\Delta V$, $\Delta A$ and $\Delta B$ are forced by the linear disturbances and are independent of $\phi$. Substituting (2.32a–d) into (2.5)–(2.10), applying the zonal average

(2.33)\begin{equation} \bar{f}(\theta) =\frac{1}{2{\rm \pi}}\int_0^{2{\rm \pi}} f (\theta, \phi)\, \mathrm{d}\phi, \end{equation}

and noting the spatial periodicity of disturbances in $\phi$, we find that $\Delta V$ and $\Delta B$ are zero and $\Delta U$ and $\Delta A$ are governed by

(2.34)\begin{gather} \frac{\partial \Delta U}{\partial t}=\frac{1}{\sin^2 \theta}\frac{\partial}{\partial \theta} \overline{\sin^2\theta (a_\ell b_\ell-u_\ell v_\ell)},\end{gather}

(2.35)\begin{gather} \frac{\partial\Delta A}{\partial t}=\frac{\partial}{\partial \theta} \overline{(u_\ell b_\ell -v_\ell a_\ell)} . \end{gather}

In the spherical system, angular momentum and the mean toroidal field are conserved, namely

(2.36)\begin{gather} \int_0^{\rm \pi} 2{\rm \pi} \sin^2\theta \frac{\partial \Delta U}{\partial t}\,\mathrm{d}\theta=0, \end{gather}

(2.37)\begin{gather} \int_0^{\rm \pi} \frac{\partial\Delta A}{\partial t}\,\mathrm{d}\theta=0, \end{gather}

as is guaranteed by (2.34) and (2.35).

3. The semicircle rules

In this section we derive semicircle rules that provide general bounds for the eigenvalue, in the style of the celebrated theory of Howard (1961). To tackle the spherical geometry, we will mainly follow the theory of Watson (1981) who derived semicircle rules for hydrodynamic instability on a sphere. We will also provide alternative bounds to his theory, which could be tighter for certain types of flows.

We proceed by multiplying (2.29) by the complex conjugate $H^*$ of $H$ and integrating from $\mu =-1$ to $\mu =1$. Applying integration by parts, we derive the integral formula

(3.1)$$\begin{gather} \int_{{-}1}^1\left[(\varOmega-c)^2-\beta^2\right]\left[(1-\mu^2)|H'|^2+\frac{m^2}{1-\mu^2}|H|^2\right]\,\mathrm{d}\mu \nonumber\\ =\int_{{-}1}^1 2\left[(\varOmega-c)(\mu\varOmega)'-\beta(\mu\beta )'\right]|H|^2\, \mathrm{d}\mu. \end{gather}$$

Writing the complex phase velocity as $c=c_{r}+\mathrm {i}c_{i}$ and taking $c_{i}>0$ for an unstable mode, the imaginary part of (3.1) yields

(3.2)\begin{equation} \int_{{-}1}^1\varOmega G\, \mathrm{d}\mu=\int_{{-}1}^1 [c_{r} G+(\mu\varOmega)'|H|^2\,]\, \mathrm{d}\mu, \end{equation}

with

(3.3)\begin{equation} G=(1-\mu^2)\left|H'\right|^2+\frac{m^2}{1-\mu^2}\, |H|^2\ge 0. \end{equation}

This result has been derived by Gilman & Fox (1997). The real part of (3.1) is

(3.4)\begin{equation} \int_{{-}1}^1 (\varOmega^2-2\varOmega c_{r}+c_{r}^2-c_{i}^2-\beta^2)G\, \mathrm{d}\mu=\int_{{-}1}^1 2 \left[(\varOmega-c_{r}) (\mu\varOmega)'-\beta(\mu\beta)'\right]|H|^2 \, \mathrm{d}\mu. \end{equation}

Using (3.2) to replace the second term on the left-hand side of (3.4), we have

(3.5)\begin{equation} (c_{r}^2+c_{i}^2)\int_{{-}1}^1G\, \mathrm{d}\mu=\int_{{-}1}^1(\varOmega^2-\beta^2)G\, \mathrm{d}\mu+ \int_{{-}1}^1 2\left[\beta(\mu\beta)'-\varOmega(\mu\varOmega)'\right]|H|^2\, \mathrm{d}\mu. \end{equation}

Now, following Howard (1961) and Watson (1981), we quote the inequality

(3.6)\begin{equation} \int_{{-}1}^1(\varOmega-\varOmega_{max})(\varOmega-\varOmega_{min})G\, \mathrm{d}\mu\leq 0, \end{equation}

where $\varOmega _{max}$ and $\varOmega _{min}$ are the maximum and minimum values of $\varOmega$ for all $\mu$, which gives

(3.7)\begin{equation} \int_{{-}1}^1\varOmega^2 G\, \mathrm{d}\mu\leq \int_{{-}1}^1\left[\varOmega(\varOmega_{\max}+\varOmega_{min})-\varOmega_{max}\varOmega_{min}\right] G\, \mathrm{d}\mu. \end{equation}

Substituting (3.7) into (3.5) to replace the $\varOmega ^2G$ term, we derive

(3.8)$$\begin{gather} (c_{r}^2+c_{i}^2)\int_{{-}1}^1G\, \mathrm{d}\mu\leq \int_{{-}1}^1 [(\varOmega_{max}+\varOmega_{min})\varOmega-\varOmega_{\max}\varOmega_{min}-\beta^2 ]G\, \mathrm{d}\mu \nonumber\\ + \int_{{-}1}^1 2 [\beta(\mu\beta)'-\varOmega(\mu\varOmega)']|H|^2\,\mathrm{d}\mu. \end{gather}$$

Finally, applying (3.2)–(3.8) again leads to the inequality

(3.9)$$\begin{gather} [(c_{r}-\bar{\varOmega})^2+c_{i}^2]\int_{{-}1}^1G\, \mathrm{d}\mu \leq\int_{{-}1}^1\left(\Delta \varOmega^2-\beta^2\right)G\, \mathrm{d}\mu \nonumber\\ +\int_{{-}1}^1 2 \left[ \beta(\mu\beta)'-(\varOmega-\bar{\varOmega})(\mu \varOmega)' \right]|H|^2\, \mathrm{d}\mu, \end{gather}$$

where

(3.10a,b)\begin{equation} \bar{\varOmega}=\frac{\varOmega_{max}+\varOmega_{min}}{2} ,\quad \Delta \varOmega=\frac{\varOmega_{max}-\varOmega_{min}}{2}. \end{equation}

We have obtained two relations for $c_{r}^2+c_{i}^2$ and $(c_{r}-\bar {\varOmega })^2+c_{i}^2$, namely (3.5) and (3.9). In the case of MHD instability in Cartesian geometry, e.g. the study of Hughes & Tobias (2001), two similar equations hold but the integrals involving $|H|^2$ are not present. In that case, one can derive two semicircle rules straightforwardly by bounding the integrals of $G$. In our case of spherical geometry, however, we need to consider how to bound the two $|H|^2$ integrals in order to find semicircle rules. First, we note that we have one bound for $|H|^2$ from the definition of $G$ in (3.3), namely

(3.11)\begin{equation} 0\le |H|^2\le \frac{1-\mu^2}{m^2} G. \end{equation}

This holds at each location of $\mu$, and so we refer to it as the pointwise bound. An alternative is to bound the integral of $|H|^2$. For this task, we invoke the theorem of Rayleigh's quotient. Let $L$ be a linear Sturm–Liouville operator and $\lambda$ be the smallest eigenvalue for the corresponding Sturm–Liouville problem $L H+\lambda H=0$, with homogeneous boundary conditions at $\mu =a$ and $\mu =b$. Then, for arbitrary smooth functions $H(\mu )$, the Rayleigh quotient $R$ satisfies

(3.12)\begin{equation} R={-} \frac{\displaystyle\int_{a}^b H^*L H\, \mathrm{d}\mu}{\displaystyle\int_a^b|H|^2\, \mathrm{d}\mu}\ge \lambda. \end{equation}

The smallest eigenvalue of the Legendre operator $L$ defined in (2.27) is $\lambda =m(m+1)$. Thus, applying integration by parts to the numerator of (3.12), we obtain

(3.13)\begin{equation} 0\le \int_{{-}1}^1 |H|^2\, \mathrm{d}\mu\le \frac{1}{m(m+1)}\int_{{-}1}^1 G\, \mathrm{d}\mu. \end{equation}

We refer to (3.13) as the integral bound for $|H|^2$.

Now we apply our two bounds (3.11) and (3.13) to (3.5) and (3.9). For the pointwise bound, substituting (3.11) into (3.5) gives

(3.14)\begin{align} (c_{r}^2+c_{i}^2)\int_{{-}1}^1G\, \mathrm{d}\mu&\le \int_{{-}1}^1\left\{\varOmega^2-\beta^2+\frac{1-\mu^2}{m^2} 2\left[\beta(\mu\beta)'-\varOmega(\mu\varOmega)'\right]^+\right\}G\, \mathrm{d}\mu, \nonumber\\ &\le \left\{\varOmega^2-\beta^2+\frac{1-\mu^2}{m^2} 2\left[\beta(\mu\beta)'-\varOmega(\mu\varOmega)'\right]^+\right\}_{max}\int_{{-}1}^1G\, \mathrm{d}\mu, \end{align}

where the plus sign superscript is defined by

(3.15)\begin{equation} f^+=\max(f,0). \end{equation}

We have therefore arrived at the semicircle rule

(3.16)\begin{equation} c_{r}^2+c_{i}^2\le \left\{\varOmega^2-\beta^2+\frac{1-\mu^2}{m^2} 2\left[\beta(\mu\beta)'-\varOmega(\mu\varOmega)'\right]^+\right\}_{max}. \end{equation}

Similarly, if we apply (3.11) to (3.9), we derive another semicircle rule using the pointwise bound

(3.17)\begin{equation} (c_{r}-\bar{\varOmega})^2+c_{i}^2\le\left\{\Delta\varOmega^2-\beta^2+\frac{1-\mu^2}{m^2} 2\left[\beta(\mu\beta)'-(\varOmega-\bar{\varOmega})(\mu\varOmega)'\right]^+\right\}_{max}. \end{equation}

To apply instead the integral bound, we first take the functions multiplying $G$ and $|H|^2$ out of the integral in (3.5) and then use (3.13). We obtain

(3.18)\begin{align} (c_{r}^2+c_{i}^2)\int_{{-}1}^1G\, \mathrm{d}\mu &\le (\varOmega^2-\beta^2)_{max}\int_{{-}1}^1G\, \mathrm{d}\mu+2\left[\beta(\mu\beta)'-\varOmega(\mu\varOmega)'\right]_{max}\int_{{-}1}^1|H|^2\, \mathrm{d}\mu \nonumber\\ &\le (\varOmega^2-\beta^2)_{max}\int_{{-}1}^1G\, \mathrm{d}\mu+\frac{2\left[\beta(\mu\beta)'-\varOmega(\mu\varOmega)'\right]_{max}^+}{m(m+1)}\int_{{-}1}^1G\, \mathrm{d}\mu. \end{align}

Hence, we have another semicircle rule

(3.19)\begin{equation} c_{r}^2+c_{i}^2\le (\varOmega^2-\beta^2)_{max}+\frac{2}{m(m+1)} \left[\beta(\mu\beta)'-\varOmega(\mu\varOmega)'\right]_{max}^+. \end{equation}

Similarly, applying the integral bound to (3.9) yields

(3.20)\begin{align} (c_{r}-\bar{\varOmega})^2+c_{i}^2\le (\Delta\varOmega^2-\beta^2)_{max}+\frac{2}{m(m+1)} \left[\beta(\mu\beta)'-\left(\varOmega-\bar{\varOmega}\right)(\mu \varOmega)'\right]_{max}^+. \end{align}

In summary, we have derived four semicircle rules: (3.16), (3.17), (3.19) and (3.20). Watson's (1981) semicircle rule corresponds to (3.17) with the field switched off (though the latter is tighter due to a more careful treatment of the geometric term). In the limit of small scale, $m\rightarrow \infty$, the geometric terms proportional to $(1-\mu ^2)/m^2$ or ${1}/{m(m+1)}$ vanish and we recover the semicircle rules derived by Gilman (1967), Cally (2000) and Hughes & Tobias (2001) for MHD instability in Cartesian geometry, which can be further reduced to the theory of Howard & Gupta (1962) when the magnetic field is uniform. When both limits are applied, (3.17) and (3.20) both reduce to the semicircle rule of Howard (1961). In that problem, the semicircle centred at $(c_{r},c_{i})=(0,0)$ completely includes the other centred at $(\bar {\varOmega },0)$ and becomes redundant.

A distinguishing feature of the current theory is that, for each possible semicircle centre, i.e. $(c_{r}, c_{i})=(0,0)$ and $(\bar {\varOmega },0)$, there are two possible radii resulting from the two different methods used to bound $|H|^2$. The smaller radius will represent a tighter bound and, thus, be the effective one. We have found that, depending on the profiles of $\varOmega$ and $\beta$, either of the two bounding methods can be tighter. In general, the results of the pointwise bounding, (3.16) and (3.17), are tighter when the shear of $\varOmega$ or $\beta$ are prominent, because these bounds take the maximum of the sum of functions, in contrast to (3.19) and (3.20) that take the sum of the maxima of two functions. On the other hand, when the shears of $\varOmega$ and $\beta$ are weak, the integral bounds (3.19) and (3.20) can be tighter due to the smaller coefficient of ${1}/{m(m+1)}$.

In an effort to make the results more compact and uniform, we present an alternative way to write the semicircle rules using functional expressions,

(3.21a,b)\begin{equation} c_{r}^2+c_{i}^2\le E[ f_1,g_1],\quad (c_{r}-\bar{\varOmega})^2+c_{i}^2\le E[ f_2,g_2], \end{equation}

where $E$ is a functional defined by

(3.22)\begin{equation} E\left[f,g\right]=\min\left(C[f,g],D[f,g]\right), \end{equation}

with

(3.23a,b)\begin{equation} C[f,g]=\left\{f+\frac{1-\mu^2}{m^2}\, g^+\right\}_{max},\quad D[f,g]=f_{max}+\frac{1}{m(m+1)} g_{max}^+, \end{equation}

corresponding to the pointwise or integral bounds, respectively. The functions in (3.21) are

(3.24a,b)\begin{gather} f_1=\varOmega^2-\beta^2,\quad g_1=2\left[\beta(\mu\beta)'-\varOmega(\mu\varOmega)'\right], \end{gather}

(3.25a,b)\begin{gather} f_2=\Delta \varOmega^2-\beta^2,\quad g_2=2\left[\beta(\mu\beta)'-(\varOmega-\bar{\varOmega})(\mu\varOmega)'\right]. \end{gather}

We now demonstrate the application of these semicircle rules to instability problems. Following Hughes & Tobias (2001), we first study the instability criterion. Since all the rules need to be satisfied by the complex value of $c$ for an arbitrary unstable mode, if the radius of either semicircle disappears, or if the semicircles for the two centres become disjoint, then instability is impossible. Using the notation in (3.21), these conditions are

(3.26a–c)\begin{equation} E[f_1,g_1] \le0 \quad \mathrm{or}\quad E[f_2,g_2] \le0 \quad \mathrm{or} \quad \sqrt{E[f_1,g_1]}+\sqrt{E[f_2,g_2]}\le|\bar{\varOmega}|, \end{equation}

any of which is a sufficient condition for stability. A straightforward example of an application of (3.26b) is the case of constant $\varOmega$ and $\beta$: the semicircle rule (3.20) for this flow becomes

(3.27)\begin{equation} (c_{r}-\bar{\varOmega})^2+c_{i}^2\le \left[\frac{2}{m(m+1)}-1\right]\beta^2\le 0, \end{equation}

given that $m\ge 1$. Hence, the flow is linearly stable when $\varOmega$ and $\beta$ are constants. In fact, the phase velocities for this flow can be solved analytically, since the waves are spherical harmonics (Márquez-Artavia, Jones & Tobias 2017), but through this example we have shown the power of the semicircle rule: for this specific flow, it can give a bound that is tight enough to exclude the possibility of instability. We note that this analysis is based on (3.20) which is found using the integral bound (itself tight for spherical harmonics); use of the pointwise bound is not tight enough to rule out instability for this flow.

We recall that, for the MHD problem in Cartesian geometry, for example, Hughes & Tobias (2001), the magnetic field may only reduce the radii of the semicircles, or keep them the same. Hence, they explored how the field profile may realise one of (3.26) to guarantee stability, and studied the tightness of these bounds. In our spherical problem, on the other hand, the important feature is that the magnetic field may increase the radii of the semicircles through the geometric terms (those proportional to $(1-\mu ^2)/m^2$ or ${1}/{m(m+1)}$), and this is often accompanied by the destabilising effect of the magnetic field. It is therefore of interest to show examples of the semicircles compared with the actual eigenvalues for sample profiles. Motivated by the solar differential rotation profile, for the basic state angular velocity, we consider the typical differential shear given by (1.1). The parameters are chosen as $r=1$ and $s=0.24$ as in Gilman & Fox (1997). For the basic state magnetic field, we consider two examples: the first is a ‘linear shear’ profile considered by Gilman & Fox (1997) and Cally (2001),

(3.28)\begin{equation} \beta=\sigma\mu, \end{equation}

where $\sigma$ is a constant. The second is a profile with a pair of opposite Gaussian distributions, studied by Dikpati & Gilman (1999) and Cally (2001),

(3.29)\begin{equation} \beta=\frac{\sigma}{2\sqrt{1-d^2}}\left[\exp\left( -\frac{4 (\mu-d)^2}{w^2(1-d^2)}\right)-\exp\left( -\frac{4 (\mu+d)^2}{w^2(1-d^2)}\right)\right], \end{equation}

where $\pm d$ are the centres of two Gaussian distributions and $w$ is a width parameter. These profiles are idealised models for the solar magnetic field. From solar magnetogram observations, the magnetic field is antisymmetric about the equator, and (3.28) is the simplest antisymmetric profile while (3.29) further models the belt patterns seen from observations. We plot the semicircles for the field profiles (3.28) and (3.29) in figure 1(a,b) in solid lines, respectively. Mode $m=1$ is chosen, which is the only wavenumber where the flow is unstable. We choose $\sigma =1$ for (3.28) corresponding to a strong field, and $\sigma =0.1$ for (3.29) corresponding to a weak field, and take $w = {\rm \pi}/6$ and $d=1/\sqrt {2}$ for (3.29). Due to the strong shear in $\beta$, the semicircles (3.16) and (3.17) from use of the pointwise bound (3.22a) have smaller radii, and so are plotted. We also plot the semicircle (3.17) in the absence of a magnetic field in dashed lines; the other semicircle (3.16) is the same with or without the magnetic field. There is an unstable mode for each field profile, which is solved numerically and we plot as a star in each panel. As mentioned earlier, in the absence of a magnetic field, the flow (1.1) is unstable when $s/r>0.29$ (Watson 1981); hence, the flow with $r=1$, $s=0.24$ that is considered here is hydrodynamically stable, and the instability is induced by the magnetic field.

Figure 1. Semicircle rules for the zonal flow (1.1) with $r=1$, $s=0.24$ and the magnetic field of (a) profile (3.28) with $\sigma =1$ and (b) profile (3.29) with $\sigma =0.1$, $d=1/\sqrt {2}$ and $w = {\rm \pi}/6$. The wavenumber is $m=1$ for both figures. Solid lines represent the result of (3.16) and (3.17), and dashed lines represent the result of (3.17) without the magnetic field. All semicircles plotted are the tightest, here from the pointwise bound. The star represents the numerical solution for the unstable mode: $c=0.95+0.023\mathrm {i}$ for panel (a) and $c=0.87+0.0074\mathrm {i}$ for panel (b).

The region where the two semicircles with solid lines overlap is the domain for all possible values of $c$. The star lies inside this region, as we expect. For the strong field case of figure 1(a), the magnetic field significantly enlarges the original, hydrodynamic semicircle (shown dashed) centred at $\varOmega =\bar {\varOmega }$, though the actual unstable mode still lies inside this semicircle. For the weak field case of figure 1(b), the magnetic field slightly reduces the purely hydrodynamic semicircle, which is a little surprising since the field has a destabilising effect (the $\beta =0$ system is stable for this flow). In general, the situation of figure 1(a) (in which the magnetic field enlarges the semicircle centred at $\varOmega =\bar {\varOmega }$) is typical for field-induced instabilities. The opposite case (in which the field destabilises the flow but reduces the radius, as in figure 1b) is relatively rare: it only happens for a weak field with certain profiles. As suggested by Gilman & Fox (1997), the maximum possible growth rate from the semicircle rule can be much larger than the actual growth rate, but the rules are still powerful in giving rigorous bounds on $c$ in the complex plane. Note that in figure 1(a) the star is very close to the edge of the semicircle centred at $\varOmega =0$, suggesting a tight bound in this case.

Further remarks on the situation seen in figure 1(a) may be helpful here. We note that the magnetic field increases one of the semicircles, but the eigenvalue still lies inside the original hydrodynamic semicircle. The explanation is that the bounding of $H$ in terms of $|G|^2$ is over all possible functions, not necessarily the solutions of the eigenvalue problem; thus, the bound can give a much larger domain than that attained by an actual eigenvalue. We note that we have not found unstable modes outside the hydrodynamic semicircle, and so we cannot answer the question of whether the larger semicircle induced by the magnetic field really represents a larger possible domain, or just a looser bound. We leave this issue for future consideration.

We give another example of the prediction of the semicircle rules for profiles that are related to the analysis of the clamshell instability in the subsequent sections. We consider a magnetic field with $\beta =0$ at a certain latitude, which coincides with the location where $|\varOmega |$ reaches its maximum. The standard profiles (1.1) with (3.28) or (3.29) studied above have this property when $rs>0$ (the case for solar differential rotation). It can be shown from (3.21)–(3.25a,b) that

(3.30a,b)\begin{equation} E[f_1, g_1] \ge ({\varOmega^2})_{max},\quad E[f_2,g_2] \ge \Delta\varOmega^2. \end{equation}

If $\varOmega$ is not a constant function of latitude, in other words the fluid flow is sheared, then using (3.30a,b) it is evident that the sufficient conditions for stability given in (3.26) are never satisfied. In fact, the lower limits given by (3.30a,b) are the results of Howard (1961), where there is a finite-area semicircle for any sheared flow in the plane. A key observation here is that the semicircle rules we have obtained always allow the possibility of instability in a sheared flow on the sphere provided the magnetic field vanishes at the latitude where the rotation rate is greatest. We will demonstrate that such an instability does indeed exist, in the limit of weak shear, by our analysis in § 5.1.

We finally note that in addition to the theory we present, there could be other versions of semicircle rules. Our derivation is based on the equation for $H=\psi /(\varOmega -c)$ in Sturm–Liouville form, following the approach of Watson (1981), but Thuburn & Haynes (1996) and Sasaki et al. (2012) worked on the Sturm–Liouville equation of $\eta =\psi /[(\varOmega -c)\sqrt {1-\mu ^2}]$ and derived semicircle rules different from Watson's. At present, we do not know which version is tighter. Indeed, there could more versions that result from other substitutions. For MHD instability in Cartesian geometry, Deguchi (2021) improved the traditional semicircle rules by finding the inner envelope of a family of semicircles, and so a yet tighter bound. These approaches could provide avenues to extend our theory, and we leave them for future study.

4. Analysis of the clamshell instability

Cally (2001) and Cally et al. (2003) have shown that when the magnetic field is strong and its profile is broad, MHD instability on a sphere features a clamshell pattern in which the field lines on the two hemispheres are tilted in opposite directions. In figure 2 we give an example of such a clamshell pattern, which corresponds to the basic state (1.1) and (3.28) plus the unstable mode with its eigenfunction shown in figure 3. Without loss of generality, we normalise the eigenfunction by $\max |H|=0.07$ and $\mathrm {Im}\, H=0$ as $\mu \rightarrow 1$ throughout the paper; other normalisations would yield a similar pattern to figure 2 as long as $|H|$ remains small.

Figure 2. Magnetic field lines for the instability of the basic state $\varOmega =1-0.1\mu ^2$, $\beta =\mu$ with wavenumber $m=1$. The eigenvalue is obtained numerically as $c=0.982+0.00847\mathrm {i}$. The plotted results correspond to the basic state plus the unstable mode with a small amplitude shown in figure 3.

Figure 3. Eigenfunctions for the instability of the profile given in figure 2. The tilting mode (4.5a–d) is plotted using circles. The eigenfunctions are normalised by $\max |H|=0.07$ and $\mathrm {Im} H=0$ as $\mu \rightarrow 1$.

In this section we will undertake an asymptotic analysis of this instability in the limit of weak shear $\varOmega '$ of the basic zonal flow $\varOmega$. We will derive an asymptotic solution for the eigenvalue $c$ in this limit, which will provide insights into the instability mechanism applying to a wide range of profiles of $\varOmega$ and $\beta$.

4.1. The tilting mode

While (2.25), (2.26) and (2.29) look quite complicated, for azimuthal wavenumber $m=1$, they admit a very simple solution for arbitrary $\varOmega$ and $\beta$, namely

(4.1a–d)\begin{equation} c=0,\quad \psi=\sqrt{1-\mu^2}\, \varOmega,\quad \chi=\sqrt{1-\mu^2}\, \beta,\quad H=\sqrt{1-\mu^2} . \end{equation}

The physical meaning of this solution lies in the rotational invariance of the spherical geometry: if we apply a solid body rotation to the entire zonal flow and magnetic field, the tilted flow and field are still solutions of the governing equations. When the angle of tilting is small, the solution given in (4.1a–d) gives the difference between the tilted and original states. To see this, first note that the rotation through an infinitesimal angle $\alpha$ about the $x$ axis is given by $(x,y,z) \to (x,y - \alpha z, z+\alpha y)$ and corresponds to the spherical polar coordinate change

(4.2a,b)\begin{equation} \theta \to \theta - \alpha \sin \phi, \quad \phi \to \phi - \alpha \cot \theta \cos \phi. \end{equation}

Thus, any function $f(\theta, \phi )$ is mapped under such a rotation by

(4.3)\begin{equation} f(\theta, \phi) \to f(\theta + \alpha \sin \phi, \phi + \alpha \cot \theta \cos \phi) \simeq f (\theta,\phi) + \alpha [ \sin \phi \, f_\theta + \cot\theta \cos \phi \, f_\phi] . \end{equation}

Apply this to a streamfunction $\varPsi (\theta )$ for the basic state, satisfying $U = \varOmega \sin \theta = - \varPsi _\theta$ (cf. (2.20a–d)), and we find that the tilted streamfunction is

(4.4)\begin{align} \varPsi \to \varPsi + \alpha \sin \phi \, \varPsi_\theta = \varPsi - \alpha \sin \phi \sin \theta\, \varOmega = \varPsi + \left(\frac{1}{2} \mathrm{i} \alpha \sqrt{1-\mu^2}\, \varOmega \, {\rm e}^{\mathrm{i} \phi} + \mathrm{c.c.}\right). \end{align}

The difference corresponds to a multiple $\tfrac {1}{2}\mathrm {i}\alpha$ of the steady solution (4.1a–d) to the linear problem (and likewise a multiple $\tfrac {1}{2}\alpha$ is a rotation around the $y$ axis). For hydrodynamic stability, this neutral mode was noticed by Watson (1981) (although there appears to be a mistake in the form of the eigenfunction given). A well-known analogue is the neutral mode arising from the translational invariance of a vortex in the plane: shifting the entire $m=0$ vortex gives an $m=1$ solution to the linear problem with zero eigenvalue, as noted by Bernoff & Lingevitch (1994).

Interestingly, for the MHD problem, (2.25), (2.26) and (2.29) also admit a slightly different mode that involves the tilting of the basic state. If $\varOmega$ is a constant, i.e. the zonal flow is solid body rotation, then for arbitrary $\beta$, we also have an exact solution for $m=1$,

(4.5a–d)\begin{equation} c=\varOmega, \quad \psi=0, \quad \chi=\sqrt{1-\mu^2} \beta, \quad H=\sqrt{1-\mu^2} . \end{equation}

In this solution there is no velocity disturbance and the magnetic field lines are slightly tilted, similar to (4.1a–d), but now they rotate with the solid body rotation $\varOmega (\mu )=c$, frozen into the flow. Furthermore, if the magnetic field vanishes somewhere, say

(4.6)\begin{equation} \beta=0,\quad {\mathrm{at}} \ \mu=\mu_\star, \end{equation}

then (2.31) is satisfied and $\mu _\star$ is a critical level where (2.29) becomes singular. However, the solution (4.5a–d) remains regular there. The mathematical picture is that (2.29) has two linearly independent solutions: one is singular at $\mu =\mu _\star$, the other is regular there and is the solution (4.5a–d).

The solution (4.5a–d) is closely related to the clamshell instability. From now on, we refer to it as the tilting mode. For the example shown in figures 2 and 3, the shear of $\varOmega$ is weak, and so we expect that the tilting mode (4.5a–d), derived for constant $\varOmega$, is relevant. In figure 3 we plot the eigenfunctions that correspond to the clamshell pattern in figure 2 together with (4.5a–d). We see that in most of the region the tilting mode agrees very well with the unstable mode; however, interestingly, the latter is reversed across $\mu =0$. It turns out that the weak shear in $\varOmega$ excites the singular solution, which becomes dominant at the critical level; this is at $\mu _\star =0$ for the magnetic profile $\beta =\mu$. The singularity has a dramatic impact on the eigenfunction: it reverses the sign of the tilting mode, opening up the clamshell in figure 2, and importantly it destabilises the flow. Based on this intuition, we will undertake an asymptotic analysis to solve the eigenvalue problem in § 4.2. Readers who are not interested in the technical details of the matched asymptotic expansions may jump to the final result given by (5.1), noting that a star subscript represents the value of a function evaluated at $\mu =\mu _\star$.

It is noted that we have not yet found a normal mode instability that is induced by the mode (4.1a–d). This is perhaps because this mode is ‘too stable’: the eigenvalue $c=0$ remains unchanged, no matter how $\varOmega$ and $\beta$ vary. Note that when $\varOmega =0$, (4.1a–d) and (4.5a–d) become identical, and then our analysis indicates that the resonance between the two modes may induce algebraic growth (instead of exponential growth) of the magnetic field. But given $\varOmega =0$ is not very common in astrophysics, this is of limited interest, and we will not discuss it further.

4.2. The matched asymptotic expansions

We consider profiles of $\varOmega (\mu )$ and $\beta (\mu )$ such that the variation in $\varOmega$ is weak and $\beta$ is zero at a location denoted by $\mu =\mu _\star$. For simplicity, in the present study we assume that there is only one such $\mu _\star$ in the domain (although the case of multiple $\mu _\star$ may be studied following the same method). We also take $\beta$ to pass through zero at $\mu =\mu _\star$ with a gradient that is not small. Thus, close to $\mu _\star$, we have the approximation

(4.7a,b)\begin{equation} \beta\approx \beta_\star'(\mu-\mu_\star),\quad \mu\approx \mu_\star, \end{equation}

where $\beta _\star '$ is $\beta '$ evaluated at $\mu =\mu _\star$, and we assume that $\beta _\star '$ is of order unity (or larger). Apart from these prescriptions, $\varOmega (\mu )$ and $\beta (\mu )$ are general functions. We perform matched asymptotic expansions that combine the bulk of the flow, where the disturbance is mainly represented by the tilting mode, and the critical layer near $\mu =\mu _\star$, where the solution changes rapidly. The azimuthal wavenumber will be fixed at $m=1$, since this is the only wavenumber that admits the tilting mode.

We express $\varOmega$ and $c$ by

(4.8a,b)\begin{equation} \varOmega=\varOmega_0+\varepsilon \varOmega_1 (\mu),\quad c=\varOmega_0+\varepsilon c_1 + \cdots. \end{equation}

Here the leading-order rotation is $\varOmega _0\neq 0$, which is a constant representing the angular velocity of the solid body rotation, while $\varepsilon$ is a small number representing the amplitude of the weak shear and $\varOmega _1$ is an arbitrary function representing the shear profile. The eigenvalue $c$ is $\varOmega _0$ at leading order, following the tilting mode (4.5a–d). The weak shear $\varepsilon \varOmega _1$ induces a small correction $\varepsilon c_1$ to the eigenvalue, and our goal is to determine this. Also, if $\varOmega _0=0$ in (4.8a,b), but $\varOmega$ as a whole is weak compared with $\beta$, we can still follow a similar asymptotic analysis. This situation is perhaps less relevant for astrophysical applications and so the derivation is consigned to Appendix A; it is different in detail, but the final result (5.1) is the same.

4.2.1. Outer solution

Away from the critical level $\mu =\mu _\star$, we expand $H$ as

(4.9)\begin{equation} H=H_0(\mu)+\varepsilon H_1(\mu) + \cdots, \end{equation}

where $H_0$ has the profile of the tilting mode but is discontinuous across the critical layer, as we see in figure 3. Thus,

(4.10) \begin{equation} H_0(\mu)=\left\{\begin{array}{@{}ll} A_-\sqrt{1-\mu^2} , \quad \mu<\mu_\star, \\ A_+\sqrt{1-\mu^2} ,\quad \mu>\mu_\star, \end{array}\right. \end{equation}

where $A_+$ and $A_-$ are constants. Here $H_1$ represents the small correction that is induced by the weak shear in $\varOmega$. Substituting (4.8a,b) and (4.9) into (2.29) and collecting terms at $O(\varepsilon )$, we obtain the equation for $H_1$,

(4.11)\begin{equation} \left[\beta^2(1-\mu^2)H_1'\right]'+\left[2\beta(\mu\beta)'-\frac{\beta^2}{1-\mu^2}\right]H_1 =2\varOmega_0(\varOmega_1-c_1)H_0. \end{equation}

Using the method of variation of constants, noting that $H_1$ must be zero at $\mu =\pm 1$, the solution for $H_1$ can be obtained in the form

(4.12) \begin{align} H_1(\mu)=\left\{\begin{array}{@{}ll} h_{s}(\mu)\displaystyle \int_{ 1}^\mu2\varOmega_0\left[\varOmega_1(\nu)-c_1\right]H_0(\nu)h_{r}(\nu)\, \mathrm{d}\nu & \\ \quad -h_{r}(\mu)\displaystyle\int_1^\mu 2\varOmega_0[\varOmega_1(\nu)-c_1]H_0(\nu)h_{s}(\nu)\, \mathrm{d}\nu+C_a h_{r}(\mu), & \mu_\star<\mu<1,\\ h_{s}(\mu)\displaystyle\int_{- 1}^\mu2\varOmega_0\left[\varOmega_1(\nu)-c_1\right]H_0(\nu)h_{r}(\nu)\, \mathrm{d}\nu & \\ \quad -h_{r}(\mu)\displaystyle\int_{{-}1}^\mu 2\varOmega_0[\varOmega_1(\nu)-c_1]H_0(\nu)h_{s}(\nu)\, \mathrm{d}\nu+C_bh_{r}(\mu), & -1<\mu<\mu_\star. \end{array}\right. \end{align}

Here $h_{r}$ and $h_{s}$ are the regular and singular solutions of the corresponding homogeneous equation, given by

(4.13a,b) \begin{align} h_{r}(\mu)=\sqrt{1-\mu^2} ,\quad h_{s}(\mu)=\left\{\begin{array}{@{}ll} \sqrt{1-\mu^2}\displaystyle\int_{\mu_a}^\mu\dfrac{1}{\beta(\nu)^2(1-\nu^2)^2}\, \mathrm{d}\nu, & \mu_\star<\mu,\mu_a<1, \\ \sqrt{1-\mu^2}\displaystyle \int_{\mu_b}^\mu\dfrac{1}{\beta(\nu)^2(1-\nu^2)^2}\, \mathrm{d}\nu, & \; -1<\mu,\mu_b<\mu_\star,\end{array} \right. \end{align}

where $\mu _a$ and $\mu _b$ may be arbitrarily chosen in the given range, and $C_a$ and $C_b$ are undetermined constants, a consequence of $h_{r}$ satisfying both boundary conditions. Using (4.7a,b), we can deduce that $H_1, h_{s}\sim (\mu -\mu _\star )^{-1}$ as $\mu \rightarrow \mu _\star$. Hence, unlike the leading-order term $H_0$, the correction $H_1$ is essentially singular at $\mu =\mu _\star$. In other words, combining the tilting mode $H_0$ (in either hemisphere) with weak shear excites the singular mode of the system. When $\mu -\mu _\star =O(\varepsilon )$, $\varepsilon H_1$ becomes as large as $H_0$ and the expansion (4.9) becomes disordered. We note the singular behaviour of $H_1'$,

(4.14)$$\begin{gather} H_1'=\frac{A_\pm}{(1-\mu_\star^2)^{{3}/{2}}\beta_{{\star}}'^2(\mu-\mu_\star)^2}\int_{{\pm} 1}^{\mu_\star} \!2\varOmega_0\left[\varOmega_1(\mu)-c_1\right](1-\mu^2)\, \mathrm{d}\mu\nonumber\\ +O\left(\frac{1}{\mu-\mu_\star}\right),\quad \mu\rightarrow \mu_\star^{{\pm}}, \end{gather}$$

which we will use later.

4.2.2. Inner solution

The region $\mu -\mu _\star =O(\varepsilon )$ is the critical layer where $H$ varies significantly. We therefore introduce a local stretched coordinate, writing

(4.15a,b)\begin{equation} \eta=\frac{\mu-\mu_\star}{\varepsilon} ,\quad H=\mathcal{H}(\eta) + \cdots. \end{equation}

Substituting (4.7a,b), (4.8a,b) and (4.15a,b) into (2.29) we find that the leading-order local equation is

(4.16)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}\eta}\left\{\left[(\varOmega_{1\star}-c_1)^2-\beta_\star'^2\eta^2\right] \frac{\mathrm{d}\mathcal{H}}{\mathrm{d }\eta}\right\}=0, \end{equation}

where $\varOmega _{1\star }$ is $\varOmega _1$ evaluated at $\mu =\mu _\star$. Integrating with respect to $\eta$, we obtain

(4.17)\begin{equation} \frac{\mathrm{d} \mathcal{H}}{\mathrm{d}\eta}=\frac{\alpha_1}{(\varOmega_{1\star}-c_1)^2-\beta_\star'^2\eta^2 } , \end{equation}

and, hence,

(4.18)\begin{equation} \mathcal{H}=\frac{\alpha_1}{2(\varOmega_{1\star}-c_1)|\beta_\star'|}\, \left[\log\left(|\beta_\star'|\eta+\varOmega_{1\star}-c_1\right)- \log\left(|\beta'_\star|\eta-\varOmega_{1\star}+c_1\right)\right]+\alpha_2 \end{equation}

for some constants $\alpha _1$ and $\alpha _2$.

A striking property of the critical layer is the presence of a significant jump in $H$ from one side to the other, which is shown in figure 3(c) (also see the later figure 5). Such a jump can be understood from the asymptotic solution (4.18) as follows. If we select the branch cuts of the logarithm functions to lie on the negative real axis then, for $\mathrm {Im}\, c_1>0$, their large-variable limits are

(4.19a,b) \begin{equation} \log\left[|\beta_\star'|\eta\pm(\varOmega_{1\star}-c_1)\right]\sim\left\{\begin{array}{@{}ll} \log|\beta_\star'\eta|,& \eta\rightarrow +\infty, \\ \log|\beta_\star'\eta|\mp \mathrm{i}{\rm \pi},& \eta\rightarrow -\infty. \end{array}\right. \end{equation}

Hence,

(4.20a,b) \begin{equation} \mathcal{H}\sim\left\{\begin{array}{@{}ll} \alpha_2, & \eta\rightarrow+\infty, \\ \displaystyle \dfrac{-\mathrm{i}{\rm \pi}\alpha_1}{(\varOmega_{1\star}-c_1)|\beta_\star'|}+\alpha_2, & \eta\rightarrow-\infty. \end{array}\right. \end{equation}

The first term of (4.20b) clearly indicates the jump of $\mathcal {H}$ across the critical layer. This jump is important because it reverses the tilting direction to form the clamshell pattern, and induces instability through the presence of the imaginary unit $\mathrm {i}$, as we will see later. According to (4.19), the jump of the logarithm functions is contingent on the existence of a non-zero $|\beta _\star '|$, highlighting the role of the gradient of the magnetic field at the critical level.

The magnetic field lines in the critical layer, rendered by the local solution (4.18) plus the basic magnetic field, are shown in figure 4. It may be seen that the critical layer induces a pair of closed loops in the field line pattern, which is also visible in figure 2 near the equator. Note that such a pattern was not shown in the corresponding figure of Cally (2001) (first and second panel of his figure 4), since he did not draw field lines in the critical layer in the early stage of the evolution. Also note that at later times, Cally's (2001) simulation has shown that the field lines on the two hemispheres will reconnect as a result of dissipation. This is quite different from the ideal MHD instability that we currently study: in our figure 4, the field lines on the two hemispheres are separated.

Figure 4. The magnetic field lines in the critical layer, corresponding to figure 2 near the equator. We have superimposed the basic magnetic field with the asymptotic local solution (4.18), where $\alpha _1$ and $\alpha _2$ are found by matching to $A_-$ and $A_+$ via (4.22) and (4.24), and $A_-$ and $A_+$ are obtained by fitting to the numerical solution shown in figure 5.

Figure 5. The comparison of the eigenfunction $H$ between the asymptotic solution (solid lines) and numerical solution (circles) for the unstable mode of figure 3. The asymptotic solution of $H$ consists of the outer solution (4.9), (4.10), (4.12) and the inner solution (4.18), with the matching condition (4.22) and (4.24). The eigenvalue $c$ of the asymptotic solution (5.1) is $0.984+0.00778\mathrm {i}$, while its numerical solution is $0.982+0.00847\mathrm {i}$.

4.2.3. Matching and eigenvalue

Matching the inner and outer solution provides relations between the constants $\alpha _1$, $\alpha _2$, $A_-$ and $A_+$, and so determines the eigenvalue $c_1$. We first match $H'$ from the inner and outer solution in an intermediate region $\mu -\mu _\star =O(\varepsilon ^{{1}/{2}})$. Here $H_0'\ll \varepsilon H'$, so that we can neglect the former in the outer solution. Hence, matching (4.17) and (4.14) via

(4.21)\begin{equation} \frac{1}{\varepsilon}\frac{\mathrm{d}\mathcal{H}}{\mathrm{d}\eta} \left|_{\eta\rightarrow \pm \infty}= \varepsilon H_1' \right|_{\mu\rightarrow \mu_\star^\pm} , \end{equation}

we find that

(4.22a,b)\begin{align} \alpha_1 &=\frac{A_-}{(1-\mu_\star^2)^{{3}/{2}}}\int^{{-}1}_{\mu_\star}2\varOmega_0[\varOmega_1(\mu)-c_1](1-\mu^2)\, \mathrm{d}\mu \nonumber\\ &=\frac{A_+}{(1-\mu_\star^2)^{{3}/{2}}}\int^{1}_{\mu_\star}2\varOmega_0[\varOmega_1(\mu)-c_1](1-\mu^2)\, \mathrm{d}\mu, \end{align}

this providing two relations between $\alpha _1$, $A_+$ and $A_-$. Next, we match $H$: at $\mu -\mu _\star =O(\varepsilon ^{{1}/{2}})$, the outer solution of $H$ is dominated by the tilting mode $H_0$, so that the matching condition is

(4.23)\begin{equation} \mathcal{H} \big|_{\eta\rightarrow \pm \infty}=H_0 \, \big|_{\mu\rightarrow \mu_\star^\pm}. \end{equation}

According to (4.10) and (4.20),

(4.24a,b)\begin{equation} \alpha_2=A_+\sqrt{1-\mu_\star^2} , \quad \frac{-\mathrm{i}{\rm \pi}\alpha_1}{(\varOmega_{1\star}-c_1)|\beta_\star'|}+\alpha_2=A_-\sqrt{1-\mu_\star^2} . \end{equation}

Combining (4.22) and (4.24), a non-trivial solution for $\alpha _1$, $\alpha _2$, $A_-$ and $A_+$ yields an equation that determines the eigenvalue $c_1$ for $\mathrm {Im}\, c_1>0$,

(4.25)$$\begin{gather} \frac{1}{\displaystyle\int_{{-}1}^{\mu_\star} 2\varOmega_0 [c_1-\varOmega_1(\mu)](1-\mu^2)\,\mathrm{d}\mu}+ \frac{1}{\displaystyle\int_{\mu_\star}^1 2\varOmega_0 [c_1-\varOmega_1(\mu)](1-\mu^2)\,\mathrm{d}\mu} \nonumber\\ =\frac{\mathrm{i}{\rm \pi}}{|\beta_\star'|(1-\mu_\star^2)^2(c_1-\varOmega_{1\star})} . \end{gather}$$

An example of the comparison between the asymptotic and numerical solutions is shown in figure 5. Note that if (4.25) yields a solution with $\mathrm {Im}\, c_1<0$, it is not a valid normal mode solution, since it contradicts our branch cut selection (4.19) based on $\mathrm {Im}\, c_1>0$.

It is useful to express (4.25) using the original variables $c$ and $\varOmega$ instead of $c_1$ and $\varOmega _1$ via (4.8a,b). The resulting equation is presented as (5.1) at the beginning of the next section, where we also discuss its implications. This is an equation for the eigenvalue $c$ with $c_{i}>0$, and with $\varOmega _\star$ equal to $\varOmega (\mu_\star)$. We have replaced $\varOmega _0=\varOmega _\star$ for the leading-order solid body rotation; as the shear is weak, $\varOmega$ is approximately constant everywhere. Equation (5.1) is derived under the condition that $\varOmega$ and $\beta$ are of the same order and the shear of $\varOmega$ is small compared with both of these. In Appendix A we also present the analysis for the situation where $\varOmega$ as a whole is small compared with $\beta$. The derivation is a little different but the final result remains the same as (5.1), and so this equation is generally applicable as long as the shear $\varOmega '$ is weak compared with the magnetic field $\beta$.

5. Results and discussion

The result of our analysis is the following implicit equation for the complex wave speed $c$, taken to have a positive imaginary part $c_{i}>0$ that gives the growth rate of a mode:

(5.1)\begin{equation} \frac{c-\varOmega_\star}{\displaystyle\int_{{-}1}^{\mu_\star} [c-\varOmega(\mu)] (1-\mu^2)\,\mathrm{d}\mu}+ \frac{c-\varOmega_\star}{\displaystyle\int_{\mu_\star}^1 [c-\varOmega(\mu)](1-\mu^2)\,\mathrm{d}\mu}=\frac{2\mathrm{i} {\rm \pi}\varOmega_\star}{|\beta_\star'|(1-\mu_\star^2)^2} . \end{equation}

Here we recall that the magnetic field profile $\beta (\mu )$ has a single, simple zero at the critical latitude given by $\mu = \mu _\star$, where the gradient $\beta '_\star \neq 0$ and the angular velocity is $\varOmega_\star$. The equation is valid provided the shear $\varOmega '(\mu )$ of the angular velocity profile is small compared with the magnetic field.

5.1. General results

Equation (5.1) has a relatively simple form and we can use it to gain significant insights into the instability properties for general profiles of $\varOmega$ and $\beta$. First, we observe that $\beta$ only enters this equation through $\beta '_\star =\beta '(\mu _\star )$, where $\mu _\star$ is the location where $\beta =0$. The other properties of $\beta$ (e.g. the value of $\beta$ at other latitudes) do not affect (5.1). This is a curious property, because the MHD instability is global, yet the magnetic field only affects the instability through its local behaviour in the critical layer. To test this finding, we consider three different profiles: $\beta =\mu$, $\sin \mu$ and $e^\mu -1$. They all have $\mu _\star =0$ and $\beta '_\star =1$, and so the same asymptotic result for $c$ given by (5.1). For the zonal flow, we select $\varOmega =1-s\mu ^2$ with $s=0.12$ and $s=0.06$; the case of $s=0.12$ has been used by Cally (2001) as a model for solar differential rotation. The results of the eigenvalues are displayed in table 1. For comparison, the asymptotic solution (5.1) is given in the last row. We see that the eigenvalues $c$ for the various profiles are indeed close, and interestingly, to a much higher degree than the precision of the asymptotic solution. We summarise this conclusion as follows: the magnetic field profile only affects the instability through the location where it passes though zero and the value of its gradient there.

Table 1. Numerical solutions for the eigenvalue $c=c_{r}+\mathrm {i}c_{i}$ for $m=1$, $\varOmega =1-s\mu ^2$, $s=0.12$ and 0.06, and three profiles of $\beta$ with $\beta =0$ and $\beta '=1$ at $\mu _\star =0$. The solutions are computed by a shooting method. The three profiles have the same asymptotic prediction for $c$, given by (5.1) and shown in the last row of the table.

To proceed further, we rewrite (5.1) as

(5.2) \begin{equation} \frac{c-\varOmega_\star}{(c-\varOmega_\star)I_-+J_-}+\frac{c-\varOmega_\star}{(c-\varOmega_\star)I_++J_+}=\mathrm{i}Q, \end{equation}

where

(5.3a–e) \begin{align} I_-&=\int_{{-}1}^{\mu_\star} (1-\mu^2)\,\mathrm{d}\mu=\mu_\star{-}\frac{1}{3}\mu_\star^3+\frac{2}{3},\quad I_+{=}\int_{\mu_\star}^1(1-\mu^2)\,\mathrm{d}\mu={-}\mu_\star+\frac{1}{3}\mu_\star^3+\frac{2}{3},\nonumber\\ J_-&=\int_{{-}1}^{\mu_\star}(\varOmega_\star{-}\varOmega)(1-\mu^2)\,\mathrm{d}\mu,\quad J_+{=}\int_{\mu_\star}^1(\varOmega_\star{-}\varOmega)(1-\mu^2)\,\mathrm{d}\mu,\nonumber\\ Q &=\frac{2{\rm \pi} \varOmega_\star}{|\beta_\star'|(1-\mu_\star^2)^2} . \end{align}

In general, (5.2) is a quadratic equation for $c$ and its solution is

(5.4)\begin{align} c&=\varOmega_\star+\frac{1} {2(I-\mathrm{i}QI_+I_-)}\left\{\vphantom{{}^{{1}/{2}}}-J+\mathrm{i}Q(I_-J_++I_+J_-)\right. \nonumber\\ &\left. \quad \pm \left[J^2-Q^2(I_+J_- -I_-J_+)^2+2\mathrm{i}Q(I_+J_- -I_-J_+)(J_+ - J_-)\right]^{{1}/{2}}\right\}, \end{align}

where

(5.5a,b)\begin{equation} I=I_++I_-{=}\frac{4}{3}, \quad J=J_++J_-{=}\int_{{-}1}^1(\varOmega_\star - \varOmega)(1-\mu^2)\, \mathrm{d}\mu. \end{equation}

Of the two solutions given by (5.4), only those with $c_{i}>0$ are valid. At this point, it is not straightforward to obtain an exact condition for $c_{i}>0$ to hold, but it is easy to find a sufficient condition as follows. Given that the second line of (5.4) has both positive and negative signs, if its first line already has a positive imaginary part, then at least one of the solutions has positive $c_{i}$. Therefore, a sufficient condition for instability is

(5.6)\begin{equation} \mathrm{Im} \left(\frac{\mathrm{i}Q(I_-J_++I_+J_-)-J}{2(I-\mathrm{i}QI_+I_-)} \right)=\frac{{\rm \pi}\varOmega_\star(I_-^2J_++I_+^2J_-)}{|\beta_\star'|(1-\mu_\star^2)^2\left(I^2+Q^2I_+^2I_-^2\right)}>0. \end{equation}

A simple example is the situation where $|\varOmega _\star |$ is the maximum of $|\varOmega |$, then according to (5.3c,d), both $\varOmega _\star J_+$ and $\varOmega _\star J_-$ are positive, and (5.6) is guaranteed. This leads to the conclusion: if the angular velocity of the zonal flow is greatest at the critical level, then the flow is unstable. This agrees with the statement earlier that the semicircle rules always allow instability for such flows, discussed at the end of § 3. For model solar differential rotation profiles, $\varOmega$ is indeed largest at the equator, and so provided $\beta$ passes (transversely) through zero there, the flow is always unstable. The instability induced by the magnetic profile $\beta =\sigma \mu$ that we showed in figures 2 and 3 belongs to this category.

5.2. The solution for $\mu _\star =0$ and even $\varOmega$

The solutions (5.4) may be further simplified and yield transparent results when the critical level is located at the equator ($\mu _\star =0$) and $\varOmega$ has even or odd symmetry, which we discuss in §§ 5.2 and 5.3, respectively.

The case in which $\mu _\star =0$ and $\varOmega$ is an even function of $\mu$ is perhaps most relevant to the Sun, and therefore, most studies on the clamshell instability focus on this case (for example, Gilman & Fox 1997; Cally 2001; Cally et al. 2003; Miesch 2007; Miesch et al. 2007). We then have $J_+=J_-$ and $I_+=I_-=2/3$ and the two solutions of (5.4) are

(5.7a,b)\begin{equation} c=\varOmega_\star+\frac{3\mathrm{i}Q J_+}{6-2\mathrm{i}Q}\quad \mathrm{and} \quad c=\varOmega_\star{-}\frac{3}{2} J_+. \end{equation}

The first solution (5.7a) is complex, and may give an unstable mode. Using the original variables, (5.7a) becomes

(5.8)\begin{equation} c=\varOmega_\star{-} \frac{3{\rm \pi} \varOmega_\star\displaystyle\int_0^1 (\varOmega_\star{-}\varOmega)(1-\mu^2)\,\mathrm{d}\mu}{2{\rm \pi}\varOmega_\star+3\mathrm{i}|\beta_\star'|}. \end{equation}

Its imaginary part is

(5.9)\begin{equation} c_{i}=\frac{9{\rm \pi}\varOmega_\star|\beta_\star'|\displaystyle\int_0^1(\varOmega_\star - \varOmega)(1-\mu^2)\,\mathrm{d}\mu}{4{\rm \pi}^2\varOmega_\star^2+9{\beta_\star'}^2} . \end{equation}

For the standard profile

(5.10a,b)\begin{equation} \varOmega=r-s\mu^2,\quad \beta=\sigma\mu, \end{equation}

we find that

(5.11a,b)\begin{equation} c_{r}=r-\frac{4{\rm \pi}^2r^2s}{20{\rm \pi}^2r^2+45\sigma^2} ,\quad c_{i}=\frac{6{\rm \pi}|\sigma|rs}{20{\rm \pi}^2 r^2+45 \sigma^2}. \end{equation}

As discussed in the last paragraph of § 4.2.3, this solution is accurate in the limit when $s$ is small or $|\sigma |$ is large. Note that large $|\sigma |$ at fixed $r$ and $s$ corresponds to the situation where $\varOmega$ as a whole is weak compared with $\beta$. The detailed analysis for this case is shown in Appendix A. The results of (5.11) are plotted in figure 6 by dashed lines, and may be compared with the numerical solutions shown by solid lines. We see that the asymptotic solution gives very good predictions in general, and that these become more precise as $s$ decreases or as $\sigma$ increases.

Figure 6. Solutions of $c=c_{r}+\mathrm {i}c_{i}$ vs $\sigma$ for the profiles $\varOmega =r-s\mu ^2$ and $\beta =\sigma \mu$ with $r=1$, $s=0.05,0.1,0.2$ and wavenumber $m=1$. The asymptotic solution (5.11) is plotted by dashed lines, and the numerical solutions are plotted by solid lines.

With the asymptotic solution for the growth rate $c_{i}$ given by (5.11b), we may address the question of whether instability persists when the parameters approach limiting values. Gilman & Fox (1997) raised the question of whether there is a lower limit of positive $s$ and an upper limit of $\sigma$ for the instability to take place. These thresholds do not seem to exist according to their numerical solutions, but the unstable mode becomes more and more singular at the critical level as $s$ decreases or $\sigma$ increases, causing numerical difficulties. Our asymptotic solution can easily address this problem: (5.11b) indicates that such limits indeed do not exist: instead, as $s\rightarrow 0^+$ or $\sigma \rightarrow \infty$, $c_{i}$ remains positive at $O(s)$ or $O(\sigma ^{-1})$ provided that $rs>0$. The clamshell instability is therefore quite different from hydrodynamic shear instability on a sphere, which requires the shear to exceed a threshold ($s/r>0.29$, Watson 1981) to overcome the stabilising effect of the rotation. The fact that the instability survives for an arbitrarily strong magnetic field is also surprising, but we note that regardless of the strength of $\sigma$, the magnetic field always vanishes at $\mu _\star =0$, and it is this feature that plays a fundamental role in inducing the instability.

The solution (5.9) can also provide insights into the instability for general flow profiles, not only those related to solar differential rotation. From the condition of $c_{i}>0$, we have: for $\mu _\star =0$ and an even profile of $\varOmega (\mu )$, the condition for instability is

(5.12)\begin{equation} \varOmega_\star\int_{0}^1(\varOmega_\star - \varOmega)(1-\mu^2)\,\mathrm{d}\mu>0. \end{equation}

This indicates that the flow is prone to instability when the angular velocity $\varOmega _\star$ at the critical level (the equator in this case) is large compared with $\varOmega (\mu )$ on the rest of the sphere. If $|\varOmega _\star |$ is the largest among all $|\varOmega |$, then the flow is definitely unstable. Interestingly, (5.12) only involves the hydrodynamic shear, and $\beta$ does not affect this condition once $\mu _\star =0$ is set. Equation (5.12) is also quite different from conditions for hydrodynamic shear instability: the latter usually involve constraints on the curvature of the basic-flow profile (cf. Rayleigh's inflection-point theorem), but (5.12) does not involve $\varOmega ''$ at all.

We can also derive a bound for the growth rate from (5.9), namely

(5.13)\begin{equation} c_{i}< \frac{9{\rm \pi} |\varOmega_\star|\,|\beta_\star'|\max|\varOmega_\star - \varOmega|\displaystyle\int_0^1(1-\mu^2)\,\mathrm{d}\mu}{12{\rm \pi}|\varOmega_\star\beta_\star'|}\leq \frac{1}{2}\max |\varOmega_\star{-}\varOmega|, \end{equation}

where the inequality $x^2 + y^2 \geq 2|xy|$ has been used in the denominator. Again, once $\mu _\star =0$ is set by the magnetic field, this bound only involves the hydrodynamic shear. We note that the semicircle rules studied in § 3 suggest that the magnetic field may increase the radii of the semicircles and, hence, the bound for the growth rate, but this does not happen in (5.13). However, we also note that the semicircle rules apply to general velocity and field profiles, and (5.13) is the result for the more specific situation in which the rotation profile is even with weak shear, and the magnetic profile passes through zero at one location.

Clearly when $\varOmega$ is even, $\varOmega _1$ is also even, and from (4.22) we have $A_-=-A_+$. Thus, the critical layer makes the tilting modes opposite on the two hemispheres (as also shown in figure 3c), which explains the typical clamshell pattern shown in figure 2.

Finally, we comment on the other solution (5.7b). The physical meaning of this solution is that it makes the singularity of $H_1'$ given in (4.14) vanish. Thus, to leading order, the weak shear does not trigger the singularity of the critical level of the tilting mode. One may need to go to higher orders in the asymptotic expansion, which may contain potential singularities and yield an even smaller $c_{i}$. Our numerical solution suggests that, for the basic state profiles (5.10a,b), the higher-order solution of (5.7b) has $c_{i}<0$ and is therefore not a valid normal mode solution. However, when the critical level $\mu _\star$ is slightly off the equator, (5.7b) becomes an unstable normal mode, as we will show in § 5.4.

5.3. The solution for $\mu _\star =0$ and odd $\varOmega -\varOmega _\star$

The situation in which $\mu _\star =0$ and the shear profile $\varOmega -\varOmega _\star$ is an odd function of $\mu$ is less relevant to the Sun, but as a basic model it is still of interest to fluid mechanics and we may draw useful general conclusions in this case. Here we have $J_+=-J_-$, $I_-=I_+=2/3$ and the imaginary part of (5.4) can be simplified to

(5.14)\begin{equation} c_{i}={\pm}\sqrt{\frac{81|Q|J_+^2}{8(Q^2+9)\left(\sqrt{Q^2+9}+|Q|\right)}}. \end{equation}

Except for the special case of $J_+=0$ (for which we would need to pursue higher orders of the asymptotic expansion), there is always a positive $c_{i}$ solution and so, surprisingly, we may conclude that if $\mu _\star =0$ and $\varOmega (\mu )-\varOmega _\star$ is odd, the flow is always unstable. Bounding the denominator of (5.14) from below via $(Q^2+9)(\sqrt {Q^2+9}+|Q|) > 9 \times 2 |Q|$, we obtain a bound for the unstable growth rate,

(5.15)\begin{equation} c_{i}< \tfrac{3}{4}|J_+| \le\tfrac{1}{2}\max |\varOmega_\star{-}\varOmega|. \end{equation}

Interestingly, this bound is the same as (5.13), but we expect it to be looser since we have bounded positive terms by zero. To give a concrete example for this instability, we consider

(5.16a,b)\begin{equation} \varOmega=r+s\mu,\quad \beta=\sigma\mu, \end{equation}

where $\varOmega$ features a ‘linear shear’ profile, analogous to Couette flow. Then we have

(5.17a,b)\begin{equation} Q=\frac{2{\rm \pi} r}{|\sigma|} , \quad J_+{=}-\frac{s}{4} . \end{equation}

The results of (5.14) with (5.17a,b) are shown figure 7(a), where they are compared with the numerical solution. Again good agreement is found and the agreement improves as $s$ becomes smaller or $\sigma$ becomes larger. The behaviour of $c_{i}$ is similar to the previous case of even $\varOmega$, as is $c_{r}$ (not shown). An example for the solution of $H$ is shown in figure 7(b). Because $\varOmega$ as a whole is neither even nor odd, $H$ has no symmetry property either.

Figure 7. Instability of the profiles $\varOmega =r+s\mu$ and $\beta =\sigma \mu$ with $r=1$, $m=1$. (a) Curves of $c_{i}$ vs $\sigma$ for $s=0.05,0.1,0.2$. The asymptotic solution (5.14) with (5.17a,b) is plotted by dashed lines, and the numerical solution is plotted by solid lines. (b) Numerical solution showing $H$ for $r=\sigma =1$ and $s=0.1$, with eigenvalue $c=0.96+0.0087\mathrm {i}$.

5.4. An example for $\mu _\star \neq 0$

When the critical level is off the equator (i.e. $\mu _\star \neq 0$), there is no obvious symmetry property that can simplify the asymptotic solution for $c$ given by (5.4). We have not been able to obtain general conclusions regarding the condition for instability in this case, but the asymptotic solution can still be helpful in understanding numerical results. As an example, we consider

(5.18a,b)\begin{equation} \varOmega=1-0.1\mu^2,\quad \beta=\mu-d. \end{equation}

The zonal flow features the solar differential rotation as before, and $\beta$ is a linear profile with the critical level located at $\mu _\star =d$. When $d=0$, we recover the standard configuration (5.10a,b).

The asymptotic solution for $c$ computed from (5.4) versus $d$ is shown in figure 8 (dashed lines), where it is compared with the numerical solution (solid lines). In figure 8(a) we also plot the value of $\varOmega$ at the critical level, $\varOmega _\star =1-0.1 d^2$ (dotted line). The main feature of this asymmetric case is that both solutions of (5.4) can have positive $c_{i}$, and thus, there are two branches of unstable modes. When $d$ slightly departs from zero, the mode that corresponds to (5.7b) becomes unstable (red dashed line), and its growth rate dominates over the other unstable mode for a large range of $d$. On the other hand, as $d$ increases, the unstable mode that corresponds to (5.7a) (blue dashed line) is weakened significantly, and disappears at $d\approx 0.37$. Again good agreement is found between the asymptotic and numerical results, but interestingly, there is a topological difference between them: the asymptotic solution predicts that when the two eigenvalues are close, the curves of $c_{i}$ intersect while those of $c_{r}$ avoid the intersection, while the opposite is true of the numerical solution.

Figure 8. Eigenvalue $c=c_{r}+\mathrm {i}c_{i}$ vs $d$ when $\varOmega =1-0.1\mu ^2$, $\beta =\mu -d$ and $m=1$. The solid lines represent the numerical solution, and the dashed lines are the results of the asymptotic solution (5.4). The dotted line in panel (a) represents the rotation rate $\varOmega$ at the critical level $\mu _\star =d$, i.e. $\varOmega _\star =1-0.1d^2$. The inset in panel (a) shows the region where the two eigenvalues become close.

It appears from figure 8(a) that $c_{r}$ is always smaller than $\varOmega _\star$, which demands an explanation. When the mode of the blue dashed line (or the red solid line) has $c_{r}$ approach $\varOmega _\star$ at $d\approx 0.37$ (left panel), the corresponding $c_{i}$ (right panel) approaches zero. Thus, $\varOmega _\star$ appears to be an upper bound for $c_{r}$ for unstable modes. There is an underlying reason for this phenomenon, related to the conservation of angular momentum, as we will explain in § 5.5.

The eigenfunctions of the two unstable modes at $d=0.36$ are plotted in figure 9. Figure 9(a) is the mode with the smaller growth rate $c_{i}$. In fact, we have chosen $d$ such that this mode is almost as close to a neutral mode as we can compute numerically. The very small $c_{i}$ makes the critical layer have a very fine structure. There is a significant difference in the amplitudes of the tilting modes on the two sides of the critical layer. Figure 9(a) has a much larger, stronger tilting mode to the left of the critical level, while figure 9(b) has the opposite feature.

Figure 9. Two numerical solutions of the eigenfunction $H$ corresponding to figure 8 at $d=0.36$. The eigenvalues are (a) $c=0.985+2.05\times 10^{-4}\mathrm {i}$, (b) $c=0.965+8.36\times 10 ^{-3}\mathrm {i}$. The case of panel (a) is the one with almost the smallest growth rate that we can compute.

5.5. The conservation of angular momentum

The asymptotic analysis clearly indicates that the critical layer plays a fundamental role in making the flow unstable. In our previous studies of instability induced by critical layers (Riedinger & Gilbert 2014; Wang & Balmforth 2018; Wang, Gilbert & Mason 2022), conservation of momentum provides a useful tool for understanding the mechanism of the instability. Indeed, it has been found that the critical layer provides a source of mean-flow momentum, which drives the exponential growth of the outer flow. In the current problem in spherical geometry, the relevant conservation law is that of angular momentum (2.36), namely,

(5.19)\begin{equation} \int_0^{\rm \pi} 2{\rm \pi} \sin^2\theta\, \frac{\partial \Delta U}{\partial t}\,\mathrm{d}\theta=0. \end{equation}

It is of interest to understand how this conservation is achieved, i.e. how different regions contribute to the integral and balance each other.

Substituting (2.20a–d), (2.24) and (2.28) into (2.34), we can derive the rate of change of angular momentum per latitude,

(5.20)\begin{equation} 2{\rm \pi}\sin^2\theta \frac{\partial \Delta U}{\partial t}= \frac{\partial}{\partial \theta}\frac{\partial L}{\partial t}, \end{equation}

with

(5.21)\begin{equation} \frac{\partial L}{\partial t}=\left[4{\rm \pi}(1-\mu^2)\left\{ \left[|\varOmega-c|^2-\beta^2\right] \mathrm{Im}(H{H'}^*)-c_{i}\varOmega' |H|^2\right\}\right] {\rm e}^{2c_{i} t}. \end{equation}

Here $L(\theta,t)$ represents the total mean-flow angular momentum between the north pole and co-latitude $\theta$. Recall that the primes denote derivatives with respect to $\mu =\cos \theta$. For the clamshell instability studied above, $\varOmega -c$, $\varOmega '$ and $c_{i}$ are all small at order $O(\varepsilon )$, so that outside, or at the edge of the critical layer, we have

(5.22)\begin{equation} \frac{\partial L}{\partial t}={-}4{\rm \pi}(1-{\mu}^2)\beta^2 \mathrm{Im}(H{H'}^*) {\rm e}^{2c_{i}t}, \end{equation}

to leading order of $\varepsilon$. This result corresponds to the fact that the Maxwell stress, i.e. $\overline {a_\ell b_\ell }$ in (2.34), has the dominant contribution to the mean-flow response, whilst the Reynolds stress $\overline {u_\ell v_\ell }$ has a minor effect due to the weak shear. We can then study the integral of (5.20) over $\theta$ in different regions. We define $\theta _\star$ as the value of $\theta$ at the critical level (i.e. $\cos \theta _\star =\mu _\star$) and set $\Delta =O(\varepsilon ^{{1}/{2}})$ as the half-thickness of the critical layer. Then, outside the critical layer the integrals are

(5.23a)\begin{gather} \left.\int_0^{\theta_\star{-}\Delta}2{\rm \pi}\sin^2\theta \frac{\partial \Delta U}{\partial t}\, \mathrm{d}\theta=\frac{\partial L}{\partial t}\right|_{\mu=\mu_\star+\Delta{\sin\theta_\star}}=8{\rm \pi} c_{i}\varOmega_0|A_+|^2{\rm e}^{2c_{i}t}I_+, \end{gather}

(5.23b)\begin{gather} \left.\int_{\theta_\star+\Delta}^{\rm \pi} 2{\rm \pi}\sin^2\theta \frac{\partial \Delta U}{\partial t}\, \mathrm{d}\theta={-}\frac{\partial L}{\partial t} \right|_{\mu=\mu_\star{-}\Delta {\sin\theta_\star}}=8{\rm \pi} c_{i}\varOmega_0|A_-|^2{\rm e}^{2c_{i}t}I_-, \end{gather}

where $I_-$ and $I_+$ are the positive quantities defined in (5.3) and we have used $H=H_0$ and $H'=\varepsilon H_1'$ from § 4.2.1 as the leading-order approximation on the edge of the critical layer. Inside the critical layer, using the inner solution given in § 4.2.2, we find the integral to be

(5.24)\begin{align} \left.\int_{\theta_\star{-}\Delta}^{\theta_\star+\Delta}2{\rm \pi}\sin^2\theta \frac{\partial \Delta U}{\partial t}\, \mathrm{d}\theta={-}\frac{\partial L}{\partial t} \right|_{\eta\rightarrow -\infty}^{\infty}=\frac{4\varepsilon{\rm \pi}^2(1-\mu_\star^2)|\alpha_1|^2}{|\beta_\star'|}\mathrm{Re} \left(\frac{1}{c_1-\varOmega_{1\star}}\right){\rm e}^{2c_{i}t}. \end{align}

Applying the relations between the constants $\alpha _1$, $\alpha _2$, $A_-$ and $A_+$ given in (4.22) and (4.24), we can show that the value of (5.24) exactly cancels the sum of (5.23a) and (5.23b), and results in (5.19) being satisfied. The critical layer thus provides a source of angular momentum that balances that of the outer flow. We note that without the critical-layer angular momentum (5.24), the only possibility that (5.23a) and (5.23b) could add up to zero is if $c_{i}=0$, i.e. the tilting modes by themselves have to be neutral modes. This demonstrates how the angular momentum provided by the critical layer is necessary to drive the instability.

We may gain some further insights by considering the sign of the mean angular momentum inside and outside the critical layer. Without loss of generality, we consider the case in which $\varOmega _0>0$. Then both (5.23a) and (5.23b), giving the mean angular momentum of the tilting components, are positive, so that the contribution from the critical layer (5.24) must be negative to make the conservation law (5.19) possible. Since

(5.25)\begin{equation} \mathrm{Re} \left(\frac{1}{c_1-\varOmega_{1\star}}\right)=\frac{c_{1\mathrm{r}}-\varOmega_{1\star}}{(c_{1\mathrm{r}}-\varOmega_{1\star})^2+c_{1i}^2}, \end{equation}

we require

(5.26)\begin{equation} c_{1{r}}<\varOmega_{1\star}\quad \mathrm{or} \quad c_{r}<\varOmega_\star. \end{equation}

This means that, for any unstable mode, the real part of the phase velocity must be smaller than the velocity of the zonal flow at the critical level. We can verify that all of the solutions we have showed previously satisfy this condition. For example, in (5.11a) we have $c_{r}< r$ for $s>0$, and in figure 8(a) the curve of $\varOmega _\star$ is always above that of $c_{r}$. In the situation of figure 8(a), we may also view (5.26) as a necessary condition for the existence of an unstable mode: when $c_{r}$ is about to exceed $\varOmega _\star$ at $d\approx 0.37$, the unstable mode disappears.

We may undertake a similar analysis for the conservation of the mean toroidal field as shown by (2.37), but we were not able to obtain straightforward general conclusions. This is mainly because the local integral of $\partial _t \Delta A$ in the critical layer has a less transparent expression. Nevertheless, we document these results in Appendix B for the readers’ interest.

6. Conclusions

We have studied the linear instability of two-dimensional MHD flows on a sphere, a problem with potential application to the instability of the solar tachocline. We derived semicircle rules for the complex phase velocity, which provide rigorous bounds for general flow and field profiles. The terms arising purely from the spherical geometry bring new features to the problem. We used two bounding methods, which provide two versions of the semicircle rules, each of which may be tighter for certain types of flows. We also found that the magnetic field may increase the radii of the semicircles, which does not happen in the case of Cartesian geometry (Hughes & Tobias 2001).

We then undertook an analytical study of the ‘clamshell instability’. Previous studies have found that the instability tilts the basic magnetic field lines on the two hemispheres in opposite directions, giving a pattern of an opening clamshell (Cally 2001; Cally et al. 2003). We studied this instability theoretically through an asymptotic analysis in the limit of weak shear of the basic zonal flow. We found that if the basic zonal flow is a pure solid body rotation, there exists an eigenmode that slightly tilts the entire magnetic field and makes it rotate with the zonal flow. We refer to this disturbance as a ‘tilting mode’. Including an additional weak shear in the zonal flow excites the critical level of the tilting mode, located at the node of the sheared field profile. Disturbances exhibit a strong singular behaviour near the critical level, inside the critical layer. We found that the critical layer reverses the direction of tilting and makes the flow unstable.

Through matching the tilting mode and the critical layer, we derived the asymptotic solution for the complex phase velocity, from which we obtained properties of the instability for general profiles. Our investigations indicate that the magnetic field only affects the instability through the location of the critical level and its gradient at the critical level; the other details of the field profile do not matter. A sufficient condition for instability is that the critical level is located where the angular velocity of the zonal flow is greatest. When the zonal flow is even about the equator and the critical level is on the equator, we derived a simple expression for the unstable growth rate, which indicates that the flow is susceptible to instability when the angular velocity at the critical level is large compared with that on the rest of the sphere. When the shear of the zonal flow is odd and the critical level is on the equator, the flow is always unstable. A simple bound for the unstable growth rate was derived for these two types of flows with even or odd symmetry properties. In the absence of symmetry, when the critical level is off the equator, there can be two branches of unstable modes. The results of the asymptotic solution are in good agreement with the numerical solutions.

A mechanism for the instability has been provided via the conservation of angular momentum. The critical layer provides a source of angular momentum, which must be balanced by a corresponding sink for the surrounding tilting mode. In order that the angular momentum of the tilting mode and critical layer have opposite signs, the phase velocity of the unstable mode must be smaller than the velocity of the zonal flow at the critical level.

Our study reveals several problems that are left for future research. We found that the magnetic field can increase the radii of the semicircles (over those for the purely hydrodynamic flow), but we have not yet found an unstable mode that resides in this new region. It is interesting to investigate whether it can be found for different flow and field profiles. In addition, when the magnetic field is relatively strong, the destabilising effect of the field is always associated with an increase of the semicircle radius. It remains an open question as to whether there is a deeper link between these observations. The theories of Thuburn & Haynes (1996), Sasaki et al. (2012) and Deguchi (2021) that may provide different semicircles could be possible routes to approach this problem.

The clamshell instability we studied occurs for idealised MHD flows with weak shear and a strong field that vanishes at one location. There are flows with field-induced instabilities that do not belong to this category; for example, strong shear combined with a weak magnetic field (Gilman & Fox 1997; Cally 2001), magnetic field profiles with multiple zero points (Dikpati & Gilman 1999) and narrow bands of magnetic field (Dikpati & Gilman 1999; Cally et al. 2003). Dissipation may also be of potential interest: our results indicate that, for ideal MHD, the field lines on the two sides of the critical layers are separated, so it would be interesting to explore the details of the reconnection caused by diffusion as seen in the simulation of Cally (2001). It has also been found that diffusion may destabilise the flow even when the zonal flow has no shear (Sharif & Jones 2005). Beyond the incompressible MHD setting, instabilities also arise in shallow-water MHD systems in spherical geometry (Márquez-Artavia et al. 2017; Gilman & Dikpati 2002). A deeper understanding might be gained by studying whether a similar asymptotic analysis is applicable to these instabilities. We also plan to explore the theory of nonlinear critical layers, to understand the saturation of growing disturbances.