# Two-dimensional ideal magnetohydrodynamic waves on a rotating sphere under a
non-Malkus field: I. Continuous spectrum and its ray-theoretical
interpretation

Ryosuke Nakashima${{\dagger}}$∗ and Shigeo Yoshida${{\ddagger}}$
${{\dagger}}$Faculty of Science, Kyushu University, Fukuoka, Japan
${{\ddagger}}$Department of Earth and Planetary Sciences, Faculty of Science,
Kyushu University, Fukuoka, Japan
∗Corresponding author. Email<EMAIL_ADDRESS>
###### Abstract

Two-dimensional (2D) ideal incompressible magnetohydrodynamic (MHD) linear
waves at the surface of a rotating sphere are studied as a model imitating the
outermost Earth’s core or the solar tachocline. This thin conducting layer is
permeated by a toroidal magnetic field whose magnitude depends only on the
latitude. The Malkus background field, which is proportional to the sine of
the colatitude, gives two well-known groups of branches on which Alfvén waves
gradually become fast or slow magnetic Rossby (MR) waves as the field
amplitude decreases. For non-Malkus fields, we show that the associated
eigenvalue problems can yield a continuous spectrum instead of Alfvén and slow
MR discrete modes. Critical latitudes attributed to the Alfvén resonance wipe
out these discrete eigenvalues and produce an infinite number of singular
eigenmodes. The theory of slowly varying wave trains in an inhomogeneous
magnetic field shows that a wave packet related to this continuous spectrum
propagates toward a critical latitude corresponding to the wave and is
eventually absorbed there. The expected behaviour that the retrograde
propagating packets which pertain to the continuous spectrum approach the
latitudes from the equatorial side and that the prograde ones approach there
from the polar side is consistent with the profiles of their eigenfunctions
shown by our numerical calculations. Further in-depth discussions of the
Alfvén continuum would progress the theory of “wave-mean field interaction” in
the MHD system and one’s understanding of the dynamics in such thin layers.

This is a non-peer reviewed preprint that has been submitted to Geophysical &
Astrophysical Fluid Dynamics (https://www.tandfonline.com/journals/ggaf20). It
is deposited under the terms of the Creative Commons Attribution-
NonCommercial-ShareAlike License (https://creativecommons.org/licenses/by-nc-
sa/4.0/). This license requires that reusers give credit to the creator. It
allows reusers to distribute, remix, adapt, and build upon the material in any
medium or format, for noncommercial purposes only. If others modify or adapt
the material, they must license the modified material under identical terms.

Keywords: Magnetic Rossby waves; Continuous modes; Stably stratified layer;
Geomagnetic variations;

##  1 Introduction

Geomagnetic and geodetic variations may partly be accounted for by MHD waves
within the Earth’s outer core (e.g. Gillet et al., 2021; Triana et al., 2021;
Hori et al., 2023). For example, Hide (1966) suggested that the westward drift
of the geomagnetic field may originate from slow MR waves in the liquid core.
Braginsky (1970) ascribed both the $60$-year length-of-day and geomagnetic
variations to the torsional oscillations. If this type of attribution is
substantiated, the comparison between observations and the theory leads to
some inferences on the relevant physical quantities. Zatman and Bloxham (1997)
accepted Braginsky (1970)’s explanation to infer that the magnitude of the
cylindrical radial field is about $0.2\,\mathrm{mT}$. This suggestion was
questioned by Gillet et al. (2010), who associated the torsional oscillations
(or torsional waves) with six-year length-of-day signals to conclude that its
magnitude is about $2\,\mathrm{mT}$.

A stably stratified layer at the top of the outer core has been proposed on
various grounds. Simple thermal (Gubbins et al., 1982) and compositional
(Braginsky, 1984) stratification were considered first, and the proposed
physical mechanisms of stratification have been more sophisticated lately. For
the thermal stratification, a subadiabatic temperature gradient due to the
core’s high thermal conductivity (Pozzo et al., 2012; Zhang et al., 2022) has
been invoked, and for the compositional stratification, barodiffusion and
chemical interactions between the core and the mantle (Buffett and Seagle,
2010; Gubbins and Davies, 2013; Brodholt and Badro, 2017; Davies et al.,
2018), and a remnant of the Moon-creating impact (Landeau et al., 2016) were
proposed. Seismological evidence for such a layer has been controversial (e.g.
Helffrich and Kaneshima, 2010; Kaneshima, 2018; Irving et al., 2018; van Tent
et al., 2020). Regional, instead of global, stratification arising from the
core-mantle boundary (CMB) heterogeneity was also proposed (Mound et al.,
2019). The obscure properties of the layer may be able to be inferred from the
identification of sources of geomagnetic and geodetic signatures from the
core.

Magnetic-Archimedes-Coriolis (MAC) and MR waves in such a stably stratified
layer at the top of the core have often been invoked as possible causes of
geomagnetic fluctuations (Braginsky, 1993, 1998, 1999; Buffett, 2014; Chulliat
et al., 2015; Buffett et al., 2016; Knezek and Buffett, 2018; Chi-Durán et
al., 2021). If the stratified layer exists atop the core, vigorous convection
prevailing in the bulk of the core overshoots the interface between the bulk
and the layer (Takehiro and Lister, 2001; Gastine et al., 2020) and can excite
MHD waves that travel in the layer (Jaupart and Buffett, 2017; Buffett and
Knezek, 2018; Couston et al., 2017; Bouffard et al., 2022).

Two-dimensional ideal incompressible MHD linear waves within a thin conducting
fluid layer over a rotating sphere are addressed by the present paper. A non-
uniform toroidal ($\phi$-directional, $\phi$ being the longitude) magnetic
field is imposed on the fluid film as a background field, and the layer
rotates almost rigidly with the sphere. To examine Hide (1966)’s suggestion, a
similar problem was first considered by Stewartson (1967), in which the
azimuthal main field is spatially uniform. Such a toroidal field does not
vanish even at the north and south poles ($\theta=0$ and $\uppi$, where
$\theta$ is the colatitude). One should therefore express the basic field as
$B_{0\phi}(\theta)=B_{0}\mathcal{B}\sin\theta$ with a constant $B_{0}$ and a
bounded continuous function $\mathcal{B}(\cos\theta)$. Although a smooth
background field at the poles also requires the function $\mathcal{B}$ to
satisfy that
$\mathcal{B}\sin\theta=-\sum_{n=1}^{\infty}t_{n}(\ud\mathrm{P}_{n}/\ud\theta)$,
where $t_{n}$ ($n=1,2,\ldots$) are constants and $\mathrm{P}_{n}$ is the
Legendre polynomial of degree $n$, we now disregard this condition. The
simplest profile $B_{0\phi}=B_{0}\sin\theta$ in this expression, or
$\mathcal{B}=1$, is equivalent to the well-known Malkus field (e.g. Malkus,
1967; Finlay, 2008). This elementary 2D model is useful to qualitatively
understand the wave dynamics in the thin layer at the top of the core.

The dispersion relation of 2D waves under the Malkus background field should
be reviewed before moving on to our main topic. Under the Malkus field, the
local Alfvén wave velocities are
$V_{\mathrm{A}\phi}(\theta)=(B_{0}/\sqrt{\rho_{0}\mu_{\mathrm{m}}})\sin\theta$,
where $\rho_{0}$ is a uniform density and $\mu_{\mathrm{m}}$ is a constant
magnetic permeability. If other effects are absent, local Alfvén waves with
speeds proportional to $\sin\theta$ can collectively travel in concert in the
azimuthal direction like the rigid body rotation. This reduces our
perturbation equations to a regular Sturm-Liouville problem because interior
poles vanish (Boyd, 1981), as will be shown in Section 2. If the Malkus field
permeates a rigidly rotating thin layer on top of a sphere with a radius
$R_{0}$, the dispersion relation for the 2D ideal incompressible MHD wave
(Zaqarashvili et al., 2007; Márquez-Artavia et al., 2017) is given by

$\lambda\,=\,\frac{-m\pm m\sqrt{1+4\alpha^{2}n(n+1)[n(n+1)-2]}}{2n(n+1)}\,,$
(1)

where $\lambda\equiv\omega/2\varOmega_{0}$ is the angular frequency
nondimensionalised by the double of the constant rotation rate
$\varOmega_{0}$. In the above expression, $m$ denotes the (positive) zonal
wavenumber, $\alpha\equiv
B_{0}/2\varOmega_{0}R_{0}\sqrt{\rho_{0}\mu_{\mathrm{m}}}$ is the (signed)
Lehnert number, and $n$ is the degree of the associated Legendre polynomial
$\mathrm{P}_{n}^{m}(\cos\theta)$. Although the Lehnert number is often defined
as $Le\equiv 2|\alpha|$ in geophysics literature, we follow the precedents
referred to above. The double sign in (1) corresponds to prograde and
retrograde propagating Alfvén waves $\lambda\simeq\pm
m|\alpha|\sqrt{1-2/n(n+1)}$ for large $|\alpha|$, and the two classes of MR
waves as $|\alpha|\to 0$: the slow mode $\lambda\simeq m\alpha^{2}[n(n+1)-2]$
peculiar to the MHD system and the fast mode $\lambda\simeq-m/n(n+1)$. The
last is closely related to (hydrodynamic) Rossby waves celebrated in
meteorology and oceanography (e.g. Longuet-Higgins, 1968). Figure 1 shows that
the exponent of $|\alpha|$ in a branch of $\lambda$, or
$(\upartial\log|\lambda|/\upartial\log|\alpha|)$, can be used to distinguish
among the three types of waves derived from (1).

Figure 1: Nondimensional angular frequency $\lambda$ calculated from (1) as a
function of the absolute value $|\alpha|$ of the Lehnert number when the zonal
wavenumber $m=1$. The left panel shows retrograde propagating waves
($\lambda<0$) and the right corresponds to prograde ones ($\lambda>0$). Fast
MR waves for $n=1,2,\ldots,10$ are represented as the warm colour curves,
while slow MR waves for $n=2,\ldots,10$ are depicted as the cold colour
curves. As $|\alpha|\to\infty$, these modes approach the lines
$|\lambda|\propto|\alpha|^{1}$ of Alfvén waves except for $m=n=1$.

The simplest equatorially-antisymmetric non-Malkus field
$B_{0\phi}=B_{0}\sin\theta\cos\theta$, or $\mathcal{B}=\cos\theta$, should be
more appropriate for the Earth’s core than the pure Malkus field. In the
current study, we devote our attention to the case of this basic field profile
because MHD waves under a non-Malkus field are poorly understood. Note that
the field is eastward in the northern hemisphere when $B_{0}>0$, and westward
when $B_{0}<0$. It is important that $|\alpha|$ is much smaller than unity in
the Earth’s core. If the field strength $|B_{0}|$ inside the core is several
millitesla (Gillet et al., 2010), one can estimate that $|\alpha|\approx
10^{-4}$ with $\varOmega_{0}\approx+0.729\times 10^{-4}\,\mathrm{s^{-1}}$,
$R_{0}\approx 3480\,\mathrm{km}$, $\rho_{0}\approx
10^{4}\,\mathrm{kg\,m^{-3}}$ and $\mu_{\mathrm{m}}\approx 1.26\times
10^{-6}\,\mathrm{H\,m^{-1}}$.

Non-Malkus fields complicate our linear wave problem because regular singular
points appear in their equations. At the singular colatitudes
$\theta=\theta_{\mathrm{c}}$ resulting from the Alfvén resonance, the
azimuthal phase velocity $\omega/(m/R_{0}\sin\theta_{\mathrm{c}})$ of a wave
is equal to the local Alfvén velocity
$V_{\mathrm{A}\phi}(\theta_{\mathrm{c}})$ at the latitude (e.g. Uberoi, 1972;
Goedbloed and Poedts, 2004). For example, such latitudes ordinarily exist for
low-frequency neutral modes under a non-Malkus configuration whose function
$\mathcal{B}$ has at least one zero in $0<\theta<\uppi$. Even if such points
exist for a given angular frequency $\lambda$, one can, however, construct its
eigenmode when its eigenfunction is permitted to have singular profiles only
at the points, or rather such singular modes are required for completeness
(Van Kampen, 1955; Case, 1959; Barston, 1964). This means that its spectrum
should include a continuous range,
$\min(\mathcal{B}^{2})\leq\lambda^{2}/m^{2}\alpha^{2}\leq\max(\mathcal{B}^{2})$.
Even though it makes no sense to extract an eigenmode (hereinafter referred to
as a “continuous eigenmode” or “continuous mode”) from this continuous
spectrum as its eigenfunction is singular, the integration (or the
superposition) with respect to $\lambda$ in which the integrand is the
eigenfunctions weighted by a coefficient depending only on $\lambda$, can
constitute a physically relevant solution. Recall that it is not always true
that the Fourier transform of a well-behaved function is not pathological.
Solutions composed only of the continuous eigenmodes do not show the behaviour
being typical of collective oscillations by normal discrete eigenvalues but
represent transient growth of initial disturbances (Farrell, 1982) and non-
diffusive attenuation (e.g. Adam, 1986). The latter involves the phase mixing
and the Landau damping, which cause algebraic decay (Case, 1960; Balmforth and
Morrison, 1995b) and exponential decay (Landau, 1946; Briggs et al., 1970;
Sedláček, 1971; Tataronis and Grossmann, 1973), respectively. These
unintuitive aspects due to the advent of continuous spectra have long been
discussed in various research areas: inviscid shear flows, collisionless
plasmas, ideal MHD systems, 2D vortices and differentially rotating disks,
self-gravitating systems, and the Kuramoto model (e.g. Adam, 1986; Balmforth
and Morrison, 1995a; Strogatz, 2000; Barré et al., 2015).

The reasons that we chose the azimuthal field
$B_{0\phi}=B_{0}\mathcal{B}\sin\theta$ are as follows, although the dominant
field may be radial in the Earth’s outermost core due to the small
conductivity in the mantle (e.g. Knezek and Buffett, 2018). First of all, this
model is simpler and can be straightforwardly extended into the “MHD shallow
water” system (Zaqarashvili et al., 2007, 2009; Heng and Spitkovsky, 2009;
Zaqarashvili et al., 2011; Márquez-Artavia et al., 2017), which was first
introduced by Gilman (2000). In contrast with the toroidal field, radial ones
normally bring vertical derivatives into their governing equations, and the
operator is obviously incompatible with a 2D model, which is computationally
undemanding. Secondly, critical latitudes often reside in equations of similar
eigenvalue problems even when a radial background field depending on $\theta$
passes through a thin layer (e.g. see the equations in Buffett and Matsui,
2019). Although approximations that are apparently appropriate for each target
of study sometimes obliterate the latitudes (Zaqarashvili, 2018; Buffett and
Matsui, 2019), it is not clear whether the simplifications that drop the
Alfvén resonance are always valid or not. Thus, we need to scrutinise the
influences of the critical latitudes on eigenmodes. Lastly, Hardy et al.
(2020) reported that slight vertical motions under strong stratification can
enhance toroidal fields because of the “Malkus constraint.” This suggestion
supports our choice as a reasonable main field.

Related instability problems may pertain to the dynamics of the solar
tachocline underlying the convection zone. Because unstable modes can easily
be picked up even in the presence of continuous spectra, this circumstance is
significantly different from that of the study of linear waves. Gilman and Fox
(1997, 1999a); Gilman and Fox (1999b) and subsequent work (Dikpati and Gilman,
1999; Gilman and Dikpati, 2000; Zaqarashvili et al., 2010a) examined the
“joint instability” or “magnetic Rossby wave instability” arising in the
tachocline, which can occur in the MHD system accompanied by latitudinal
differential rotation. According to them, some non-axisymmetric infinitesimal
perturbations are likely to destabilise the coexistence of the solar-like
angular velocity profile deduced from helioseismic observations and a variety
of plausible toroidal field configurations such as broad profiles written as
$B_{0\phi}=(B_{0}+B_{1}\cos^{2}\theta)\sin\theta\cos\theta$ (in which $B_{1}$
is also constant) or latitudinally localised field bands expressed by Gaussian
functions. These global unstable modes may play an important role in the
persistence of such a thin shear layer via latitudinal transport of angular
momentum (Spiegel and Zahn, 1992) and put an upper limit on the strength of a
toroidal field stored within the layer through the $\omega$-effect (Arlt et
al., 2007a, b). The nonlinear evolution of the modes has been developed by
Cally (2001) and Cally et al. (2003); Cally et al. (2004), who found the novel
“clamshell” and “tipping” patterns. Additionally, viscosity and magnetic
diffusion were also introduced in the radial (Dikpati et al., 2004) and
horizontal (Sharif and Jones, 2005) directions in an anisotropic manner.

Even beyond the strict 2D model, stability analyses of the tachocline have
been conducted. The MHD shallow water equations (Gilman, 2000) and the 3D thin
shell model or “MHD hydrostatic primitive” equations (Miesch and Gilman, 2004)
were newly proposed for evaluating the impacts of subadiabatic stratification
and weak vertical displacement of fluid particles. Gilman and Dikpati (2002)
and Dikpati et al. (2003) demonstrated that the combinations of the
differential rotation and the toroidal fields in the shallow model become
unstable easily again and that the growing perturbations have non-zero kinetic
helicities, which are related to the $\alpha$-effect (Moffatt, 1978).
Furthermore, unstable MR waves in the layer may have been causing some
periodicities detected in the solar activity (Zaqarashvili et al. 2010a, b,
2015; Dikpati et al. 2017; Gachechiladze et al. 2019; see Zaqarashvili et al.
2021, for a review). For instance, nonlinear development in the MHD shallow
water system (Dikpati et al., 2017, 2018a, 2018b) indicated that MR waves
exchange angular momentum with mean fields and that their wave patterns
deformed by consequent reconstructions of the mean profiles can trigger
nonlinear quasi-periodic oscillations. In addition, searches for not only non-
axisymmetric unstable modes (Cally, 2003; Gilman et al., 2007; Kitchatinov and
Rüdiger, 2008) but also axisymmetric ones (Cally et al., 2008; Dikpati et al.,
2009) were performed in extensive studies of the 3D thin shell model. The
studies currently include nonlinear simulations of these growth modes (Miesch
et al., 2007; Hollerbach and Cally, 2009) and linear stability analyses taking
vertical profiles of the differential rotation and the background toroidal
field into account (Arlt et al., 2007a, b).

Critical lines (or levels, latitudes, and so on) and their concomitant
continuous modes have important bearings even on unstable modes which are
outside of the objective of this paper. Despite the fact that the
eigenfrequencies of unstable modes are complex numbers, their eigenfunctions
are affected by the positions of critical points, as demonstrated by, e.g.
Gilman and Fox (1999a) and Wang et al. (2022a). Moreover, a neutral mode
sometimes interacts with continuous modes when the branch of the neutral one
overlaps the continuous spectrum. This leads to the appearance of a pair of
unstable and decay modes (Iga, 1999; Taniguchi and Ishiwatari, 2006) or non-
modal growth (Heifetz et al., 2020). On the other hand, continuous spectra
sometimes cover and hide the branches of such interacting neutral modes (Iga,
2013). Therefore, we also want to detail the critical latitudes and the
neutral continuous eigenmodes deeply in advance to understand linear stability
in similar problems (e.g Wang et al., 2022b), though our current problem does
not have such unstable modes.

Non-ideality and non-linearity are also intimately related to critical points.
Small viscosity and magnetic diffusion transform a singular point on the real
axis on the “fictional” complex plane of the involved spatial coordinate into
a complex turning point, because the diffusion terms include higher-order
derivatives (e.g. Drazin and Reid, 1981; Shivamoggi, 1992). Since the
eigenfunctions that we are going to seek are functions on the real axis on the
complex plane, they then become non-singular. This seems at a glance to mean
that we have to recover these ignored damping terms. However, very weak
diffusions only give rise to thin boundary layers around the turning points
(and by any walls), and the eigenfunctions would still be similar to the
profiles of continuous eigenmodes obtained from the non-diffusive limit
sufficiently outside the layers. Although decaying normal modes stemming from
measurable dissipations (Steinolfson, 1985; Gizon et al., 2020) and nonlinear
boundary layers (e.g. Tung, 1979; Maslowe, 1986) may also be somewhat
important to our problem, they are beyond the scope of this article.

This paper is organised as follows. In Section 2, we shall derive the
governing equations for our problem, and then present a method to seek
eigenmodes numerically by the associated Legendre polynomial expansion in
Section 2.1. Section 3 gives numerical solutions for the case when the
background field is expressed as the simplest equatorially-antisymmetric non-
Malkus field $\mathcal{B}=\cos\theta$. In Sections 3.1 and 3.2, the structures
of the obtained eigenfunctions are examined outside and near critical
latitudes, respectively. In Section 4, we also conduct the numerical
integration of ray-tracing equations at large wavenumbers (e.g. Bardsley and
Davidson, 2017; Teruya et al., 2022) in order to interpret our eigenvalue
problem from a different angle. This approach tracks paths of wave packets
migrating with their group velocities. Note that this is different from
Dikpati et al. (2020)’s calculations of the Stokes drift, which are
trajectories of fluid particles advected by oscillatory flow induced by MR
waves. Finally, we conclude in Section 5.

##  2 Mathematical formulation

We now begin with a description of the governing equations for the 2D ideal
incompressible MHD on a rotating sphere. In the spherical coordinate system
$(r,\theta,\phi)$, the 2D vorticity and 2D uncurled induction equations on the
spherical surface $r=R_{0}$ (e.g. Raphaldini and Raupp, 2020) are

$\displaystyle\frac{\mathrm{D}\zeta}{\mathrm{D}t}\,-\,\frac{2\varOmega_{0}\sin\theta}{R_{0}}u_{\theta}\,$
$\displaystyle=\,\frac{(\bm{B}\bm{\cdot}\bm{\nabla}_{\mathrm{H}})(\mu_{\mathrm{m}}J)}{\rho_{0}\mu_{\mathrm{m}}}\,,$
(2a) $\displaystyle\frac{\mathrm{D}A}{\mathrm{D}t}\,$ $\displaystyle=\,0\,,$
(2b)

where
$\bm{\nabla}_{\mathrm{H}}\equiv(\hat{\bm{e}}_{\theta}/R_{0})(\upartial/\upartial\theta)+(\hat{\bm{e}}_{\phi}/R_{0}\sin\theta)(\upartial/\upartial\phi)$
is the horizontal nabla operator, and the material derivative is expressed as
$(\mathrm{D}/\mathrm{D}t)\equiv(\upartial/\upartial
t)+\bm{u}\bm{\cdot}\bm{\nabla}_{\mathrm{H}}$. The velocity field
$\bm{u}=(0,u_{\theta},u_{\phi})$ (relative to the rotating frame of reference)
and the magnetic field $\bm{B}=(0,B_{\theta},B_{\phi})$ are assumed to have
only the horizontal components within the thin layer, for the sake of
simplicity. The radial components of the vorticity and electrical current in
(2) are defined as

$\displaystyle\zeta\,$
$\displaystyle\equiv\,\frac{1}{R_{0}\sin\theta}\left[\frac{\upartial(u_{\phi}\sin\theta)}{\upartial\theta}-\frac{\upartial
u_{\theta}}{\upartial\phi}\right]\,=\,-\nabla_{\mathrm{H}}^{2}\psi\,,$ (3a)
$\displaystyle J\,$
$\displaystyle\equiv\,\frac{1}{\mu_{\mathrm{m}}}\frac{1}{R_{0}\sin\theta}\left[\frac{\upartial(B_{\phi}\sin\theta)}{\upartial\theta}-\frac{\upartial
B_{\theta}}{\upartial\phi}\right]\,=\,-\mu_{\mathrm{m}}^{-1}\nabla_{\mathrm{H}}^{2}A\,,$
(3b)

in which we introduce the stream function $\psi$ and the magnetic vector
potential $\bm{A}=(A,0,0)$. Owing to the solenoidal conditions of the fields
$\bm{u}$ and $\bm{B}$, the scalars $\psi$ and $A$ are related to the two
vectors as

$\displaystyle
u_{\theta}\,=\,\frac{1}{R_{0}\sin\theta}\frac{\upartial\psi}{\upartial\phi}\,,$
$\displaystyle\qquad
u_{\phi}\,=\,-\frac{1}{R_{0}}\frac{\upartial\psi}{\upartial\theta}\,,$ (4a)
$\displaystyle B_{\theta}\,=\,\frac{1}{R_{0}\sin\theta}\frac{\upartial
A}{\upartial\phi}\,,$ $\displaystyle\qquad
B_{\phi}\,=\,-\frac{1}{R_{0}}\frac{\upartial A}{\upartial\theta}\,.$ (4b)

Using these expressions, the governing equations (2) become

$\displaystyle\frac{\upartial(\nabla_{\mathrm{H}}^{2}\psi)}{\upartial
t}\,+\,\mathcal{J}(\nabla_{\mathrm{H}}^{2}\psi,\psi)\,+\,\frac{2\varOmega_{0}}{R_{0}^{2}}\frac{\upartial\psi}{\upartial\phi}\,$
$\displaystyle=\,\frac{1}{\rho_{0}\mu_{\mathrm{m}}}\mathcal{J}(\nabla_{\mathrm{H}}^{2}A,A)\,,$
(5a) $\displaystyle\frac{\upartial A}{\upartial t}\,+\,\mathcal{J}(A,\psi)\,$
$\displaystyle=\,0\,,$ (5b) where the operator $\mathcal{J}(f,g)$ for any two
scalar functions $f$ and $g$ is defined as
$\mathcal{J}(f,g)\,\equiv\,\frac{1}{R_{0}^{2}\sin\theta}\left(\frac{\upartial
f}{\upartial\theta}\frac{\upartial g}{\upartial\phi}-\frac{\upartial
f}{\upartial\phi}\frac{\upartial g}{\upartial\theta}\right)\,.$ (5c)

In what follows, waves of small amplitude are considered in this system. We
introduce a small positive parameter $\varepsilon$ ($\ll 1$) which represents
the amplitude of waves to rewrite $\psi$ and $A$ as

$\psi\,=\,\varepsilon\psi_{1}\,+\,\mathrm{O}(\varepsilon^{2})\,,\qquad
A\,=\,A_{0}(\theta)\,+\,\varepsilon a_{1}\,+\,\mathrm{O}(\varepsilon^{2})\,.$
(6)

The basic state is assumed to be the rigid body rotation $\psi_{0}\equiv 0$
with a latitude-dependent toroidal field $B_{0\phi}(\theta)$. With the
expression of $B_{0\phi}$ given in the previous section, one obtains $(\ud
A_{0}/\ud\theta)=-R_{0}B_{0}\mathcal{B}\sin\theta$. Then, the equations (5)
become

$\displaystyle\left(\frac{\upartial}{\upartial
t}\nabla_{\mathrm{H}}^{2}+\frac{2\varOmega_{0}}{R_{0}^{2}}\frac{\upartial}{\upartial\phi}\right)\psi_{1}\,$
$\displaystyle=\,\frac{1}{R_{0}\rho_{0}\mu_{\mathrm{m}}\sin\theta}\left[B_{0\phi}\nabla_{\mathrm{H}}^{2}-\frac{1}{R_{0}}\frac{\ud(\mu_{\mathrm{m}}J_{0})}{\ud\theta}\right]\frac{\upartial
a_{1}}{\upartial\phi}\,+\,\mathrm{O}(\varepsilon)\,,$ (7a)
$\displaystyle\frac{\upartial a_{1}}{\upartial t}\,$
$\displaystyle=\,\frac{B_{0\phi}}{R_{0}\sin\theta}\frac{\upartial\psi_{1}}{\upartial\phi}\,+\,\mathrm{O}(\varepsilon)\,,$
(7b)

where
$J_{0}(\theta)=\mu_{\mathrm{m}}^{-1}(1/R_{0}\sin\theta)[\ud(B_{0\phi}\sin\theta)/\ud\theta]$
is the background electrical current.

The normal mode approach is valid for our purposes. For a given azimuthal
wavenumber $m$ and an angular frequency $\omega$ determined later by a
dispersion relation, which becomes (1) for the Malkus field $\mathcal{B}=1$
and is going to be numerically sought in Section 3 when
$\mathcal{B}=\cos\theta$, we postulate that
$\psi_{1}(\theta,\phi,t)\equiv\mathrm{Re}[\tilde{\psi}(\mu;m,\omega)\mathrm{e}^{\mathrm{i}\varphi(\phi,t)}]=\tilde{\psi}\mathrm{e}^{\mathrm{i}\varphi}/2+\mathrm{c.c.}$
and $a_{1}\equiv\mathrm{Re}[\tilde{a}\mathrm{e}^{\mathrm{i}\varphi}]$, where
$\mu=\cos\theta$ and $\varphi\equiv m\phi-\omega t=m\phi-\lambda\tau$ is the
phase of waves with the nondimensional time $\tau\equiv 2\varOmega_{0}t$. We
use $\mathrm{c.c.}$ for the complex conjugate of the preceding terms. Upon
substituting this ansatz into (7), we immediately get

$\displaystyle(-\lambda\nabla_{\mathrm{h}}^{2}+m)\tilde{\psi}\,$
$\displaystyle=\,m|\alpha|\left\\{\mathcal{B}\nabla_{\mathrm{h}}^{2}-\frac{{}^{2}[\mathcal{B}(1-\mu^{2})]}{\ud\mu^{2}}\right\\}\left[\frac{\operatorname{sgn}(\alpha)\tilde{a}}{\sqrt{\rho_{0}\mu_{\mathrm{m}}}}\right]\,,$
(8a)
$\displaystyle-\lambda\left[\frac{\operatorname{sgn}(\alpha)\tilde{a}}{\sqrt{\rho_{0}\mu_{\mathrm{m}}}}\right]\,$
$\displaystyle=\,m|\alpha|\mathcal{B}\tilde{\psi}\,,$ (8b)

in which $\nabla_{\mathrm{h}}^{2}\equiv R_{0}^{2}\nabla_{\mathrm{H}}^{2}$ is
the dimensionless horizontal Laplacian. From now on, our interest is limited
to the case where $m\neq 0$. The above equations are easily transformed into a
single ordinary differential equation in the form

$\frac{\ud}{\ud\mu}\left[\varLambda(1-\mu^{2})\frac{\ud\tilde{\psi}}{\ud\mu}\right]\,-\,\left\\{\frac{m^{2}\varLambda}{1-\mu^{2}}+m\left[\lambda+2m\alpha^{2}\mathcal{B}\frac{\ud(\mathcal{B}\mu)}{\ud\mu}\right]\right\\}\tilde{\psi}\,=\,0\,,$
(9)

where the factor
$\varLambda(\mu)\equiv\lambda^{2}-m^{2}\alpha^{2}\mathcal{B}^{2}$ is crucial
to our problem. If there exist real values $\mu$ which satisfy
$\varLambda(\mu)=0$ within the interval $-1<\mu<1$, which are hereinafter
denoted by $\mu_{\mathrm{c}}$, those points are interior poles depending on
$\lambda$. The poles yield continuous spectra, as we stated in the
introduction section. On the other hand, the Malkus field $\mathcal{B}=1$
obviously produces no singular points except for the endpoints $\mu=\pm 1$.
Dividing (9) by the factor $\varLambda$ ($\neq 0$), one can obtain its
dispersion relation (1) without any hurdles because (9) is then reduced to the
associated Legendre differential equation. Note that if $\lambda$ and
$\tilde{\psi}$ are an eigenvalue and its corresponding eigenfunction of (9),
respectively, the same holds true for their complex conjugates.

The existence of continuous spectra can be justified when the function
$\varLambda$ has zeros in $-1<\mu<1$. Since (9) is a second-order differential
equation, it should have two linearly independent solutions. Let
$\tilde{\psi}_{\mathrm{I}}$ be the non-singular one of the solutions, and let
the other be expressed in the form
$\tilde{\varPsi}_{\mathrm{I\\!I}}=\tilde{\psi}_{\mathrm{I}}\int^{\mu}f(\mu_{*})\ud\mu_{*}$
with a function $f(\mu)$. Substituting the form of
$\tilde{\varPsi}_{\mathrm{I\\!I}}$ into (9), we find

$\tilde{\psi}\,=\,C_{\mathrm{I}}\tilde{\psi}_{\mathrm{I}}\,+\,C_{\mathrm{I\\!I}}\tilde{\varPsi}_{\mathrm{I\\!I}}\,,\qquad\tilde{\varPsi}_{\mathrm{I\\!I}}(\mu)\,=\,\tilde{\psi}_{\mathrm{I}}(\mu)\int^{\mu}\frac{\ud\mu_{*}}{(1-\mu_{*}^{2})\varLambda(\mu_{*})\tilde{\psi}_{\mathrm{I}}^{2}(\mu_{*})}\,,$
(10)

with both $C_{\mathrm{I}}$ and $C_{\mathrm{I\\!I}}$ being constants. Let us
examine what the integral becomes if the interval of integration in (10) pass
over zeros of $\varLambda$ (we assume that $\tilde{\psi}_{\mathrm{I}}^{2}\neq
0$ within the interval). With an arbitrary starting point $\mu_{0}$ of the
integration interval (in the following equation, we suppose that
$\mu_{0}<\operatorname{min}(\mu,\mu_{\mathrm{c}})$ and that at most one
singular point exists between $\mu_{0}$ and $\mu$), we can express the second
solution of the linearly independent set as the improper integral

$\tilde{\varPsi}_{\mathrm{I\\!I}}(\mu)\,=\,\begin{cases}\displaystyle\tilde{\psi}_{\mathrm{I}}(\mu)\int_{\mu_{0}}^{\mu}\frac{\ud\mu_{*}}{(1-\mu_{*}^{2})\varLambda(\mu_{*})\tilde{\psi}_{\mathrm{I}}^{2}(\mu_{*})}&(\mu<\mu_{\mathrm{c}})\\\\[14.22636pt]
\displaystyle\tilde{\psi}_{\mathrm{I}}(\mu)\left[\lim_{\varDelta_{1}\to+0}\int_{\mu_{0}}^{\mu_{\mathrm{c}}-\varDelta_{1}}\frac{\ud\mu_{*}}{(1-\mu_{*}^{2})\varLambda(\mu_{*})\tilde{\psi}_{\mathrm{I}}^{2}(\mu_{*})}\right.\\\
\displaystyle\hskip
85.35826pt\,+\,\left.\lim_{\varDelta_{2}\to+0}\int_{\mu_{\mathrm{c}}+\varDelta_{2}}^{\mu}\frac{\ud\mu_{*}}{(1-\mu_{*}^{2})\varLambda(\mu_{*})\tilde{\psi}_{\mathrm{I}}^{2}(\mu_{*})}\right]&(\mu_{\mathrm{c}}<\mu)\\\
\end{cases}\,.$ (11)

We may then deduce from (11) that (10) is also written as

$\tilde{\psi}\,=\,C_{\mathrm{I}}\tilde{\psi}_{\mathrm{I}}\,+\,C_{\mathrm{I\\!I}}\tilde{\psi}_{\mathrm{I\\!I}}\,+\,\tilde{\psi}_{\mathrm{I}}\sum_{i}C_{\mathrm{I\\!I\\!I},i}\mathrm{H}(\mu-\mu_{\mathrm{c},i})\,,$
(12a) where the integer $i$ is the index of the singular latitudes and
$\mathrm{H}(\mu)$ is the step function. In addition, we have introduced the
new second solution
$\tilde{\psi}_{\mathrm{I\\!I}}(\mu)\,\equiv\,\tilde{\psi}_{\mathrm{I}}(\mu)\,\mathcal{P}\\!\\!\int^{\mu}\frac{\ud\mu_{*}}{(1-\mu_{*}^{2})\varLambda(\mu_{*})\tilde{\psi}_{\mathrm{I}}^{2}(\mu_{*})}\,,$
(12b) in which $\mathcal{P}$ denotes the Cauchy principal value. In the
vicinity of $\mu=\mu_{\mathrm{c}}$, we know that
$\mathcal{B}^{2}=\mathcal{B}_{\mathrm{c}}^{2}+(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}(\mu-\mu_{\mathrm{c}})+(\mathcal{B}_{\mathrm{c}}^{2})^{\prime\prime}(\mu-\mu_{\mathrm{c}})^{2}/2+(\mathcal{B}_{\mathrm{c}}^{2})^{\prime\prime\prime}(\mu-\mu_{\mathrm{c}})^{3}/6+\mathrm{O}(|\mu-\mu_{\mathrm{c}}|^{4})$,
where
$\mathcal{B}_{\mathrm{c}}^{2}\equiv\mathcal{B}^{2}(\mu_{\mathrm{c}})=\lambda^{2}/m^{2}\alpha^{2}$,
$(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}\equiv(\ud\mathcal{B}^{2}/\ud\mu)|_{\mu=\mu_{\mathrm{c}}}$
and the like. Now, in this paper, we will narrow down a target only to the
case when $(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}$ does not vanish. Using
this series expansion, one obtains the magnitude $C_{\mathrm{I\\!I\\!I},i}$ of
the discontinuities at the latitudes in the form
$C_{\mathrm{I\\!I\\!I},i}\,\equiv\,-\frac{C_{\mathrm{I\\!I}}I_{i}}{(1-\mu_{\mathrm{c},i}^{2})m^{2}\alpha^{2}(\mathcal{B}_{\mathrm{c},i}^{2})^{\prime}\tilde{\psi}_{\mathrm{I}}^{2}(\mu_{\mathrm{c},i})}$
(12c) with the integral $I_{i}$ given by
$I_{i}\,\equiv\,\lim_{\varDelta,\varDelta_{1,i}\to+0}\int_{\mu_{\mathrm{c},i}-\varDelta}^{\mu_{\mathrm{c},i}-\varDelta_{1,i}}\frac{\ud\mu_{*}}{\mu_{*}-\mu_{\mathrm{c},i}}\,+\,\lim_{\varDelta,\varDelta_{2,i}\to+0}\int_{\mu_{\mathrm{c},i}+\varDelta_{2,i}}^{\mu_{\mathrm{c},i}+\varDelta}\frac{\ud\mu_{*}}{\mu_{*}-\mu_{\mathrm{c},i}}\,=\,\lim_{\varDelta_{1,i},\varDelta_{2,i}\to+0}\ln\frac{\varDelta_{1,i}}{\varDelta_{2,i}}\,.$
(12d)

This integral can be an arbitrary number, depending on how we take the limit
$\varDelta_{1,i},\varDelta_{2,i}\to+0$. Indeed, Van Kampen (1955) suggested
that, in the problem of plasma oscillations, his counterpart of $I$ may be
considered as an arbitrary parameter, which can be determined by a
normalization condition of his counterpart of $\tilde{\psi}$, or the
distribution function of plasma. As seen in (12a), the existence of these more
than two linearly independent solutions despite (9) being a second-order
differential equation can lead to an excess of arbitrary coefficients which
should be adjusted to satisfy boundary conditions. This excessive freedom
results in continuous spectra.

Another explanation for the appearance of continuous spectra is derived from
the condition that should be satisfied by the solution of (9) in the vicinity
of the singular point. On integrating (9) with respect to $\mu$ over the
narrow range sandwiched between $\mu_{\mathrm{c}}-\varDelta_{1}$ and
$\mu_{\mathrm{c}}+\varDelta_{2}$ with $\varDelta_{1}$, $\varDelta_{2}\to+0$,
one get

$\lim_{\varDelta_{1}\to+0}\varDelta_{1}\left.\frac{\ud\tilde{\psi}}{\ud\mu}\right|_{\mu=\mu_{\mathrm{c}}-\varDelta_{1}}\,+\,\lim_{\varDelta_{2}\to+0}\varDelta_{2}\left.\frac{\ud\tilde{\psi}}{\ud\mu}\right|_{\mu=\mu_{\mathrm{c}}+\varDelta_{2}}\,\to\,0\,,$
(13a) provided that the condition
$\lim_{\varDelta_{1},\varDelta_{2}\to+0}\int_{\mu_{\mathrm{c}}-\varDelta_{1}}^{\mu_{\mathrm{c}}+\varDelta_{2}}|\tilde{\psi}|\ud\mu\,\to\,0$
(13b)

is fulfilled. The fact that $\tilde{\psi}_{\mathrm{I}}$ is surely a solution
for (9) and that $\tilde{\psi}_{\mathrm{I}}\mathrm{H}(\mu-\mu_{\mathrm{c}})$
always fulfills (13a) shows that the third term in (12a) is a weak solution
for (9) and $C_{\mathrm{I\\!I\\!I}}$ should then be an undetermined parameter.

An alternative survey of their structures at the critical latitudes is
conducted with the Frobenius method (e.g. Braun, 1975). Suppose that a power
series of the form
$\tilde{\psi}^{(\mathrm{c})}\equiv\sum_{k=0}^{\infty}a_{k}(\mu-\mu_{\mathrm{c}})^{k+\varrho}$
($a_{0}\neq 0$) is a solution for (9) around a critical latitude
$\mu_{\mathrm{c}}$, in which $\varrho$ is a root of the indicial equation for
(9). On substituting this assumption for $\tilde{\psi}$ in (9), one obtains
the equation $\varrho^{2}=0$ from its leading order term. As a result, the two
linearly independent solutions near the latitude are given on the basis of the
way to deal with the repeated root $\varrho=0$ by

$\displaystyle\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}\,$
$\displaystyle\equiv\,1\,+\,\sum_{k=1}^{\infty}a_{k}(\mu-\mu_{\mathrm{c}})^{k}\,,$
(14a) $\displaystyle\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}\,$
$\displaystyle\equiv\,\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}\ln|\mu-\mu_{\mathrm{c}}|\,+\,\sum_{k=1}^{\infty}b_{k}(\mu-\mu_{\mathrm{c}})^{k}\,.$
(14b)

Some of the expansion coefficients are found with slightly tedious but
standard manipulations as $a_{1}=D_{1}$,
$a_{2}=(D_{1}^{2}+2D_{1}D_{2}+D_{3})/4$, $b_{1}=-2D_{1}+D_{2}$,
$b_{2}=(-3D_{1}^{2}-2D_{1}D_{2}+2D_{2}^{2}-D_{3}+2D_{4})/4$, where

$\displaystyle
D_{1}\,\equiv\,-\frac{\lambda/m\alpha^{2}+2\mathcal{B}_{\mathrm{c}}^{2}+(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}\mu_{\mathrm{c}}}{(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}(1-\mu_{\mathrm{c}}^{2})}\,,$
$\displaystyle\quad
D_{2}\,\equiv\,\frac{2(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}\mu_{\mathrm{c}}-(\mathcal{B}_{\mathrm{c}}^{2})^{\prime\prime}(1-\mu_{\mathrm{c}}^{2})/2}{(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}(1-\mu_{\mathrm{c}}^{2})}\,,$
$\displaystyle
D_{3}\,\equiv\,\frac{m^{2}(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}/(1-\mu_{\mathrm{c}}^{2})-3(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}-(\mathcal{B}_{\mathrm{c}}^{2})^{\prime\prime}\mu_{\mathrm{c}}}{(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}(1-\mu_{\mathrm{c}}^{2})}\,,$
$\displaystyle\quad
D_{4}\,\equiv\,\frac{(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}+(\mathcal{B}_{\mathrm{c}}^{2})^{\prime\prime}\mu_{\mathrm{c}}-(\mathcal{B}_{\mathrm{c}}^{2})^{\prime\prime\prime}(1-\mu_{\mathrm{c}}^{2})/6}{(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}(1-\mu_{\mathrm{c}}^{2})}\,.$
(15)

If one substitutes the first Frobenius solution
$\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}$ for $\tilde{\psi}_{\mathrm{I}}$
into the integral expression (12b) to calculate the second solution up to its
second-order terms, the resulting expression certainly agrees with the second
Frobenius one $\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}$ up to the same
order (strictly speaking,
$\tilde{\psi}_{\mathrm{I\\!I}}\simeq-(\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}+\text{const.}\times\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})})/(1-\mu_{\mathrm{c}}^{2})m^{2}\alpha^{2}(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}$
around the latitude). Accordingly, we can associate
$\tilde{\psi}_{\mathrm{I}}$ with $\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}$,
and $\tilde{\psi}_{\mathrm{I\\!I}}$ with
$\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}$, in the vicinity of
$\mu=\mu_{\mathrm{c}}$. The series solutions (14) also justify (13b), since
the term which becomes the largest contributor to the integral is
$\int_{\mu_{\mathrm{c}}-\varDelta_{1}}^{\mu_{\mathrm{c}}+\varDelta_{2}}\ln|\mu-\mu_{\mathrm{c}}|\ud\mu$.
The comparison between linear combinations of these linearly independent
solutions and our numerical solutions will be shown in Section 3.2.

Necessary conditions for instability often receive interest from a lot of
hydrodynamicists. As will be proved in Appendix A, one of the so-called
semicircle theorems giving an eigenvalue bound in the current problem is
written as

$\left[\frac{\mathrm{Re}(\lambda)}{m}+\alpha^{2}\max\left(2\mathcal{B}\frac{\ud(\mathcal{B}\mu)}{\ud\mu}\right)\right]^{2}+\left(\frac{\mathrm{Im}(\lambda)}{m}\right)^{2}\,\leq\,\alpha^{4}\left[\max\left(2\mathcal{B}\frac{\ud(\mathcal{B}\mu)}{\ud\mu}\right)\right]^{2}\,-\,\alpha^{2}\min(\mathcal{B}^{2})\,,$
(16a) only if $\mathrm{Im}(\lambda)\neq 0$. Additionally, we also find another
bound for the case when $\mathrm{Im}(\lambda)\neq 0$ in the form
$-\frac{1}{2m(m+1)}\leq\frac{\mathrm{Re}(\lambda)}{m}\leq 0\,.$ (16b)

These relations mean that if unstable modes exist, they must propagate in the
retrograde direction, and the value
$\max\\{2\mathcal{B}[\ud(\mathcal{B}\mu)/\ud\mu]\\}$ concerning the gradient
of an imposed magnetic field must be positive. Hughes and Tobias (2001); Mak
et al. (2016), and Wang et al. (2022a, b) derived similar theorems for the MHD
or MHD shallow water systems with a background shear flow. The theorem (16a)
indicates that magnetic shear ascribed to the spherical geometry may have a
destabilising effect. We were not, however, able to find out unstable modes
that are likely to be physically meaningful when $\mathcal{B}=\mu$, as
described in Section 3.

###  2.1 Numerical method

A numerical method to solve our eigenvalue problem is described now. In this
article, we shall focus on the simplest equatorially-antisymmetric non-Malkus
field $\mathcal{B}=\mu$. This choice prompts us to utilise the associated
Legendre polynomial expansion, since (9) exactly becomes the associated
Legendre differential equation if $\mathcal{B}=1$ and recurrence formulae of
these polynomials are useful in the present situation. For a fixed $m$, the
polynomials $\mathrm{P}_{n}^{m}$ ($n\geq m$) constitute a basis for function
expansion in the Galerkin discretization on a spherical surface. Thus, we have

$\tilde{\psi}(\mu)\,\equiv\,\sum_{n=m}^{N_{\mathrm{t}}}\tilde{\psi}^{[n]}\mathcal{N}_{n}^{m}\mathrm{P}_{n}^{m}(\mu)\,,\qquad\frac{\operatorname{sgn}(\alpha)\tilde{a}}{\sqrt{\rho_{0}\mu_{\mathrm{m}}}}\,\equiv\,\sum_{n=m}^{N_{\mathrm{t}}}\tilde{a}^{[n]}\mathcal{N}_{n}^{m}\mathrm{P}_{n}^{m}\,,$
(17a) where $N_{\mathrm{t}}$ denotes the truncation degree and the normalising
factor is written as
$\mathcal{N}_{n}^{m}\,\equiv\,(-1)^{m}\sqrt{\frac{2n+1}{2}\frac{(n-m)!}{(n+m)!}}\,.$
(17b)

Assuming (17a) to be an approximate solution for (8) and using useful
relations in the forms

$\displaystyle\nabla_{\mathrm{h}}^{2}\mathrm{P}_{n}^{m}\,$
$\displaystyle=\,-n(n+1)\mathrm{P}_{n}^{m}\,,$ (18a)
$\displaystyle\mu\mathrm{P}_{n}^{m}\,$
$\displaystyle=\,\frac{n+m}{2n+1}\mathrm{P}_{n-1}^{m}\,+\,\frac{n-m+1}{2n+1}\mathrm{P}_{n+1}^{m}\,,$
(18b)

and the orthogonality relation
$\int_{-1}^{1}\mathrm{P}_{n}^{m}\mathrm{P}_{\nu}^{m}\ud\mu=\delta_{n,\nu}/\mathcal{N}_{n}^{m}\mathcal{N}_{\nu}^{m}$
with the Kronecker delta $\delta_{n,\nu}$, we get the following simultaneous
equations for all integers $n$ which satisfy $m\leq n\leq N_{\mathrm{t}}$:

$\displaystyle\left[\lambda+\frac{m}{n(n+1)}\right]\tilde{\psi}^{[n]}\,$
$\displaystyle=\,-m|\alpha|\left[\frac{(n-3)(n+2)}{n(n+1)}k_{n}^{m}\tilde{a}^{[n-1]}\,+\,\frac{(n-1)(n+4)}{n(n+1)}k_{n+1}^{m}\tilde{a}^{[n+1]}\right]\,,$
(19a) $\displaystyle\lambda\tilde{a}^{[n]}\,$
$\displaystyle=\,-m|\alpha|\left(k_{n}^{m}\tilde{\psi}^{[n-1]}\,+\,k_{n+1}^{m}\tilde{\psi}^{[n-1]}\right)\,,$
(19b) in which the number sequence
$k_{n}^{m}\,=\,\sqrt{\frac{(n-m)(n+m)}{(2n-1)(2n+1)}}\qquad(n=m,m+1,\ldots,N_{\mathrm{t}})\,,$
(19c)

is introduced for convenience.

The system of linear equations (19) is equivalent to the eigenvalue problem
for the corresponding $2(N_{\mathrm{t}}-m+1)\times 2(N_{\mathrm{t}}-m+1)$
matrix, and arrays of the expansion coefficients $\tilde{\psi}^{[n]}$ and
$\tilde{a}^{[n]}$ are its eigenvectors. We performed numerical calculations
solving this eigenvalue problem with our Python code which is based on the
numpy.linalg.eig function of the NumPy library. In these calculations, the
truncation number $N_{\mathrm{t}}$ was set to $2000$. As can be seen from
(19), the subset of $\tilde{\psi}^{[n]}$ with $n$ being odd numbers pertains
only to the subset of $\tilde{a}^{[n]}$ with $n$ being even, and the same is
true of the relationship between $\tilde{\psi}^{[n]}$ with $n$ even and
$\tilde{a}^{[n]}$ with $n$ odd. On the basis of this dichotomy, we refer to
eigenmodes for which $\tilde{\psi}^{[n]}$’s are non-zero only when $n-m$ is
even (in other words, $\tilde{a}^{[n]}$’s do not vanish only when $n-m$ is
odd, and $u_{\theta}$, $b_{\phi}$ are equatorially-symmetric and $u_{\phi}$,
$b_{\theta}$ are equatorially-antisymmetric) as the sinuous modes (cf.
Márquez-Artavia et al., 2017). Conversely, eigenmodes for which
$\tilde{\psi}^{[n]}$’s become non-zero only when $n-m$ is odd (or,
$u_{\theta}$ is antisymmetric and $u_{\phi}$ is symmetric about the equator)
are hereinafter referred to as the varicose modes.

The validity of eigenmodes obtained numerically is diagnosed from the aspect
of convergence. This is conducted by evaluating

$\sum_{n=m}^{\lfloor
N_{\mathrm{t}}/2\rfloor}\left|\tilde{\psi}^{[n]}\right|^{2}\,>\,10^{2}\sum_{n=\lfloor
N_{\mathrm{t}}/2\rfloor+1}^{N_{\mathrm{t}}}\left|\tilde{\psi}^{[n]}\right|^{2}\quad\text{and}\quad\sum_{n=m}^{\lfloor
N_{\mathrm{t}}/2\rfloor}\left|\tilde{a}^{[n]}\right|^{2}\,>\,10^{2}\sum_{n=\lfloor
N_{\mathrm{t}}/2\rfloor+1}^{N_{\mathrm{t}}}\left|\tilde{a}^{[n]}\right|^{2},$
(20)

where $\lfloor x\rfloor$ is the integer part of $x$. Only eigenmodes that pass
the screening by this validation will be studied and illustrated in the
results section.

Normalising the amplitudes of eigenmodes is valuable to those who want to
understand their characteristics by comparing physical quantities such as
energies. We employed the manner of normalization in which we let the mean
total energies $\mathrm{MKE}+\mathrm{MME}$ of perturbations be
$(\rho_{0}/8R_{0}^{2})\mathrm{e}^{2\mathrm{Im}(\omega)t}$, where the mean
kinetic and mean magnetic energies of an eigenmode are expressed as

$\displaystyle\mathrm{MKE}\,$
$\displaystyle\equiv\,\frac{1}{4\uppi}\int_{0}^{\uppi}\ud\theta\int_{0}^{2\uppi}\sin\theta\ud\phi\frac{\rho_{0}|\bm{u}_{1}|^{2}}{2}\,=\,\frac{\rho_{0}}{8R_{0}^{2}}\mathrm{e}^{2\mathrm{Im}(\omega)t}\sum_{n=m}^{N_{\mathrm{t}}}n(n+1)\left|\tilde{\psi}^{[n]}\right|^{2}\,,$
(21a) $\displaystyle\mathrm{MME}\,$
$\displaystyle\equiv\,\frac{1}{4\uppi}\int_{0}^{\uppi}\ud\theta\int_{0}^{2\uppi}\sin\theta\ud\phi\frac{|\bm{b}_{1}|^{2}}{2\mu_{\mathrm{m}}}\,=\,\frac{\rho_{0}}{8R_{0}^{2}}\mathrm{e}^{2\mathrm{Im}(\omega)t}\sum_{n=m}^{N_{\mathrm{t}}}n(n+1)\left|\tilde{a}^{[n]}\right|^{2}\,,$
(21b)

respectively, with $\bm{u}=\varepsilon\bm{u}_{1}+\mathrm{O}(\varepsilon^{2})$
and $\bm{B}=\bm{B}_{0}+\varepsilon\bm{b}_{1}+\mathrm{O}(\varepsilon^{2})$. The
energy partitioning between MKE and MME allows us to examine the force balance
of eigenmodes and to classify the types of waves. For instance, Figure 2
represents energy partitions for various eigenmodes under the Malkus field and
shows that MKEs dominate the mean total energies for fast MR waves, MMEs are
predominant over MKEs for slow MR waves, and Alfvén waves show almost
equipartition between the two for large $|\alpha|$. This normalization of
eigenvectors is applied to all the figures of profiles of eigenmodes displayed
in the following sections. The associated Legendre polynomials
$\mathrm{P}_{n}^{m}$ employed to construct the eigenfunctions from their
corresponding eigenvectors are given by the scipy.special.lpmv function of the
SciPy library.

Figure 2: Ratio of the mean kinetic energy MKE to the mean total energy
$\text{MKE}+\text{MME}$ against the absolute value $|\alpha|$ of the Lehnert
number when the zonal wavenumber $m=1$ and the Malkus field is imposed. This
plot is obtained from the relation
$\text{MKE}/(\text{MKE}+\text{MME})=\lambda^{2}/(\lambda^{2}+m^{2}\alpha^{2})$
in which the nondimensional angular frequency $\lambda$ is calculated by (1).
As $|\alpha|\to 0$, either MKE or MME approaches zero, depending on whether
the eigenmode is a slow or fast MR wave. The colours of the curves correspond
to those in Figure 1.

##  3 Numerical results

The results section starts by presenting the dispersion relation for our
current problem. Figure 3 shows the real parts of eigenfrequencies obtained
numerically when $m=1$ as functions of $|\alpha|$ (we find that their
imaginary parts vanish, that is $\mathrm{Re}(\lambda)=\lambda$, for all the
eigenmodes except for unreliable eigenmodes which will be slightly discussed
later). The figure has four panels. The left and right columns show the
retrograde and prograde modes, respectively, and the upper and lower rows show
the sinuous and varicose modes, respectively. Each colour in the scatter plots
represents the fraction of the mean kinetic energy within the mean total
energy of an eigenmode corresponding to a point on the diagrams. Note that we
used a nonlinear colour scale which is made by utilising the arctangent
function so as to highlight whether an eigenmode is similar to the Alfvén wave
($\text{MKE}\approx\text{MME}$; the colour of its marker is greenish) or not.
To make it easier to find markers of modes dissimilar to the Alfvén wave (we
chose the range $\text{MKE}<0.49$ or $0.51<\text{MKE}$), we furthermore set
the size of them to be larger than that of markers representing the Alfvén
wave. On the basis of the knowledge learned from Figure 2, even in the present
situation, we would be justified to think of an eigenmode whose marker in the
figure is coloured reddish as a mode similar to the fast MR wave and a bluish
one as a mode similar to the slow MR wave.

Figure 3: Dispersion diagrams when the zonal wavenumber $m=1$ and the simplest
equatorially-antisymmetric non-Malkus field $\mathcal{B}=\mu$ pervades the
system. The ordinates of each panel represent the real part
$\mathrm{Re}(\lambda)$ of the dimensionless angular frequency and the
abscissas are the absolute value $|\alpha|$ of the Lehnert number. Retrograde
modes ($\mathrm{Re}(\lambda)<0$) are shown in the left column, and the right
displays prograde modes ($\mathrm{Re}(\lambda)>0$). The upper and lower rows
illustrate sinuous and varicose modes, respectively. The colours of the
markers represent the ratio of the mean kinetic energy MKE to the mean total
energy MKE + MME.

In all the panels, we can observe bands crowded with eigenmodes lying just
below the lines $|\lambda|=m|\alpha|$. The kinetic and magnetic energies of
most eigenmodes in the bands are partitioned almost equally. We conjecture
that these bands should be identified with the continuous spectrum due to the
Alfvén resonance, which is hereinafter referred to as the Alfvén continuous
spectrum or Alfvén continuum, although our numerical method yields only
approximate discrete modes even when the system has a continuous spectrum.
Even though the spectrum ought to cover the expected range
$m|\alpha|\operatorname{min}(\mu)=-m|\alpha|\leq\lambda\leq
m|\alpha|=m|\alpha|\operatorname{max}(\mu)$ without any gaps as explained in
the previous sections, our eigenmodes satisfying (20) do not have eigenvalues
smaller than some levels in terms of absolute values. This is because, in
addition to the vertical axes of the panels being logarithmic scale, the
critical latitudes $\mu_{\mathrm{c}}=\pm\lambda/m|\alpha|$ get close to the
equator as $\lambda\to 0$ for a given $m$ and a given $|\alpha|$, a fine
structure in its eigenfunction appearing around the equator. It is necessary
to calculate with a higher truncation degree in order for (17a) to express
such a fine structure. Besides, small values of $|\alpha|$ ($\ll 1$) reduce
typical meridional wavelengths of perturbations, as will be shown in Figure 9.
It follows for the same reason as above that the bands are cut off below
certain values of $|\alpha|$. If numerical calculations with an infinite
degree were performed, obtained eigenvalues would cover the entire range below
the line $|\lambda|=m|\alpha|$ completely. We also conducted calculations with
truncation degrees somewhat lower than that of Figure 3, for example,
$N_{\mathrm{t}}=1000$ (not shown). The outlines of these dispersion diagrams
look almost unchanged aside from the difference in the widths of the bands;
the higher the truncation degree, the wider the band. In the lower panels of
Figure 3, some branches of the retrograde modes look like discrete eigenvalues
that lie below the bands. For these eigenmodes, MMEs are dominant. In
particular, the lowermost two eigenvalues, whose branches overlap (since they
are a complex conjugate pair, as mentioned in Section 2), in the lower left
panel have non-zero imaginary parts (not shown), and their values are
consistent with (16a) and (16b). However, we consider the branches including
the unstable modes as unreliable eigenvalues, or a part of the Alfvén
continuous modes, because the calculations (see also Figure 4) reveal that
these eigenvalues depend strongly upon $N_{\mathrm{t}}$ as opposed to normal
discrete modes (cf. Carpenter and Guha, 2019).

Figure 4: Dependence of the real parts $\mathrm{Re}(\lambda)$ of eigenvalues
on the truncation degree $N_{\mathrm{t}}$ when the zonal wavenumber $m=1$, the
absolute value of the Lehnert number $|\alpha|=1$, and the background field is
the simplest equatorially-antisymmetric non-Malkus one ($\mathcal{B}=\mu$).
The horizontal axis is the eigenvalue number when eigenvalues are arranged in
ascending order. The red open squares and the blue circles represent the case
when $N_{\mathrm{t}}=2000$ and $1999$, respectively.

Discrete branches equivalent to slow MR waves found in Figure 1 disappear from
Figure 3 as a result of the modification of the main field. Instead, the
markers dyed blue, for which the fractions of MMEs of their eigenmodes are
close to unity like slow MR waves, are distributed within the Alfvén
continuum. Therefore, we suggest that discrete modes of slow MR waves turn
into continuous ones under a non-Malkus field. This situation is similar to
that of equatorial Rossby waves in Taniguchi and Ishiwatari (2006), who
studied eigenmodes in a linear shear flow on an equatorial $\beta$-plane. An
alternative explanation may be their transformation into quasi-modes (Spencer
and Rasband, 1997; Schecter et al., 2000; Balmforth et al., 2001; Turner and
Gilbert, 2007; Wang et al., 2022a), or the Landau damping in a broad sense.
When embedded in a continuous spectrum due to the replacement of basic fields,
a discrete real eigenvalue (on the principle Riemann sheet) may change into a
complex pole on the next Riemann sheet (e.g. Crawford and Hislop, 1989).
Although this novel pole does not produce a true eigenmode, it plays a crucial
role in the time evolution of the system as a non-diffusive decaying
oscillation whose frequency is close to the original real eigenvalue. We
suspect that blue upward wedges with the approximate slope
$\lambda\propto|\alpha|^{2}$ at $|\alpha|\approx 10^{-2}$ and $\lambda\approx
10^{-4}$ in the right panels of Figure 3 are connected with quasi-modes
originating from slow MR waves. To confirm this, we are preparing a paper that
gives the result of a treatment finding quasi-modes, which is described in the
foregoing literature (e.g. Spencer and Rasband, 1997), different from our
approach.

Outside the continuous spectrum, that is, above the lines
$|\lambda|=m|\alpha|$ in the diagrams, fast MR waves remain discrete
eigenmodes even in the non-Malkus field. Their semi-analytical solutions can
be obtained from eigenvalues of the spheroidal differential equation to which
(9) is reduced when $\mathcal{B}=\mu$ and $m^{2}\alpha^{2}/\lambda^{2}$ is
small (see Appendix B). Additionally, we find that the lowest branch of the
$m=1$ sinuous modes of fast MR waves can penetrate the band of the continuous
spectrum without interaction (see the upper left panel of Figure 3). This is
because the Lorentz force does not act on this eigenmode. This mode is
explained in detail in Appendix C.

Figure 5 depicts the dispersion diagrams for $m=2$. They roughly epitomise the
diagrams when $m\geq 3$ (not shown). Their outlines do not change much from
those of Figure 3 with the exception of the absence of the branch of the fast
MR waves that penetrates the continuous spectrum. Again, more conspicuous blue
upward wedges exist around $10^{-2}\leq|\alpha|\leq 10^{-1}$ and
$10^{-4}\leq\lambda\leq 10^{-2}$ in the right panels of Figure 5.

Figure 5: Same as Figure 3, but for $m=2$. The two white asterisks of the
upper right panel correspond to the two eigenmodes depicted in Figure 7.

###  3.1 Eigenfunctions of the Alfvén continuous modes

We shall investigate the eigenfunctions of the Alfvén continuous modes in this
and the next subsections. Figure 6 shows typical structures of the
perturbations in the stream function $\psi_{1}$ and the magnetic vector
potential $a_{1}$ for the continuous modes. The dependences of their
amplitudes $\tilde{\psi}$ and $\tilde{a}$ on the colatitude are illustrated in
Figure 6(a), and are used for making the contour maps of Figure 6(b). From
these figures, we notice that spiky singular structures appear at the critical
latitudes of the eigenmode. As shown in Section 2, the eigenfunction has a
logarithmic singularity or a step function singularity or both. In Section
3.2, we will provide the ratios of their contributions to the eigenfunctions
around their corresponding critical latitudes. The eigenfunctions of
eigenmodes which are extracted from each of the two noticeable blue upward
wedges in the upper right panel of Figure 5 and which may have something to do
with quasi-modes stemming from slow MR waves are also plotted in Figure 7.

(a) Meridional structures.

(b) Global contour maps.

Figure 6: Eigenfunction of the sinuous mode with the nondimensional angular
frequency $\lambda\approx 0.05006$ when the zonal wavenumber $m=1$, the
absolute value of the Lehnert number $|\alpha|=0.1$, and the simplest
equatorially-antisymmetric non-Malkus field $\mathcal{B}=\mu$ is imposed. The
critical colatitudes $\theta_{\mathrm{c}}\approx 59.96^{\circ}$ and
$120.04^{\circ}$. (a) Amplitudes of the stream function $\tilde{\psi}$ (red
line) and the scaled magnetic vector potential
$\operatorname{sgn}(\alpha)\tilde{a}/\sqrt{\rho_{0}\mu_{\mathrm{m}}}$ (blue
line) as functions of the colatitude. (b) Contour maps of the stream function
$\psi_{1}$ (left panel) and the scaled magnetic vector potential
$\operatorname{sgn}(\alpha)a_{1}/\sqrt{\rho_{0}\mu_{\mathrm{m}}}$ (right
panel) in the Mollweide projection.

(a) Sinuous mode with $\lambda\approx 0.00012$.

(b) Sinuous mode with $\lambda\approx 0.00025$.

Figure 7: Same as Figure 6(a), but for $m=2$, $|\alpha|=0.013$.

To get the whole picture of the eigenfunctions of the continuous modes, we
exhibit those of all the obtained continuous modes. Figures 8 and 9 are the
heatmaps of the absolute values, or $|\tilde{\psi}|$ and $|\tilde{a}|$, of
their amplitudes as functions of $\lambda$ and the colatitude when
$|\alpha|=0.1$ and $0.01$, respectively, and $m=1$. The left columns in these
figures correspond to $|\tilde{\psi}|$, and $|\tilde{a}|$ is depicted in their
right ones. The sinuous and varicose modes are shown on the upper and lower
rows, respectively. The colour gets darker as the absolute value increases.
The maps indicate that the retrograde continuous modes are evanescent on the
polar side of the critical latitudes, while the prograde ones are evanescent
on the equatorial side. However, only the case in which $|\alpha|$ is
sufficiently smaller than unity displays this behaviour (see Figure 10).
Furthermore, the comparison between eigenmodes having the same value of
$\lambda/m|\alpha|$ in Figures 8 and 9 shows that the smaller the value of
$|\alpha|$ is, the smaller the typical north-south wavelengths of their
amplitudes. In Section 4, we will therefore examine the behaviour of wave
packets possessing large wavenumbers at a small $|\alpha|$, which is similar
to the Earth’s core conditions, on the basis of the ray theory and attempt to
get a better grasp of our numerical results.

Figure 8: Amplitudes of the eigenfunctions of all the obtained continuous
modes when the zonal wavenumber $m=1$, the absolute value of the Lehnert
number $|\alpha|=0.1$, and the basic field is the simplest equatorially-
antisymmetric non-Malkus one ($\mathcal{B}=\mu$). The four panels divide into
the stream function and the magnetic vector potential on the left and the
right columns, respectively, and sinuous and varicose modes on the upper and
lower rows, respectively. The vertical axes of each panel correspond to the
colatitude and the horizontal axes represent the nondimensional angular
frequency $\lambda$. The darker the shades of the colour, the higher the
absolute values $|\tilde{\psi}|$ and $|\tilde{a}|$ of the amplitudes of the
stream function and the magnetic vector potential.

Figure 9: Same as Figure 8, but for $|\alpha|=0.01$.

For moderate or large values of $|\alpha|$, the less striking difference
exists in the evanescent property between the retrograde and prograde modes
which possess the same absolute value $|\lambda|$ of their angular frequency.
This statement is based on Figure 10, which shows the heatmaps for $m=1$ and
$|\alpha|=1$, and other experiments with several values of $m$ and $|\alpha|$
(not shown). Therefore, the contrast in the property between the retrograde
and prograde modes as demonstrated in Figures 8 and 9 would be attributed to
the planetary $\beta$ effect, that is, the effect of rotation. The ray-tracing
analysis and the local dispersion relation which we are going to discuss in
Section 4 offer a similar explanation. In addition, the fast MR mode buried in
the continuous modes at $\lambda=-1/2$ is discernible in the upper panels of
Figure 10 (see also Appendix C). In fact, we tried to utilise these plots as a
means to discover buried discrete eigenmodes other than that, though no such
eigenmodes have been found.

Figure 10: Same as Figure 8, but for $m=1$ and $|\alpha|=1$.

The evanescent property can be judged by the function

$\mathcal{L}^{2}(\mu;m,\lambda)\,\equiv\,-m^{2}\,-\,\frac{m(1-\mu^{2})}{\varLambda}\left[\lambda+2m\alpha^{2}\mathcal{B}\frac{\ud(\mathcal{B}\mu)}{\ud\mu}\right]\,-\,\frac{1-\mu^{2}}{2\sqrt{\varLambda}}\frac{\ud}{\ud\mu}\left(\frac{1-\mu^{2}}{\sqrt{\varLambda}}\frac{\ud\varLambda}{\ud\mu}\right)\,,$
(22a) which appears in an alternative form of the differential equation (9)
$\frac{\ud}{\ud\mu}\left[(1-\mu^{2})\frac{\ud(\tilde{\psi}\sqrt{\varLambda})}{\ud\mu}\right]\,+\,\frac{\mathcal{L}^{2}}{1-\mu^{2}}(\tilde{\psi}\sqrt{\varLambda})\,=\,0\,.$
(22b)

The Mercator projection transformation $y=(1/2)\ln[(1+\mu)/(1-\mu)]$ yields a
differential equation of the harmonic oscillation,

$\frac{{}^{2}(\tilde{\psi}\sqrt{\varLambda})}{\ud
y^{2}}\,+\,\mathcal{L}^{2}(\tilde{\psi}\sqrt{\varLambda})\,=\,0\,,$ (23a)
which shows that the sign of $\mathcal{L}^{2}$ determines the evanescent
property of $\tilde{\psi}\sqrt{\varLambda}$ at the latitude. One can also
rewrite (22b) into the form
$\frac{{}^{2}[\tilde{\psi}\sqrt{(1-\mu^{2})\varLambda}]}{\ud\mu^{2}}\,+\,\frac{\mathcal{L}^{2}+1}{(1-\mu^{2})^{2}}[\tilde{\psi}\sqrt{(1-\mu^{2})\varLambda}]\,=\,0\,,$
(23b)

which likewise shows that the value of $\mathcal{L}^{2}$ informs us of the
evanescent property. The equivalent of this differential equation was derived
by Gilman and Fox (1999a), though the form of $\mathcal{L}^{2}$ differs from
theirs because their equation is based on a non-rotating frame; the partial
derivatives with respect to $t$ in our equations have to be replaced by
$(\upartial/\upartial t)+\varOmega_{0}(\upartial/\upartial\phi)$ in the non-
rotating frame. The left panel in Figure 11 shows contour plots of
$\mathcal{L}^{2}$ as a function of $\lambda$ and the colatitude when
$|\alpha|=0.01$ and $\mathcal{B}=\mu$. We observe that the area where
$\mathcal{L}^{2}>0$ (or $\mathcal{L}^{2}>-1$) in the panel certainly agrees
with the wavy regions in Figure 9. Now, we shall consider the case when
$|\alpha|$ is small. Noting that $\lambda=\mathrm{O}(|\alpha|)$ and
$|\varLambda|=\mathrm{O}(|\alpha|^{2})$ for continuous modes, we have

$\mathcal{L}^{2}\,=\,-\frac{m(1-\mu^{2})}{\varLambda}\lambda\,+\,\mathrm{O}(|\varLambda|^{-2}|\alpha|^{4})$
(24)

unless the latitudinal position $\mu$ is very close to a critical latitude
($|\varLambda|\gg\mathrm{O}(|\alpha|^{3})$), since
$(\ud\varLambda/\ud\mu)=\mathrm{O}(|\alpha|^{2})$ and
$({}^{2}\varLambda/\ud\mu^{2})=\mathrm{O}(|\alpha|^{2})$ as can be seen from
the definition of the function $\varLambda$. When $\mathcal{B}=\mu$ (and
$|\mu^{2}-\lambda^{2}/m^{2}\alpha^{2}|\gg\mathrm{O}(|\alpha|)$), the
oscillatory condition $\mathcal{L}^{2}>0$ on the equatorial side
($\mu^{2}<\lambda^{2}/m^{2}\alpha^{2}$) of the critical latitudes, therefore,
requires that $\lambda<0$ (the retrograde continuous modes), and
$\mathcal{L}^{2}>0$ on the polar side ($\mu^{2}>\lambda^{2}/m^{2}\alpha^{2}$)
for the prograde modes ($\lambda>0$). This explains the contrasting evanescent
behaviour between the retrograde and prograde continuous modes for a small
value of $|\alpha|$.

(a) $\mathcal{B}=\mu$.

(b) $\mathcal{B}=\mu\sqrt{1-\mu^{2}}$.

Figure 11: Dependence of the function $\mathcal{L}^{2}$ given by (22a) on the
nondimensional angular frequency $\lambda$ and the colatitude when the zonal
wavenumber $m=1$, the absolute value of the Lehnert number $|\alpha|=0.01$,
and the basic fields are (a) $\mathcal{B}=\mu$ and (b)
$\mathcal{B}=\mu\sqrt{1-\mu^{2}}$. The solid and dashed curves illustrate the
contour lines which correspond to $\mathcal{L}^{2}=0$ and $-1$, respectively.

###  3.2 Comparison of the numerical results with the Frobenius series
solutions

Here, we confirm that the eigenfunctions obtained numerically can be
approximated by linear combinations of the linearly independent Frobenius
series solutions (14). Let $\tilde{\psi}_{\mathrm{num}}$ be one of the
numerical stream functions, and we fit it into the form
$C_{\mathrm{I}}\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}+C_{\mathrm{I\\!I}}\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}$
by adjusting the coefficients $C_{\mathrm{I}}$ and $C_{\mathrm{I\\!I}}$ on
each side of the critical latitudes. The fitting procedure is as follows. The
colatitude ($0\leq\theta\leq\uppi$) is divided into $N_{\theta}$ points at
even intervals. On each side of the nearest point of a singular latitude, the
$N_{\mathrm{data}}$ points closest to the point are chosen among the
$N_{\theta}$ points for the fitting. For these points $\theta_{i}$ ($1\leq
i\leq N_{\mathrm{data}}$) either on the equatorial or polar sides, we may
write

$\frac{\tilde{\psi}_{\mathrm{num}}(\cos\theta_{i})}{\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}(\cos\theta_{i})}\,\approx\,C_{\mathrm{I\\!I}}\frac{\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}(\cos\theta_{i})}{\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}(\cos\theta_{i})}\,+\,C_{\mathrm{I}}\quad\text{and}\quad\frac{\tilde{\psi}_{\mathrm{num}}(\cos\theta_{i})}{\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}(\cos\theta_{i})}\,\approx\,C_{\mathrm{I}}\frac{\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}(\cos\theta_{i})}{\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}(\cos\theta_{i})}\,+\,C_{\mathrm{I\\!I}}\,.$
(25)

Now, $\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}$ and
$\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}$ are approximated by the second-
order Frobenius solutions with (2). We obtain a candidate value for each of
$C_{\mathrm{I}}$ and $C_{\mathrm{I\\!I}}$ from each equation, through least
squares fittings of (25) with the numpy.polyfit function of the NumPy library.
Thereby we have four candidate values for each of $C_{\mathrm{I}}$ and
$C_{\mathrm{I\\!I}}$ for one critical latitude, two from the equatorial side,
and two from the polar side. The upper panels of Figure 12 show a result when
we perform this procedure with $N_{\theta}=7201$, and $N_{\mathrm{data}}=200$
individually for the equatorial (red circles and solid lines) and the polar
(blue circles and dashed lines) sides. These fittings demonstrate that
$C_{\mathrm{I}}$ typically has different values
($C_{\mathrm{I}}^{(\mathrm{e})}$ and $C_{\mathrm{I}}^{(\mathrm{p})}$, say)
between the two sides of a critical latitude, whilst $C_{\mathrm{I\\!I}}$ has
the same value on both sides. Since $C_{\mathrm{I\\!I\\!I}}$ in (12a) can also
be written as
$\operatorname{sgn}(\mu_{\mathrm{c}})[C_{\mathrm{I}}^{(\mathrm{p})}-C_{\mathrm{I}}^{(\mathrm{e})}]$,
our numerical eigenmodes are consistent with the general results as described
in Section 2. We accordingly adopt the mean values of the two candidate values
for each of $C_{\mathrm{I}}^{(\mathrm{e})}$ and
$C_{\mathrm{I}}^{(\mathrm{p})}$ and four candidate ones of
$C_{\mathrm{I\\!I}}$
($C_{\mathrm{I\\!I}}^{(\mathrm{e})}=C_{\mathrm{I\\!I}}^{(\mathrm{p})}$) as
their definite values, which are used in the graph comparing
$\tilde{\psi}_{\mathrm{num}}$ with
$C_{\mathrm{I}}\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}+C_{\mathrm{I\\!I}}\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}$
(the lower panel of Figure 12).

Figure 12: Comparison between the numerical eigenfunction
$\tilde{\psi}_{\mathrm{num}}$ of the sinuous mode with the dimensionless
angular frequency $\lambda\approx 0.05006$ (the critical latitude
$\theta_{\mathrm{c}}=59.96^{\circ}$ in the north hemisphere) and linear
combinations
$C_{\mathrm{I}}\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}+C_{\mathrm{I\\!I}}\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}$
of the Frobenius series solutions (14) with (2) when the zonal wavenumber
$m=1$, the absolute value of the Lehnert number $|\alpha|=0.1$, and the
simplest equatorially-antisymmetric non-Malkus field $\mathcal{B}=\mu$
permeates the system. The undetermined coefficients $C_{\mathrm{I}}$ and
$C_{\mathrm{I\\!I}}$ are estimated from least squares fittings of
$\tilde{\psi}_{\mathrm{num}}/\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}=C_{\mathrm{I\\!I}}(\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}/\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})})+C_{\mathrm{I}}$
(upper left panel) and
$\tilde{\psi}_{\mathrm{num}}/\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}=C_{\mathrm{I}}(\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}/\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})})+C_{\mathrm{I\\!I}}$
(upper right panel) with $N_{\mathrm{data}}=200$ points on the equatorial (red
circles and solid lines) and the polar (blue circles and dashed lines) sides
of the critical latitude for $N_{\theta}=7201$. The lower panel shows
$\tilde{\psi}_{\mathrm{num}}$ (black curve) and
$C_{\mathrm{I}}\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}+C_{\mathrm{I\\!I}}\tilde{\psi}_{\mathrm{I\\!I}}^{(\mathrm{c})}$
with $C_{\mathrm{I}}$ and $C_{\mathrm{I\\!I}}$ determined by the fittings (red
and blue curves) as functions of the colatitude. The vertical grey shaded area
of this panel contains the $2N_{\mathrm{data}}=400$ points used in the
fittings.

The above procedure is also applied to all the continuous modes obtained from
our numerical calculations. Figure 13 depicts their values of
$\operatorname{sgn}(\mu_{\mathrm{c}})[C_{\mathrm{I}}^{(\mathrm{p})}-C_{\mathrm{I}}^{(\mathrm{e})}]$
($=C_{\mathrm{I\\!I\\!I}}$, red circles) and $C_{\mathrm{I\\!I}}$ (blue
circles) for $m=1$ and $|\alpha|=0.1$ in the left ordinates as functions of
$\lambda$. Their results for the sinuous and varicose modes are shown in the
left and right panels, respectively. Meanwhile, the right vertical axes of
these panels represent the numerical counterpart of (12d), which is written in
the present instance as

$I_{\mathrm{num}}\,\equiv\,-\operatorname{sgn}(\mu_{\mathrm{c}})\frac{C_{\mathrm{I}}^{(\mathrm{p})}-C_{\mathrm{I}}^{(\mathrm{e})}}{C_{\mathrm{I\\!I}}}(1-\mu_{\mathrm{c}}^{2})m^{2}\alpha^{2}(\mathcal{B}_{\mathrm{c}}^{2})^{\prime}\left[\tilde{\psi}_{\mathrm{I}}^{(\mathrm{c})}(\mu_{\mathrm{c}})\right]^{2}\,.$
(26)

This outcome demonstrates that the values of $I_{\mathrm{num}}$ appear to be
compatible with our expectation stated in Section 2; the values are arbitrary
numbers and adjusted to satisfy boundary conditions. In addition, we observe
that the value of
$\operatorname{sgn}(\mu_{\mathrm{c}})[C_{\mathrm{I}}^{(\mathrm{p})}-C_{\mathrm{I}}^{(\mathrm{e})}]$
vanishes at the extremum points of $C_{\mathrm{I\\!I}}$, and that
$\operatorname{sgn}(\mu_{\mathrm{c}})[C_{\mathrm{I}}^{(\mathrm{p})}-C_{\mathrm{I}}^{(\mathrm{e})}]$
has extremums at the zeros of $C_{\mathrm{I\\!I}}$. If discrete eigenmodes
without logarithmic and step function singularities are buried in the
continuum, the two values must simultaneously vanish.

Figure 13: Values of
$\operatorname{sgn}(\mu_{\mathrm{c}})[C_{\mathrm{I}}^{(\mathrm{p})}-C_{\mathrm{I}}^{(\mathrm{e})}]$
(red circles and the left vertical axes), $C_{\mathrm{I\\!I}}$ (blue circles
and the left vertical axes), and $I_{\mathrm{num}}$ defined as (26) (black
circles and the right vertical axes), relevant to the coefficients of the
Frobenius series solutions against the nondimensional angular frequency
$\lambda$. The case when the zonal wavenumber $m=1$, the absolute value of the
Lehnert number $|\alpha|=0.1$, and the simplest equatorially-antisymmetric
non-Malkus field $\mathcal{B}=\mu$ is imposed is shown in the left and the
right panels for sinuous and varicose modes, respectively. $N_{\theta}$ and
$N_{\mathrm{data}}$ are set to $7201$ and $200$, respectively, in the same way
as Figure 12.

##  4 Interpretation in terms of the ray theory and discussion

To further comprehend the continuous modes and their eigenfunctions, we here
reduce the system investigated so far to a more restricted situation when
$|\alpha|\ll 1$. We apply the ray theory, in which an inhomogeneous background
field varies with a much larger spatial scale than typical wavelengths. We
introduce local coordinates $(\varTheta,\varPhi)$ which suitably measure the
spatial scale size of the typical wavelength of a wave train. Introducing
small parameters helps incorporate such a setting into governing equations. In
Section 3.1, we found that typical meridional wavelengths decrease as the
value of $|\alpha|$ decreases. Thus, it would be reasonable to select
$|\alpha|$ as the parameter, if the value of $|\alpha|$ is sufficiently
smaller than unity. The local coordinates are then stretched in the forms

$\varTheta\,\equiv\,|\alpha|^{-1/2}\theta\,,\qquad\varPhi\,\equiv\,|\alpha|^{-1/2}\phi\,.$
(27a) The temporal scale of the wave period
($\lambda^{-1}=\mathrm{O}(|\alpha|^{-1/2})$) is similarly far from that of the
migration of the wave train. The new shrunk time $T$ useful in measuring the
latter is given by $T\,\equiv\,|\alpha|\tau\,.$ (27b)

The values of the exponents of $|\alpha|$ of the variables above are
determined in Appendix D.

We then introduce a locally defined wavenumber and angular frequency which
depend on the global coordinates $(\theta,\phi)$ and $T$, and subsequently
derive a local dispersion relation and ray-tracing equations, which predict
the movement of a wave packet. Their derivations are based on explanations in
standard textbooks on wave dynamics (e.g. Lighthill, 1978), and we explain
their details in Appendix D. Here we summarise the results. The expression of
perturbations of the stream function postulated in Section 2 is here rewritten
as
$\psi_{1}\equiv\mathrm{Re}[M(\phi,\theta,T)\mathrm{e}^{\mathrm{i}\varphi_{\mathrm{L}}(\varPhi,\varTheta,\tau)}]$,
where $M$ is the wave amplitude and $\varphi_{\mathrm{L}}$ is the phase of the
wave packet. With this ansatz, the local wavenumber and the local
nondimensional angular frequency are expressed as

$k(\phi,\theta,T)\,\equiv\,\frac{1}{\sin\theta}\frac{\upartial\varphi_{\mathrm{L}}}{\upartial\varPhi}\,,\qquad
l(\phi,\theta,T)\,\equiv\,\frac{\upartial\varphi_{\mathrm{L}}}{\upartial(-\varTheta)}\,,\qquad\lambda(\phi,\theta,T)\,\equiv\,-\frac{\upartial\varphi_{\mathrm{L}}}{\upartial\tau}\,.$
(28)

Note that $\varphi_{\mathrm{L}}$ depends on the local coordinates
$(\varPhi,\varTheta,\tau)$, while $M$, $k$, $l$ and $\lambda$ depend on the
global ones $(\phi,\theta,T)$. The local dispersion relation for our present
problem is obtained from the leading order terms in the governing equations
(7) in the form

$\mathcal{D}(\phi,\theta,T,k,l,H)\,\equiv\,H^{2}(k^{2}+l^{2})\,+\,Hk\sin\theta\,-\,k^{2}\mathcal{B}^{2}\sin^{2}\theta(k^{2}+l^{2})\,=\,0\,,$
(29)

in which $H\equiv|\alpha|^{-1/2}\lambda$ is the scaled nondimensional angular
frequency. Replacing $\sin\theta$ and $\mathcal{B}$ in (29) with constants,
one would have an equivalent to the nondimensional dispersion relation on a
middle latitude $\beta$ plane (Zaqarashvili et al., 2007). Note that we denote
$H=\mathcal{H}(\phi,\theta,k,l,T)$ as the solution of (29) for $H$. From (29),
the components of the nondimensional local group velocity
$\bm{c}_{\mathrm{g}}$ are given by

$\displaystyle\frac{c_{\mathrm{g},\phi}}{|\alpha|}\,$
$\displaystyle\equiv\,\frac{\upartial\mathcal{H}}{\upartial
k}\,=\,-\frac{(\upartial\mathcal{D}/\upartial
k)|_{H=\mathcal{H}}}{(\upartial\mathcal{D}/\upartial
H)|_{H=\mathcal{H}}}\,=\,\frac{2k\mathcal{B}^{2}\sin^{2}\theta(2k^{2}+l^{2})-\mathcal{H}(2k\mathcal{H}+\sin\theta)}{2\mathcal{H}(k^{2}+l^{2})+k\sin\theta}\,,$
(30a) $\displaystyle\frac{c_{\mathrm{g},-\theta}}{|\alpha|}\,$
$\displaystyle\equiv\,\frac{\upartial\mathcal{H}}{\upartial
l}\,=-\frac{(\upartial\mathcal{D}/\upartial
l)|_{H=\mathcal{H}}}{(\upartial\mathcal{D}/\upartial
H)|_{H=\mathcal{H}}}\,=\,-\frac{2l(\mathcal{H}^{2}-k^{2}\mathcal{B}^{2}\sin^{2}\theta)}{2\mathcal{H}(k^{2}+l^{2})+k\sin\theta}\,.$
(30b)

We eventually find the ray-tracing equations

$\sin\theta\frac{{}_{\mathrm{g}}\phi}{\ud
T}\,=\,\frac{c_{\mathrm{g},\phi}}{|\alpha|}\,,\qquad\frac{{}_{\mathrm{g}}(-\theta)}{\ud
T}\,=\,\frac{c_{\mathrm{g},-\theta}}{|\alpha|}\,,$ (31a) and
$\displaystyle\frac{{}_{\mathrm{g}}(k\sin\theta)}{\ud T}\,$
$\displaystyle=\,-\left(\frac{\upartial\mathcal{H}}{\upartial\phi}\right)_{k,l}\,=\,0\,,$
(31b) $\displaystyle\frac{{}_{\mathrm{g}}l}{\ud
T}\,+\,\frac{c_{\mathrm{g},\phi}}{|\alpha|}k\cot\theta\,$
$\displaystyle=\,-\left[\frac{\upartial\mathcal{H}}{\upartial(-\theta)}\right]_{k,l}\,=\,\frac{k^{2}[\ud(\mathcal{B}^{2}\sin^{2}\theta)/\ud\theta](k^{2}+l^{2})-\mathcal{H}k\cos\theta}{2\mathcal{H}(k^{2}+l^{2})+k\sin\theta}\,,$
(31c) $\displaystyle\frac{{}_{\mathrm{g}}H}{\ud T}\,$
$\displaystyle=\,\left(\frac{\upartial\mathcal{H}}{\upartial
T}\right)_{k,l}\,=\,0\,,$ (31d) where the material time derivative moving with
the local group velocity is $\frac{{}_{\mathrm{g}}}{\ud
T}\,\equiv\,\frac{\upartial}{\upartial
T}\,+\,\frac{c_{\mathrm{g},\phi}}{|\alpha|\sin\theta}\frac{\upartial}{\upartial\phi}\,+\,\frac{c_{\mathrm{g},-\theta}}{|\alpha|}\frac{\upartial}{\upartial(-\theta)}\,=\,\frac{\upartial}{\upartial
T}\,+\,\frac{\bm{c}_{\mathrm{g}}}{|\alpha|}\bm{\cdot}\bm{\nabla}_{\mathrm{G}}\,,$
(31e)

with
$\bm{\nabla}_{\mathrm{G}}\equiv(\hat{\bm{e}}_{\phi}/\sin\theta)(\upartial/\upartial\phi)+\hat{\bm{e}}_{-\theta}[\upartial/\upartial(-\theta)]$.
According to these equations, a wave train migrates with its group velocity
depending on its latitudinal position and its dominant local wavenumber, which
also varies with the colatitude $\theta$. Furthermore, (31) show that
$k\sin\theta$ and $H$ (or $\lambda$) are invariant along a ray trajectory, but
$l$ is not.

We conduct the numerical time integration of the ray-tracing equations (31)
with (30) for the movement of a wave packet originating at a given initial
position $(\theta,\phi)$ with a given initial local wavenumber $(k,l)$ and a
local dimensionless angular frequency $\lambda$ determined by the local
dispersion relation (29) (e.g. Teruya et al., 2022). In our code, the initial
local longitudinal wavenumber $k_{\mathrm{init}}$ is calculated from the
relation (29) with the scipy.optimize.fsolve function of the SciPy library
after one specifies the initial values of $l$, $H$, $\phi$ and $\theta$. The
succeeding time integration of (31) is based on an explicit Runge-Kutta method
of order 8 (the DOP853 algorithm in the `scipy.integrate.solve_ivp` function
of the SciPy library). This integration is conducted without explicitly using
(29) on the way, and the numerical errors in our calculations are monitored by
the value of the function $\mathcal{D}$ of (29). Their results are also
compared to those of Section 3.1 from the aspect of the evanescent property.

Before demonstrating their trajectories obtained numerically, we examine some
properties of the local dispersion relation. In the following preliminary
considerations, we assume that the physical variables satisfy the relation
$H=\mathcal{H}(\phi,\theta,k,l,T)$ at any time. Figures 14 and 15 show contour
plots of the scaled dimensionless angular frequency $H$ as a function of $k$
and $l$ for three different latitudes, when $\mathcal{B}=\cos\theta$ and
$\mathcal{B}=1$, respectively. For ease of understanding Figure 14, we first
explain Figure 15, which corresponds to Figure 1 for global modes. These
diagrams are calculated from the equation transformed from (29) in the form

$H\,=\,\frac{-k\sin\theta\pm|k|\sin\theta\sqrt{1+4\mathcal{B}^{2}(k^{2}+l^{2})^{2}}}{2(k^{2}+l^{2})}\,.$
(32a) In Figures 14 and 15, we take the plus of the plus-minus sign so that
$H$ should be positive. This means that the sign of $k$ signifies the
longitudinal direction of the phase velocity. The nearly vertical contour
lines for large absolute values of $l$ in Figure 15 correspond to the
relations describing the propagation properties of wave packets that belong to
prograde ($H/k>0$) and retrograde ($H/k<0$) Alfvén waves. Additionally, it can
be concluded that, in Figure 15, the circular contour lines which are tangent
to the line $k=0$ represent the dispersion relation for fast MR waves, and
that the slightly curved part of the nearly vertical contour lines near the
line $l=0$ on the half plane $H/k>0$ explains how slow MR waves propagate. The
similarities between Figure 15 and the left and middle panels of Figure 14
suggest that the same is true for the case where $\mathcal{B}=\cos\theta$.
However, for $\mathcal{B}=\cos\theta$, no branches of Alfvén and slow MR waves
exist at the equator ($\theta=90^{\circ}$) as shown in the right panel of
Figure 14, since the main field vanishes there. Note that the direction of the
gradient $(\upartial\mathcal{H}/\upartial k,\upartial\mathcal{H}/\upartial l)$
at a point $(k,l)$ in these plots is identical with that of the group velocity
of a wave packet whose dominant local wavenumber is $(k,l)$ at the colatitude
$\theta$. For a wave packet belonging to either Alfvén or slow MR waves, the
sign of the azimuthal component $c_{\mathrm{g},\phi}$ of its group velocity is
the same as that of the azimuthal component $|\alpha|(H/k)$ of its
nondimensional local phase velocity, whilst those of the meridional components
($c_{\mathrm{g},-\theta}$ and $|\alpha|(H/l)$) are opposite for the retrograde
Alfvén packet.

Figure 14 illustrates the remarkable feature that the north-south component
$c_{\mathrm{g},-\theta}$ of the group velocity vanishes at $l=0$ and
$l=\pm\infty$, as can also be seen from (30b). Then, the wave train can be
refracted at or absorbed into the latitude, heading only in the $\phi$
direction there (e.g. Acheson, 1972; McKenzie, 1973; Eltayeb, 1977; Eltayeb
and Mckenzie, 1977; Grimshaw, 1979). In particular, from (29) the latter
situation $l^{2}\to\infty$ with a reasonable condition $Hk\sin\theta\neq 0$
leads to the limit $H^{2}-k^{2}\mathcal{B}^{2}\sin^{2}\theta\to 0$, which
signifies that the latitude is a critical one. It follows that the nearly
vertical lines for Alfvén waves in Figure 14 should be linked to the Alfvén
continuous modes observed in the results in Section 3. Note that, though the
nearly vertical lines in Figure 15 are similar to those in Figure 14, the
Malkus field $\mathcal{B}=1$ does not yield any continuous modes since
$H^{2}-k^{2}\mathcal{B}^{2}\sin^{2}\theta$ is constant.

Figure 14: Local dispersion relation given by (32a) when the background field
is the simplest equatorially-antisymmetric non-Malkus one
($\mathcal{B}=\cos\theta$). In these panels, the scaled nondimensional angular
frequency $H=|\alpha|^{-1/2}\lambda$ is shown as a function of the local
wavenumber $(k,l)$ for the colatitudes $\theta=30^{\circ}$ (left panel),
$60^{\circ}$ (middle), and $90^{\circ}$ (right).

Figure 15: Same as Figure 14, but for $\mathcal{B}=1$.

The local dispersion relation (29) also allows one to understand the
evanescent property of waves from the sign of the squared local meridional
wavenumber $l^{2}$. This value can be calculated from

$l^{2}\,=\,-k^{2}\,-\,\frac{Hk\sin\theta}{H^{2}-k^{2}\mathcal{B}^{2}\sin^{2}\theta}\,.$
(32b)

Figure 16 contains contour plots of $l^{2}$ as a function of $k$ and $H$ for
three different latitudes in the case where $\mathcal{B}=\cos\theta$. Waves
can propagate only when their wavenumbers fall within the parameter domains
where $\l^{2}>0$ in these panels, and these regions are classified into three
groups. The two thin regions near the lines $H=\pm k|\mathcal{B}|\sin\theta$
in the left and middle panels correspond to the relations for prograde Alfvén
and slow MR waves ($H/k>0$) and retrograde Alfvén waves ($H/k<0$). Again, as
shown in the right panel, no areas for Alfvén and slow MR waves are found at
the equator ($\theta=90^{\circ}$), since the background field vanishes there.
The propagation properties of fast MR waves are represented as the domain near
the line $k=0$ on the half plane $H/k<0$. Note that again the curved lines
$l=0$ and the lines $H=\pm k|\mathcal{B}|\sin\theta$ show that wave packets
can be refracted at or absorbed into the latitude since the latter lines are
equivalent to the situation $l^{2}\to\infty$.

Figure 16: Local dispersion relation written as (32b) when the simplest
equatorially-antisymmetric non-Malkus field $\mathcal{B}=\cos\theta$ is
imposed. In these panels, the square $l^{2}$ of the local meridional
wavenumber is displayed as a function of the local longitudinal wavenumber $k$
and the scaled dimensionless angular frequency $H=|\alpha|^{-1/2}\lambda$ for
the colatitudes $\theta=30^{\circ}$ (left panel), $60^{\circ}$ (middle), and
$90^{\circ}$ (right).

Figure 16 are also helpful in the short-term prediction of the migration of a
wave packet. We will now consider the two cases: (i) when the wave packet
heads toward its corresponding critical latitude, and (ii) when the packet
proceeds in the direction away from the critical latitude. Since $k\sin\theta$
and $H$ remain constant during the movement of a wave packet, its migration in
the north-south direction can be converted into the movement of the point
$(k,H)$ in the horizontal direction of the panels in Figure 16 (unless the
outlines of their contour plots change significantly depending on the
latitude).

  1. 1.

When a wave train moves equatorward with $H/k$ and $l^{2}$ positive (then
$H/l<0$ in the north hemisphere from (30b) or Figure 14), $k$ decrease and the
point $(k,H)$ approaches the line $H=k|\mathcal{B}|\sin\theta$ from the right
on the plots of Figure 16. This means that the train belonging to either
prograde Alfvén or slow MR waves approaches its corresponding critical
latitude from the polar side and is refracted or absorbed there. If the train
went beyond the latitude, $l^{2}$ would become negative, hence evanescent
waves. Note that, although plots for the Malkus field $\mathcal{B}=1$, which
are illustrated in Figure 17, are similar to Figure 16 except for its right
panel ($\theta=90^{\circ}$), then the point $(k,H)$ never arrives near the
line $H=\pm k|\mathcal{B}|\sin\theta$ due to the constancy of $k\sin\theta$
and $H$ (unless the initial condition has already approximately satisfied this
equality). When $H/k$ is negative (under $\mathcal{B}=\cos\theta$), a wave
train that travels poleward ($H/l<0$ in the northern hemisphere) and that does
not belong to fast MR waves approaches its corresponding critical latitude,
because $k$ increases while it migrates and the point $(k,H)$ approaches the
line $H=-k|\mathcal{B}|\sin\theta$ from the left. It follows that the train
belonging to retrograde Alfvén waves approaches the critical latitude from the
equatorial side.

  2. 2.

A wave packet moving in the opposite direction from its corresponding critical
latitude is realised by a local meridional phase velocity $|\alpha|(H/l)$
which is opposite in sign to the case (i). In other words, we here focus on a
packet that travels poleward ($H/l>0$ in the northern hemisphere) with $H/k$
positive and one that moves equatorward ($H/l>0$ in the north hemisphere) with
$H/k$ negative. Then, the point $(k,H)$ approaches the curved line $l=0$,
resulting in its refraction or absorption. It can be concluded that the packet
that belongs to either prograde Alfvén or slow MR waves approaches the
latitude where $l=0$ from the equatorial side and that the packet belonging to
retrograde Alfvén waves approaches there from the polar side.

Figure 17: Same as Figure 16, but for $\mathcal{B}=1$.

Whether a wave packet is to be refracted at or absorbed into a latitude where
$l^{2}\to\infty$, or $H^{2}=k^{2}\mathcal{B}^{2}\sin^{2}\theta$, is considered
now. For a general profile of $\mathcal{B}$, near such a colatitude
$\theta=\theta_{\mathrm{c}}$, (32b) becomes

$l^{2}\,\simeq\,\frac{H/(k\sin\theta)}{(\ud\mathcal{B}^{2}/\ud\theta)|_{\theta=\theta_{\mathrm{c}}}(\theta-\theta_{\mathrm{c}})}\,,$
(33)

if $(\ud\mathcal{B}^{2}/\ud\theta)|_{\theta=\theta_{\mathrm{c}}}\neq 0$.
Specifically, when $\mathcal{B}=\cos\theta$, the oscillatory condition
$l^{2}>0$ requires that $(\theta-\theta_{\mathrm{c}})\cos\theta_{\mathrm{c}}$
has the sign opposite to $H/k$; the packet which belongs to either prograde
Alfvén or slow MR waves ($H/k>0$) approaches the critical colatitude
$\theta_{\mathrm{c}}$ from the polar side
($(\theta-\theta_{\mathrm{c}})\cos\theta_{\mathrm{c}}<0$), while the
retrograde Alfvén one ($H/k<0$) does from the equatorial side
($(\theta-\theta_{\mathrm{c}})\cos\theta_{\mathrm{c}}>0$). This is just a
mathematical paraphrase of the consideration of the ways that the packets
approach the critical latitudes in the previous paragraph. Since the term
$Hk\sin\theta$ in (29), which has the same sign as the numerator of (33),
represents the planetary $\beta$ effect, the aforementioned distinction
between the prograde and retrograde waves is caused by the $\beta$ effect.
Additionally, on using (30) and (33), one can obtain an asymptotic expression
of the group velocity

$\displaystyle\frac{c_{\mathrm{g},\phi}}{|\alpha|}\,$
$\displaystyle\simeq\,\frac{H}{k}\,=\,\mathrm{O}(|\theta-\theta_{\mathrm{c}}|^{0})\,,$
(34a) $\displaystyle\frac{c_{\mathrm{g},-\theta}}{|\alpha|}\,$
$\displaystyle\simeq\,\frac{(\ud\mathcal{B}^{2}/\ud\theta)|_{\theta=\theta_{\mathrm{c}}}(k\sin\theta)^{2}(\theta-\theta_{\mathrm{c}})}{Hl}\,=\,\mathrm{O}(|\theta-\theta_{\mathrm{c}}|^{3/2})\,.$
(34b)

This expression gives the travel time of the packet from a given latitude
$\theta$ in the vicinity of the critical one to the latter in the form

$\int_{\theta}^{\theta_{\mathrm{c}}}\frac{|\alpha|}{-c_{\mathrm{g},-\theta}(\theta_{*})}\ud\theta_{*}\,\simeq\,-\frac{H}{(\ud\mathcal{B}^{2}/\ud\theta)|_{\theta=\theta_{\mathrm{c}}}(k\sin\theta)^{2}}\int_{\theta}^{\theta_{\mathrm{c}}}\frac{l(\theta_{*})}{\theta_{*}-\theta_{\mathrm{c}}}\ud\theta_{*}\,=\,\mathrm{O}(|\theta-\theta_{\mathrm{c}}|^{-1/2})\,.$
(35)

Any packets therefore never reach their corresponding critical latitudes in a
finite time and are absorbed there. In contrast, wave packets are refracted at
the latitudes where $l$ vanishes, as will be explained in a similar fashion in
Appendix E. The latitudes are often referred to as “turning latitudes.”

At last, we shall demonstrate the ray trajectories obtained by the numerical
time integration of the ray-tracing equations, though rough ones can be
sketched even from the above examination. Figures 19 and 19 show two of the
trajectories for wave trains belonging to prograde and retrograde Alfvén
waves, respectively. Each of the trains is injected at the black asterisk in
each of their upper panels, its position evolving in accordance with (31a).
The colours in the trajectories represent their local wavenumbers, or the
directions of their local phase velocities, by hue (see their lower left
panels for the colour scale). In their lower right panels, the longitudes
$\phi$ of their positions at time $T$ are recorded in the same colouring
scheme as their upper panels. The numerical errors for the results shown in
Figures 19 and 19 are $\mathcal{D}\lesssim 10^{-3}$ throughout their numerical
integrations. The trajectories agree with the above predictions in terms of
the refraction at turning latitudes, the absorption into the critical ones,
and the incident directions to those. From the asymptotic expression (34a) of
the group velocity near a critical latitude, one obtains the period for a wave
packet to circle along the latitude around a sphere as
$\operatorname{sgn}(\varOmega_{0})T=|2\uppi k\sin\theta/H|$, which
approximates the periods read from the lower right panels of Figures 19 and
19.

(a) Initial local meridional wavenumber $l_{\mathrm{init}}=2$.

(b) Initial local meridional wavenumber $l_{\mathrm{init}}=-2$.

Figure 18: Ray trajectories in the Mollweide projection for wave packets
belonging to prograde Alfvén waves (the scaled zonal wavenumber
$k\sin\theta\approx 1.22338$, and the critical colatitude
$\theta_{\mathrm{c}}\approx 35.17^{\circ}$ in the north hemisphere) when the
simplest equatorially-antisymmetric non-Malkus field $\mathcal{B}=\cos\theta$
pervades the system (upper panels). The scaled dimensionless angular frequency
$H=|\alpha|^{-1/2}\lambda=1$. The black asterisks correspond to their initial
position
$(\theta_{\mathrm{init}},\phi_{\mathrm{init}})=(30^{\circ},0^{\circ})$. The
directions $(k,l)$ of their local phase velocities are represented by hue, the
colour scales of which are the left lower panels. The lower right panels show
the longitudes $\phi$ of their positions against the scaled nondimensional
time $T=|\alpha|\tau$ with the colour denoting the directions of their local
phase velocities. Note that both $\phi$ and $T$ are multiplied by
$\operatorname{sgn}(\varOmega_{0})$ so that the “time”
$\operatorname{sgn}(\varOmega_{0})T$ always become positive even when the
rotation rate $\varOmega_{0}$ of the sphere is negative.

(a) Initial local meridional wavenumber $l_{\mathrm{init}}=2$.

(b) Initial local meridional wavenumber $l_{\mathrm{init}}=-2$.

Figure 19: Same as Figure 19, but for retrograde Alfvén waves (the scaled
zonal wavenumber $k\sin\theta\approx-1.77036$, and the critical colatitude
$\theta_{\mathrm{c}}\approx 55.61^{\circ}$ in the north hemisphere). The
initial colatitude $\theta_{\mathrm{init}}=60^{\circ}$.

The facts that wave packets can not cross their corresponding critical
latitudes and that whether the packets approach there from the polar or
equatorial sides depends on the sign of $H/k$ are consistent with the features
observed from the numerical results for the continuous modes (Figures 8 and 9)
in the eigenvalue problem for global modes when $|\alpha|$ is small. Although
the spatial scale of waves focused on in this section is smaller than that of
global modes which were thematised in Section 3, this implies that, in the
case when $|\alpha|\ll 1$, ray trajectories for Alfvén waves enable one to
roughly predict the behaviours of the continuous modes without actually
solving the eigenvalue problem. For instance, four of the trajectories for
another equatorially-antisymmetric field $\mathcal{B}=\sin\theta\cos\theta$
are displayed in Figures 21 and 21. These figures indicate that, if
$\mathcal{B}=\sin\theta\cos\theta$, the prograde continuous modes should be
evanescent in the polar and equatorial regions, while the retrograde ones
should become evanescent in the mid-latitudes. The value of the function
$\mathcal{L}^{2}$ of (22a) also gives the information about the evanescent
property for global modes, as with the sign of $l^{2}$. This function
explicitly includes the effects of the gradient of the background field, parts
of which have been indirectly ignored in the derivation of the ray-tracing
equations (see Appendix D). However, its approximate formula (24) for a small
$|\alpha|$ agrees with (33). Figure 11(b) depicts the values of
$\mathcal{L}^{2}$ when $\mathcal{B}=\mu\sqrt{1-\mu^{2}}$, and our prediction
using the ray theory is consistent with this plot. Based on these facts, we
can consider that wave packets belonging to Alfvén waves pertain to the
continuous modes in Section 3 as expected.

(a) Initial local meridional wavenumber $l_{\mathrm{init}}=2$.

(b) Initial local meridional wavenumber $l_{\mathrm{init}}=-2$.

Figure 20: Same as Figure 19, but for $\mathcal{B}=\sin^{2}\theta\cos\theta$.
Ray trajectories for wave packets that belong to prograde Alfvén waves (the
scaled zonal wavenumber $k\sin\theta\approx 2.15606$, and the critical
colatitudes $\theta_{\mathrm{c}}\approx 34.03^{\circ},55.97^{\circ}$ in the
north hemisphere) are shown. The initial colatitude
$\theta_{\mathrm{init}}=45^{\circ}$.

(a) Scaled zonal wavenumber $k\sin\theta\approx-2.19789$, and the critical
colatitudes $\theta_{\mathrm{c}}\approx 32.75^{\circ},57.25^{\circ}$ in the
north hemisphere. The initial colatitude $\theta_{\mathrm{init}}=30^{\circ}$.

(b) Scaled zonal wavenumber $k\sin\theta\approx-2.04809$, and the critical
colatitudes $\theta_{\mathrm{c}}\approx 38.78^{\circ},51.22^{\circ}$ in the
north hemisphere. The initial colatitude $\theta_{\mathrm{init}}=60^{\circ}$.

Figure 21: Same as Figure 21, but for retrograde Alfvén waves. The initial
local meridional wavenumber $l_{\mathrm{init}}=2$.

Moreover, we here learn why discrete branches of slow MR waves disappear under
non-Malkus fields, as shown in the dispersion diagrams in Section 3, from the
ray-tracing approach. Figure 22 depicts one of the trajectories for a wave
packet that (at least initially) belongs to the slow MR wave, because the
condition

$H^{2}(k^{2}+l^{2})\,\ll\,Hk\sin\theta$ (36)

has been satisfied by its initial condition. This inequality is obtained from
the comparison between the first two terms in the local dispersion relation
(29) by analogy with slow MR waves in the case of the Malkus field. The
illustrated trajectory is similar to those of prograde Alfvén waves (see
Figure 19). Since the propagation properties for slow MR waves and those of
prograde Alfvén waves are continuous, as seen in Figures 14 and 16, wave
packets for slow MR waves transform into prograde Alfvén waves as they migrate
in an inhomogeneous background field, changing their dominant local meridional
wavenumber. The time evolution of $l$ can reverse the inequality sign of the
condition (36) as $l^{2}\to\infty$, when the Alfvén balance $H^{2}\approx
k^{2}\mathcal{B}^{2}\sin^{2}\theta$ is reached. From (32b), the latitudinal
variation of $l^{2}$ during the migration of a wave packet is written as

$\left(\frac{\upartial
l^{2}}{\upartial\theta}\right)_{k\sin\theta,H}\,=\,2k^{2}\cot\theta\,-\,\frac{Hk^{3}\sin^{3}\theta}{(H^{2}-k^{2}\mathcal{B}^{2}\sin^{2}\theta)^{2}}\frac{\ud\mathcal{B}^{2}}{\ud\theta}\,,$
(37)

and one finds that $(\upartial l^{2}/\upartial\theta)_{k\sin\theta,H}$ is
always positive (negative) in the northern (southern) hemisphere when
$\mathcal{B}=\cos\theta$ and $H/k>0$. Therefore, mode conversion from slow MR
into prograde Alfvén waves should have occurred between the initial colatitude
$\theta_{\mathrm{init}}$ and the critical colatitude $\theta_{\mathrm{c}}$ in
the example of Figure 22, since $0\leq l^{2}\leq l_{\mathrm{init}}^{2}$ within
the interval between $\theta_{\mathrm{init}}$ and the turning latitude.
Although mode conversions can cause the valve effect (e.g. Acheson, 1972;
McKenzie, 1973; Eltayeb, 1977; Grimshaw, 1979), the effect does not occur in
our system because there exists only one kind of critical latitudes: the
Alfvén resonance $H^{2}=k^{2}\mathcal{B}^{2}\sin^{2}\theta$. On the other
hand, wave trains belonging to fast MR waves, which have discrete branches
when $\lambda<-m|\alpha|$ for global modes, move back and forth between two
turning latitudes without transforming into Alfvén waves even though its
dominant local wavenumber evolves (see Appendix B).

Figure 22: Same as Figure 19, but for slow MR waves (the scaled zonal
wavenumber $k\sin\theta\approx 0.14633$, and the critical colatitude
$\theta_{\mathrm{c}}\approx 86.08^{\circ}$ in the north hemisphere). The
scaled nondimensional angular frequency $H=|\alpha|^{-1/2}\lambda=0.01$, the
initial local meridional wavenumber $l_{\mathrm{init}}=0.5$, and the initial
colatitude $\theta_{\mathrm{init}}=60^{\circ}$.

The invariants including the products of perturbations are useful for
understanding waves and their associated phenomena. As will be derived in
Appendix D, our approximation leads to a conservation law in the form

$\frac{\upartial}{\upartial T}\left(\frac{\upartial\mathcal{D}}{\upartial
H}|M|^{2}\right)\,+\,\bm{\nabla}_{\mathrm{G}}\bm{\cdot}\left(\frac{\bm{c}_{\mathrm{g}}}{|\alpha|}\frac{\upartial\mathcal{D}}{\upartial
H}|M|^{2}\right)\,=\,0\,.$ (38)

It follows from the constancy of $k\sin\theta$ and $H$ that the equation in
which $(\upartial\mathcal{D}/\upartial H)|M|^{2}$ is replaced by
$(k\sin\theta/H^{2})(\upartial\mathcal{D}/\upartial H)|M|^{2}$ is also
correct. Our subsequent paper will prove that the latter quantity is
equivalent to the pseudomomentum density (up to constant factor) for the 2D
ideal incompressible MHD system. If weak dissipations are introduced, a packet
that is related to the Alfvén continuous modes would attenuate near its
corresponding critical latitude, or within its corresponding thin inner
boundary layer (see Section 1), due to its long travel time (35). The fact
that the packet carries the pseudomomentum in line with (38) implies that the
mean flow would be accelerated there because the damping of waves can cause
angular momentum exchange between waves and a mean flow in accordance with the
wave-mean flow interaction theory (e.g. Bühler, 2009). This possibly may
induce nonlinear oscillations such as the quasi-biennial oscillation (QBO) in
the Earth’s equatorial stratosphere (e.g. Baldwin et al., 2001).

At the end of this section, we shall confirm the validity of the ray theory
when $l^{2}\to\infty$. The approximation requires that the spatial scales at
which the dominant local wavenumber $(k,l)$ and the amplitude $M$ of a wave
packet vary are sufficiently larger than its wavelength. This condition can be
written as

$\min\left(\left|\frac{1}{l}\frac{\upartial
l}{\upartial\theta}\right|^{-1},\left|\frac{1}{M}\frac{\upartial
M}{\upartial\theta}\right|^{-1}\right)\,\gg\,\frac{2\uppi}{|\alpha|^{-1/2}|l|}\,.\qquad$
(39)

To get an asymptotic expression of $|M|$ near its corresponding critical
latitude $\theta_{\mathrm{c}}$, we take advantage of the conservation law (38)
in the form

$\iint_{S_{\mathrm{g}}}\left(\frac{\upartial\mathcal{D}}{\upartial
H}|M|^{2}\right)\sin\theta\ud\theta\ud\phi\,=\,\text{const.}\,,$ (40)

where $S_{\mathrm{g}}$ is a region moving and deforming with the local group
velocity. Let $\theta_{1}(T)$ and $\theta_{2}(T)$ be the latitudinal
positions, at an arbitrary time $T$, of the rear and front of an isolated wave
packet, respectively. Since the local group velocity $\bm{c}_{\mathrm{g}}$ of
the packet depends only on the latitudinal position $\theta$, we have

$T\,=\,\int_{\theta_{1}(0)}^{\theta_{1}(T)}\frac{|\alpha|}{-c_{\mathrm{g},-\theta}}\ud\theta\,=\,\int_{\theta_{2}(0)}^{\theta_{2}(T)}\frac{|\alpha|}{-c_{\mathrm{g},-\theta}}\ud\theta\,.$
(41a) Thus, the latitudinal length $|\theta_{2}-\theta_{1}|$ of the packet
satisfies
$\int_{\theta_{1}(T)}^{\theta_{2}(T)}\frac{|\alpha|}{-c_{\mathrm{g},-\theta}}\ud\theta\,=\,\text{const.}\,.$
(41b)

Using these relations, we can estimate that $(\upartial\mathcal{D}/\upartial
H)|M|^{2}(\Delta\phi\sin\theta)\propto c_{\mathrm{g},-\theta}^{-1}$ with the
longitudinal length $\Delta\phi$ of the packet and achieve
$M=\mathrm{O}(|\theta-\theta_{\mathrm{c}}|^{-1/4})$. Note that the same
asymptotic expression as this is also obtained from the steady problem in the
last two paragraphs of Appendix D. Since
$l=\mathrm{O}(|\theta-\theta_{\mathrm{c}}|^{-1/2})$ (from (33)) and
$M=\mathrm{O}(|\theta-\theta_{\mathrm{c}}|^{-1/4})$, we can not properly
discuss the behaviour of a wave packet if it approaches within the distance
$|\theta-\theta_{\mathrm{c}}|=\mathrm{O}(|\alpha|)$ from its corresponding
critical latitudes. The case where $l=0$ is addressed in Appendix E.

##  5 Conclusion

In the present paper, we numerically scrutinised 2D ideal incompressible MHD
linear waves within a thin layer on a rotating sphere with latitudinally
varying toroidal magnetic fields
$B_{0\phi}=B_{0}\mathcal{B}(\cos\theta)\sin\theta$. From the eigenvalue
problem for the simplest equatorially-antisymmetric non-Malkus field
$\mathcal{B}=\cos\theta$, we did not find Alfvén and slow MR discrete branches
but a continuous spectrum. In particular, the slow MR waves turn into parts of
the Alfvén continuous modes owing to the imposition of the non-Malkus field.
The eigenfunctions of the continuous modes and their concomitant critical
latitudes were investigated in unprecedented detail and compared with the
Frobenius series solutions. The observed difference in the evanescent property
between the prograde and retrograde continuous modes results from the
planetary $\beta$ effect.

The theory of slowly varying wave packets, or the ray theory, in an
inhomogeneous magnetic field for a small absolute value of the Lehnert number
$\alpha$ implies that a wave packet that is related to the continuous spectrum
moves toward its corresponding critical latitude and is ultimately absorbed
there. The fact that whether the packet approaches the latitude from the polar
or equatorial sides depends on the sign of the azimuthal component
$|\alpha|(H/k)$ of its nondimensional local phase velocity is consistent with
the evanescent property of the global modes obtained from the eigenvalue
problem, and the planetary $\beta$ effect still causes this distinction of the
evanescent property. Additionally, this theory strongly corroborates the idea
that slow MR waves transform into Alfvén continuous modes under a non-Malkus
field. The novel conservation law (38) derived from this approximation will
provide insights into the interaction between waves and the mean flow and
magnetic fields, thereby into the weakly nonlinear evolution of the background
fields. Accordingly, our results could act as a stepping stone to a deeper
understanding of the dynamics of the outermost Earth’s core and the solar
tachocline.

Other considerations for our problem should be put off until future works.
Since slow MR waves occupy an important position as possible causes of
geomagnetic fluctuations, it is a pivotal issue whether such waves can remain
as discrete modes even under a non-Malkus field by recovering additional
effects which we have omitted here or not.

## References

  * Abramowitz and Stegun [1964] Abramowitz, M. and Stegun, I.A., Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Vol. 55, 1964 (US Government printing office).
  * Acheson [1972] Acheson, D.J., The critical level for hydromagnetic waves in a rotating fluid. Journal of Fluid Mechanics, 1972, 53, 401–415.
  * Adam [1986] Adam, J.A., Critical layer singularities and complex eigenvalues in some differential equations of mathematical physics. Physics Reports, 1986, 142, 263–356.
  * Arlt et al. [2007a] Arlt, R., Sule, A. and Filter, R., Stability of the solar tachocline with magnetic fields. Astronomische Nachrichten, 2007a, 328, 1142–1145.
  * Arlt et al. [2007b] Arlt, R., Sule, A. and Rüdiger, G., Stability of toroidal magnetic fields in the solar tachocline. Astronomy and Astrophysics, 2007b, 461, 295–301.
  * Baldwin et al. [2001] Baldwin, M.P., Gray, L.J., Dunkerton, T.J., Hamilton, K., Haynes, P.H., Randel, W.J., Holton, J.R., Alexander, M.J., Hirota, I., Horinouchi, T., Jones, D.B.A., Kinnersley, J.S., Marquardt, C., Sato, K. and Takahashi, M., The quasi-biennial oscillation. Reviews of Geophysics, 2001, 39, 179–229.
  * Balmforth et al. [2001] Balmforth, N.J., Llewellyn Smith, S.G. and Young, W.R., Disturbing vortices. Journal of Fluid Mechanics, 2001, 426, 95–133.
  * Balmforth and Morrison [1995a] Balmforth, N.J. and Morrison, P.J., Normal modes and continuous spectra. Annals of the New York Academy of Sciences, 1995a, 773, 80–94.
  * Balmforth and Morrison [1995b] Balmforth, N.J. and Morrison, P.J., Singular eigenfunctions for shearing fluids I. 1995b, Technical report, University of Texas at Austin. Institue for Fusion Studies.
  * Bardsley and Davidson [2017] Bardsley, O.P. and Davidson, P.A., The dispersion of magnetic-Coriolis waves in planetary cores. Geophysical Journal International, 2017, 210, 18–26.
  * Barré et al. [2015] Barré, J., Olivetti, A. and Yamaguchi, Y.Y., Landau damping and inhomogeneous reference states. Comptes Rendus Physique, 2015, 16, 723–728.
  * Barston [1964] Barston, E.M., Electrostatic oscillations in inhomogeneous cold plasmas. Annals of Physics, 1964, 29, 282–303.
  * Bouffard et al. [2022] Bouffard, M., Favier, B., Lecoanet, D. and Le Bars, M., Internal gravity waves in a stratified layer atop a convecting liquid core in a non-rotating spherical shell. Geophysical Journal International, 2022, 228, 337–354.
  * Boyd [1981] Boyd, J.P., Sturm-Liouville eigenproblems with an interior pole. Journal of Mathematical Physics, 1981, 22, 1575–1590.
  * Braginsky [1970] Braginsky, S.I., Torsional magnetohydrodynamic vibrations in the Earth’s core and variations in day length. Geomagn. Aeron., 1970, 10, 3–12.
  * Braginsky [1984] Braginsky, S.I., Short-period geomagnetic secular variation. Geophysical & Astrophysical Fluid Dynamics, 1984, 30, 1–78.
  * Braginsky [1993] Braginsky, S.I., MAC-oscillations of the hidden ocean of the core. Journal of geomagnetism and geoelectricity, 1993, 45, 1517–1538.
  * Braginsky [1998] Braginsky, S.I., Magnetic Rossby waves in the stratified ocean of the core, and topographic core-mantle coupling. Earth, Planets and Space, 1998, 50, 641–649.
  * Braginsky [1999] Braginsky, S.I., Dynamics of the stably stratified ocean at the top of the core. Physics of the Earth and Planetary Interiors, 1999, 111, 21–34.
  * Braun [1975] Braun, M., Differential Equations and Their Applications, 1975 (Springer-Verlag, New York).
  * Bretherton [1966] Bretherton, F.P., The propagation of groups of internal gravity waves in a shear flow. Quarterly Journal of the Royal Meteorological Society, 1966, 92, 466–480.
  * Briggs et al. [1970] Briggs, R.J., Daugherty, J.D. and Levy, R.H., Role of Landau damping in crossed-field electron beams and inviscid shear flow. The Physics of Fluids, 1970, 13, 421–432.
  * Brodholt and Badro [2017] Brodholt, J. and Badro, J., Composition of the low seismic velocity E’ layer at the top of Earth’s core. Geophysical Research Letters, 2017, 44, 8303–8310.
  * Buffett [2014] Buffett, B., Geomagnetic fluctuations reveal stable stratification at the top of the Earth’s core. Nature, 2014, 507, 484–487.
  * Buffett and Knezek [2018] Buffett, B. and Knezek, N., Stochastic generation of MAC waves and implications for convection in Earth’s core. Geophysical Journal International, 2018, 212, 1523–1535.
  * Buffett et al. [2016] Buffett, B., Knezek, N. and Holme, R., Evidence for MAC waves at the top of Earth’s core and implications for variations in length of day. Geophysical Journal International, 2016, 204, 1789–1800.
  * Buffett and Matsui [2019] Buffett, B. and Matsui, H., Equatorially trapped waves in Earth’s core. Geophysical Journal International, 2019, 218, 1210–1225.
  * Buffett and Seagle [2010] Buffett, B.A. and Seagle, C.T., Stratification of the top of the core due to chemical interactions with the mantle. Journal of Geophysical Research: Solid Earth, 2010, 115.
  * Bühler [2009] Bühler, O., Waves and Mean Flows, Cambridge Monographs on Mechanics 2009 (Cambridge University Press).
  * Cally [2001] Cally, P.S., Nonlinear evolution of 2D tachocline instabilities. Solar Physics, 2001, 199, 231–249.
  * Cally [2003] Cally, P.S., Three-dimensional magneto-shear instabilities in the solar tachocline. Monthly Notices of the Royal Astronomical Society, 2003, 339, 957–972.
  * Cally et al. [2004] Cally, P.S., Dikpati, M. and Gilman, P.A., The solar tachocline: Limiting magneto-tipping instabilities. Symposium - International Astronomical Union, 2004, 219, 541–545.
  * Cally et al. [2003] Cally, P.S., Dikpati, M. and Gilman, P.A., Clamshell and tipping instabilities in a two-dimensional magnetohydrodynamic tachocline. The Astrophysical Journal, 2003, 582, 1190–1205.
  * Cally et al. [2008] Cally, P.S., Dikpati, M. and Gilman, P.A., Three-dimensional magneto-shear instabilities in the solar tachocline - II. Axisymmetric case. Monthly Notices of the Royal Astronomical Society, 2008, 391, 891–900.
  * Carpenter and Guha [2019] Carpenter, J.R. and Guha, A., Instability of a smooth shear layer through wave interactions. Physics of Fluids, 2019, 31, 081701.
  * Case [1959] Case, K.M., Plasma oscillations. Annals of Physics, 1959, 7, 349–364.
  * Case [1960] Case, K.M., Stability of inviscid plane Couette flow. The Physics of Fluids, 1960, 3, 143–148.
  * Chi-Durán et al. [2021] Chi-Durán, R., Avery, M.S. and Buffett, B.A., Signatures of high-latitude waves in observations of geomagnetic acceleration. Geophysical Research Letters, 2021, 48, e2021GL094692.
  * Chulliat et al. [2015] Chulliat, A., Alken, P. and Maus, S., Fast equatorial waves propagating at the top of the Earth’s core. Geophysical Research Letters, 2015, 42, 3321–3329.
  * Couston et al. [2017] Couston, L.A., Lecoanet, D., Favier, B. and Le Bars, M., Dynamics of mixed convective–stably-stratified fluids. Physical Review Fluids, 2017, 2, 094804.
  * Crawford and Hislop [1989] Crawford, J.D. and Hislop, P.D., Application of the method of spectral deformation to the Vlasov-poisson system. Annals of Physics, 1989, 189, 265–317.
  * Davies et al. [2018] Davies, C.J., Pozzo, M., Gubbins, D. and Alfè, D., Partitioning of oxygen between ferropericlase and Earth’s liquid core. Geophysical Research Letters, 2018, 45, 6042–6050.
  * Dikpati et al. [2018a] Dikpati, M., Belucz, B., Gilman, P.A. and McIntosh, S.W., Phase speed of magnetized Rossby waves that cause solar seasons. The Astrophysical Journal, 2018a, 862, 159.
  * Dikpati et al. [2004] Dikpati, M., Cally, P.S. and Gilman, P.A., Linear analysis and nonlinear evolution of two-dimensional global magnetohydrodynamic instabilities in a diffusive tachocline. The Astrophysical Journal, 2004, 610, 597–615.
  * Dikpati et al. [2017] Dikpati, M., Cally, P.S., McIntosh, S.W. and Heifetz, E., The origin of the “seasons” in space weather. Scientific Reports, 2017, 7, 14750.
  * Dikpati and Gilman [1999] Dikpati, M. and Gilman, P.A., Joint instability of latitudinal differential rotation and concentrated toroidal fields below the solar convection zone. The Astrophysical Journal, 1999, 512, 417–441.
  * Dikpati et al. [2009] Dikpati, M., Gilman, P.A., Cally, P.S. and Miesch, M.S., Axisymmetric MHD instabilities in solar/stellar tachoclines. The Astrophysical Journal, 2009, 692, 1421–1431.
  * Dikpati et al. [2020] Dikpati, M., Gilman, P.A., Chatterjee, S., McIntosh, S.W. and Zaqarashvili, T.V., Physics of Magnetohydrodynamic Rossby Waves in the Sun. The Astrophysical Journal, 2020, 896, 141.
  * Dikpati et al. [2003] Dikpati, M., Gilman, P.A. and Rempel, M., Stability analysis of tachocline latitudinal differential rotation and coexisting toroidal band using a shallow-water model. The Astrophysical Journal, 2003, 596, 680–697.
  * Dikpati et al. [2018b] Dikpati, M., McIntosh, S.W., Bothun, G., Cally, P.S., Ghosh, S.S., Gilman, P.A. and Umurhan, O.M., Role of interaction between magnetic Rossby waves and tachocline differential rotation in producing solar seasons. The Astrophysical Journal, 2018b, 853, 144.
  * Drazin and Reid [1981] Drazin, P.G. and Reid, W.H., Hydrodynamic Stability, Cambridge Mathematical Library 1981 (Cambridge University Press).
  * Eltayeb [1977] Eltayeb, I.A., On linear wave motions in magnetic-velocity shears. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1977, 285, 607–636.
  * Eltayeb and Mckenzie [1977] Eltayeb, I.A. and Mckenzie, J.F., Propagation of hydromagnetic planetary waves on a beta-plane through magnetic and velocity shear. Journal of Fluid Mechanics, 1977, 81, 1–23.
  * Farrell [1982] Farrell, B.F., The initial growth of disturbances in a baroclinic flow. Journal of Atmospheric Sciences, 1982, 39, 1663 – 1686.
  * Finlay [2008] Finlay, C.C., Course 8 Waves in the presence of magnetic fields, rotation and convection. In Dynamos, Les Houches, edited by P. Cardin and L. Cugliandolo, Vol. 88 of Les Houches, pp. 403–450, 2008 (Elsevier: Amsterdam).
  * Gachechiladze et al. [2019] Gachechiladze, T., Zaqarashvili, T.V., Gurgenashvili, E., Ramishvili, G., Carbonell, M., Oliver, R. and Ballester, J.L., Magneto-Rossby waves in the solar tachocline and the annual variations in solar activity. The Astrophysical Journal, 2019, 874, 162.
  * Gastine et al. [2020] Gastine, T., Aubert, J. and Fournier, A., Dynamo-based limit to the extent of a stable layer atop Earth’s core. Geophysical Journal International, 2020, 222, 1433–1448.
  * Gillet et al. [2021] Gillet, N., Gerick, F., Angappan, R. and Jault, D., A dynamical prospective on interannual geomagnetic field changes. Surveys in Geophysics, 2021\.
  * Gillet et al. [2010] Gillet, N., Jault, D., Canet, E. and Fournier, A., Fast torsional waves and strong magnetic field within the Earth’s core. Nature, 2010, 465, 74–77.
  * Gilman [2000] Gilman, P.A., Magnetohydrodynamic “shallow water” equations for the solar tachocline. The Astrophysical Journal, 2000, 544, L79–L82.
  * Gilman and Dikpati [2000] Gilman, P.A. and Dikpati, M., Joint instability of latitudinal differential rotation and concentrated toroidal fields below the solar convection zone. II. Instability of narrow bands at all latitudes. The Astrophysical Journal, 2000, 528, 552–572.
  * Gilman and Dikpati [2002] Gilman, P.A. and Dikpati, M., Analysis of instability of latitudinal differential rotation and toroidal field in the solar tachocline using a magnetohydrodynamic shallow-water model. I. Instability for broad toroidal field profile. The Astrophysical Journal, 2002, 576, 1031–1047.
  * Gilman et al. [2007] Gilman, P.A., Dikpati, M. and Miesch, M.S., Global MHD instabilities in a three-dimensional thin-shell model of solar tachocline. The Astrophysical Journal Supplement Series, 2007, 170, 203–227.
  * Gilman and Fox [1997] Gilman, P.A. and Fox, P.A., Joint instability of latitudinal differential rotation and toroidal magnetic fields below the solar convection zone. The Astrophysical Journal, 1997, 484, 439–454.
  * Gilman and Fox [1999a] Gilman, P.A. and Fox, P.A., Joint instability of latitudinal differential rotation and toroidal magnetic fields below the solar convection zone. II. Instability for toroidal fields that have a node between the equator and pole. The Astrophysical Journal, 1999a, 510, 1018–1044.
  * Gilman and Fox [1999b] Gilman, P.A. and Fox, P.A., Joint instability of latitudinal differential rotation and toroidal magnetic fields below the solar convection zone. III. Unstable disturbance phenomenology and the solar cycle. The Astrophysical Journal, 1999b, 522, 1167–1189.
  * Gizon et al. [2020] Gizon, L., Fournier, D. and Albekioni, M., Effect of latitudinal differential rotation on solar Rossby waves: Critical layers, eigenfunctions, and momentum fluxes in the equatorial $\beta$ plane. Astronomy and Astrophysics, 2020, 642, A178.
  * Goedbloed and Poedts [2004] Goedbloed, J.P.H. and Poedts, S., Principles of Magnetohydrodynamics: With Applications to Laboratory and Astrophysical Plasmas, 2004 (Cambridge University Press).
  * Grimshaw [1979] Grimshaw, R., A general theory of critical level absorption and valve effects for linear wave propagation. Geophysical and Astrophysical Fluid Dynamics, 1979, 14, 303–326.
  * Gubbins and Davies [2013] Gubbins, D. and Davies, C.J., The stratified layer at the core-mantle boundary caused by barodiffusion of oxygen, sulphur and silicon. Physics of the Earth and Planetary Interiors, 2013, 215, 21–28.
  * Gubbins et al. [1982] Gubbins, D., Thomson, C.J. and Whaler, K.A., Stable regions in the Earth’s liquid core. Geophysical Journal International, 1982, 68, 241–251.
  * Hardy et al. [2020] Hardy, C.M., Livermore, P.W. and Niesen, J., Enhanced magnetic fields within a stratified layer. Geophysical Journal International, 2020, 222, 1686–1703.
  * Heifetz et al. [2020] Heifetz, E., Guha, A. and Carpenter, J.R., Wave interactions in neutrally stable shear layers: Regular and singular modes, and non-modal growth. Physics of Fluids, 2020, 32, 074106.
  * Helffrich and Kaneshima [2010] Helffrich, G. and Kaneshima, S., Outer-core compositional stratification from observed core wave speed profiles. Nature, 2010, 468, 807–810.
  * Heng and Spitkovsky [2009] Heng, K. and Spitkovsky, A., Magnetohydrodynamic shallow water waves: Linear analysis. The Astrophysical Journal, 2009, 703, 1819–1831.
  * Hide [1966] Hide, R., Free hydromagnetic oscillations of the Earth’s core and the theory of the geomagnetic secular variation. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1966, 259, 615–647.
  * Hollerbach and Cally [2009] Hollerbach, R. and Cally, P.S., Nonlinear evolution of axisymmetric twisted flux tubes in the solar tachocline. Solar Physics, 2009, 260, 251–260.
  * Hori et al. [2023] Hori, K., Nilsson, A. and Tobias, S.M., Waves in planetary dynamos. Reviews of Modern Plasma Physics, 2023, 7, 5.
  * Hughes and Tobias [2001] Hughes, D.W. and Tobias, S.M., On the instability of magnetohydrodynamic shear flows. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 2001, 457, 1365–1384.
  * Iga [1999] Iga, K., Critical layer instability as a resonance between a non-singular mode and continuous modes. Fluid Dynamics Research, 1999, 25, 63–86.
  * Iga [2013] Iga, K., Shear instability as a resonance between neutral waves hidden in a shear flow. Journal of Fluid Mechanics, 2013, 715, 452–476.
  * Irving et al. [2018] Irving, J.C.E., Cottaar, S. and Lekić, V., Seismically determined elastic parameters for Earth’s outer core. Science Advances, 2018, 4, eaar2538.
  * Jaupart and Buffett [2017] Jaupart, E. and Buffett, B., Generation of MAC waves by convection in Earth’s core. Geophysical Journal International, 2017, 209, 1326–1336.
  * Kaneshima [2018] Kaneshima, S., Array analyses of SmKS waves and the stratification of Earth’s outermost core. Physics of the Earth and Planetary Interiors, 2018, 276, 234–246.
  * Kitchatinov and Rüdiger [2008] Kitchatinov, L.L. and Rüdiger, G., Stability of toroidal magnetic fields in rotating stellar radiation zones. Astronomy and Astrophysics, 2008, 478, 1–8.
  * Knezek and Buffett [2018] Knezek, N. and Buffett, B., Influence of magnetic field configuration on magnetohydrodynamic waves in Earth’s core. Physics of the Earth and Planetary Interiors, 2018, 277, 1–9.
  * Landau [1946] Landau, L., On the vibration of the electronic plasma. Journal of Physics U.S.S.R., 1946, 10, 25.
  * Landeau et al. [2016] Landeau, M., Olson, P., Deguen, R. and Hirsh, B.H., Core merging and stratification following giant impact. Nature Geoscience, 2016, 9, 786–789.
  * Lighthill [1978] Lighthill, J., Waves in Fluids, 1978 (Cambridge University Press).
  * Longuet-Higgins [1968] Longuet-Higgins, M.S., The eigenfunctions of Laplace’s tidal equation over a sphere. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1968, 262, 511–607.
  * Mak et al. [2016] Mak, J., Griffiths, S.D. and Hughes, D.W., Shear flow instabilities in shallow-water magnetohydrodynamics. Journal of Fluid Mechanics, 2016, 788, 767–796.
  * Malkus [1967] Malkus, W.V.R., Hydromagnetic planetary waves. Journal of Fluid Mechanics, 1967, 28, 793–802.
  * Márquez-Artavia et al. [2017] Márquez-Artavia, X., Jones, C.A. and Tobias, S.M., Rotating magnetic shallow water waves and instabilities in a sphere. Geophysical and Astrophysical Fluid Dynamics, 2017, 111, 282–322.
  * Maslowe [1986] Maslowe, S.A., Critical layers in shear flows. Annual Review of Fluid Mechanics, 1986, 18, 405–432.
  * McKenzie [1973] McKenzie, J.F., On the existence of critical levels, with applications to hydromagnetic waves. Journal of Fluid Mechanics, 1973, 58, 709–726.
  * Miesch and Gilman [2004] Miesch, M.S. and Gilman, P.A., Thin-shell magnetohydrodynamic equations for the solar tachocline. Solar Physics, 2004, 220, 287–305.
  * Miesch et al. [2007] Miesch, M.S., Gilman, P.A. and Dikpati, M., Nonlinear evolution of global magnetoshear instabilities in a three-dimensional thin-shell model of the solar tachocline. The Astrophysical Journal Supplement Series, 2007, 168, 337–361.
  * Moffatt [1978] Moffatt, H.K., Magnetic field generation in electrically conducting fluids. Cambridge Monographs on Mechanics and Applied Mathematics, 1978.
  * Mound et al. [2019] Mound, J., Davies, C., Rost, S. and Aurnou, J., Regional stratification at the top of Earth’s core due to core–mantle boundary heat flux variations. Nature Geoscience, 2019, 12, 575–580.
  * Pozzo et al. [2012] Pozzo, M., Davies, C., Gubbins, D. and Alfè, D., Thermal and electrical conductivity of iron at Earth’s core conditions. Nature, 2012, 485, 355–358.
  * Raphaldini and Raupp [2020] Raphaldini, B. and Raupp, C.F.M., Nonlinear MHD Rossby wave interactions and persistent geomagnetic field structures. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2020, 476, 20200174.
  * Schecter et al. [2000] Schecter, D.A., Dubin, D.H.E., Cass, A.C., Driscoll, C.F., Lansky, I.M. and O’Neil, T.M., Inviscid damping of asymmetries on a two-dimensional vortex. Physics of Fluids, 2000, 12, 2397–2412.
  * Sedláček [1971] Sedláček, Z., Electrostatic oscillations in cold inhomogeneous plasma I. Differential equation approach. Journal of Plasma Physics, 1971, 5, 239–263.
  * Sharif and Jones [2005] Sharif, B.W. and Jones, C.A., Rotational and magnetic instability in the diffusive tachocline. Geophysical and Astrophysical Fluid Dynamics, 2005, 99, 493–511.
  * Shivamoggi [1992] Shivamoggi, B.K., Ideal and resistive magnetohydrodynamic modes. International Journal of Theoretical Physics, 1992, 31, 2121–2141.
  * Spencer and Rasband [1997] Spencer, R.L. and Rasband, S.N., Damped diocotron quasi-modes of non-neutral plasmas and inviscid fluids. Physics of Plasmas, 1997, 4, 53–60.
  * Spiegel and Zahn [1992] Spiegel, E.A. and Zahn, J.P., The solar tachocline. Astronomy and Astrophysics, 1992, 265, 106–114.
  * Steinolfson [1985] Steinolfson, R.S., Resistive wave dissipation on magnetic inhomogeneities Normal modes and phase mixing. The Astrophysical Journal, 1985, 295, 213–219.
  * Stewartson [1967] Stewartson, K., Slow oscillations of fluid in a rotating cavity in the presence of a toroidal magnetic field. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 1967, 299, 173–187.
  * Strogatz [2000] Strogatz, S.H., From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D: Nonlinear Phenomena, 2000, 143, 1–20.
  * Takehiro and Lister [2001] Takehiro, S. and Lister, J.R., Penetration of columnar convection into an outer stably stratified layer in rapidly rotating spherical fluid shells. Earth and Planetary Science Letters, 2001, 187, 357–366.
  * Taniguchi and Ishiwatari [2006] Taniguchi, H. and Ishiwatari, M., Physical interpretation of unstable modes of a linear shear flow in shallow water on an equatorial beta-plane. Journal of Fluid Mechanics, 2006, 567, 1–26.
  * Tataronis and Grossmann [1973] Tataronis, J. and Grossmann, W., Decay of MHD waves by phase mixing. Zeitschrift für Physik A Hadrons and nuclei, 1973, 261, 203–216.
  * Teruya et al. [2022] Teruya, A.S.W., Raphaldini, B. and Raupp, C.F.M., Ray tracing of MHD Rossby waves in the solar tachocline: Meridional propagation and implications for the solar magnetic activity. Frontiers in Astronomy and Space Sciences, 2022, 9.