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 (Reference Gilman and Fox1997). 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 (Reference Watson1981), 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 (Reference Gilman and Dikpati2002). 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 Reference Gilman and Fox1999), and the other is a shooting method (Dikpati & Gilman Reference Dikpati and Gilman1999). The method of Legendre polynomial expansion computes all of the eigenvalues but it can be expensive. Similarly to Gilman & Fox (Reference Gilman and Fox1999), 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 (Reference Howard1961). To tackle the spherical geometry, we will mainly follow the theory of Watson (Reference Watson1981) 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 (Reference Gilman and Fox1997). 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 (Reference Howard1961) and Watson (Reference Watson1981), 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 (Reference Hughes and Tobias2001), 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 (Reference Watson1981) 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 (Reference Gilman1967), Cally (Reference Cally2000) and Hughes & Tobias (Reference Hughes and Tobias2001) for MHD instability in Cartesian geometry, which can be further reduced to the theory of Howard & Gupta (Reference Howard and Gupta1962) when the magnetic field is uniform. When both limits are applied, (3.17) and (3.20) both reduce to the semicircle rule of Howard (Reference Howard1961). 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 (Reference Hughes and Tobias2001), 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 Reference Márquez-Artavia, Jones and Tobias2017), 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 (Reference Hughes and Tobias2001), 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 (Reference Gilman and Fox1997). For the basic state magnetic field, we consider two examples: the first is a ‘linear shear’ profile considered by Gilman & Fox (Reference Gilman and Fox1997) and Cally (Reference Cally2001),

(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 (Reference Dikpati and Gilman1999) and Cally (Reference Cally2001),

(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 Reference Watson1981); 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 (Reference Gilman and Fox1997), 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 (Reference Howard1961), 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 (Reference Watson1981), but Thuburn & Haynes (Reference Thuburn and Haynes1996) and Sasaki et al. (Reference Sasaki, Takehiro and Yamada2012) 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 (Reference Deguchi2021) 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 (Reference Cally2001) and Cally et al. (Reference Cally, Dikpati and Gilman2003) 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 (Reference Watson1981) (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 (Reference Bernoff and Lingevitch1994).

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 (Reference Cally2001) (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 (Reference Cally2001) 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$.