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

,


Introduction
Geomagnetic and geodetic variations may be partially accounted for by magnetohydrodynamic (MHD) waves within the Earth's outer core (e.g.Triana et al., 2022;Gillet et al., 2022;Hori et al., 2023).For example, Hide (1966) suggested that the westward drift of the geomagnetic field may originate from slow magnetic Rossby (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 results in some inferences on the relevant physical quantities.Zatman and Bloxham (1997) accepted Braginsky's (1970) explanation and inferred that the magnitude of the cylindrical radial field is approximately  0.2 mT.This suggestion was challenged by Gillet et al. (2010), who associated the torsional oscillations (or torsional waves) with 6-year length-of-day signals to conclude that its magnitude is approximately 2 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 first considered, and the recent proposed physical mechanisms of stratification have become more sophisticated.For the thermal stratification, a subadiabatic temperature gradient due to the high thermal conductivity of the core (Pozzo et al., 2012;Zhang et al., 2022) has been invoked, and for the compositional stratification, barodiffusion and chemical interactions between the core and 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.The seismological evidence reported for such a layer has been a controversial subject (e.g.Helffrich and Kaneshima, 2010;Kaneshima, 2018;Irving et al., 2018;van Tent et al., 2020).Regional stratification -instead of global stratification -owing to the core-mantle boundary (CMB) heterogeneity was also proposed (Mound et al., 2019).The obscure properties of the layer could be inferred from the identification of the 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(Braginsky, , 1998(Braginsky, , 1999;;Buffett, 2014;Chulliat et al., 2015;Buffett et al., 2016;Knezek and Buffett, 2018;Chi-Durán et al., 2021).If a stratified layer exists at the top of the core, the 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;Couston et al., 2017;Buffett and Knezek, 2018;Bouffard et al., 2022).
Two-dimensional (2D) ideal incompressible MHD linear waves within a thin conducting fluid layer over a rotating sphere are focused on in the present study.A non-uniform toroidal (ϕ-directional, ϕ being the longitude) magnetic field is imposed on the fluid film as a background field, and the layer rotates rigidly with the sphere.To investigate the suggestion of Hide (1966), 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 (θ = 0 and π, where θ is the colatitude).One should therefore express the basic field as B 0ϕ (θ) = B 0 B sin θ with a constant B 0 and a bounded continuous function B(cos θ) (e.g.Gilman and Fox 1997; see Figure 1).The simplest profile B 0ϕ = B 0 sin θ in this expression, or B = 1, is equivalent to the wellknown Malkus field (e.g.Malkus 1967;Finlay 2008; see Figure 1(a)).This elementary 2D model is useful for qualitatively understanding 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 the main focus of this study.Under the Malkus field, the local Alfvén wave velocities are V Aϕ (θ) = (B 0 / √ ρ 0 µ m ) sin θ, which is proportional to sin θ, where ρ 0 is a uniform density, and µ m is the constant magnetic permeability.If local Alfvén waves were excited at the same longitude but at different latitudes simultaneously, they would travel with the same longitudinal angular velocity in a manner similar to rigid body rotation, as the length of the circle of colatitude θ is proportional to sin θ, as shown in Figure 1 (a).This indicates that the Malkus field at the surface of a sphere behaves like a uniform horizontal magnetic field and reduces our perturbation equations to a regular Sturm-Liouville problem because interior poles vanish (Boyd, 1981), as shown in Section 2. If the Malkus field permeates a rigidly rotating thin layer on top of a sphere having 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 λ = −m ± m 1 + 4α 2 n(n + 1)[n(n + 1) − 2] 2n(n + 1) , where λ ≡ ω/2Ω 0 is the angular frequency that is nondimensionalised by twice the constant (signed) rotation rate Ω 0 .In the above expression, m denotes the (positive) zonal wavenumber, α ≡ B 0 /2Ω 0 R 0 √ ρ 0 µ m is the (signed) Lehnert number, and n is the degree of the associated Legendre polynomial P m n (cos θ).Although the Lehnert number is often defined as Le ≡ 2|α| in geophysics literature, we follow the precedents mentioned above.The double sign in (1) corresponds to prograde and retrograde propagating Alfvén waves λ ≃ ±m|α| 1 − 2/n(n + 1) for a large value of |α|, and the two classes of MR waves as |α| → 0: the slow mode λ ≃ mα 2 [n(n + 1) − 2] peculiar to the MHD system and the fast mode λ ≃ −m/n(n + 1).The fast mode is closely related to the well-known (hydrodynamic) Rossby waves in meteorology and oceanography (e.g.Longuet-Higgins, 1968).Figure 2 shows that the exponent of |α| in a branch of λ, or (∂ log |λ|/∂ log |α|), can be used to distinguish the three types of waves derived from (1): Alfvén (λ ∝ ±|α| 1 ), fast MR (λ ∝ −|α| 0 ), and slow MR (λ ∝ |α| 2 ) waves.
The simplest equatorially antisymmetric non-Malkus field B 0ϕ = B 0 sin θ cos θ, or B = cos θ, should be more appropriate for the Earth's core than the pure Malkus field (Figure 1(b)).In the current study, we focus our attention on the case of this basic field profile because MHD waves under a non-Malkus field are poorly understood.It should be noted that the field is eastward in the northern hemisphere when B 0 > 0 and westward when B 0 < 0. It is important that |α| is significantly less 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 |α| ≈ 10 −4 with Ω 0 ≈ +0.729 × 10 −4 s −1 , R 0 ≈ 3480 km, ρ 0 ≈ 10 4 kg m −3 , and µ m ≈ 1.26 × 10 −6 H m −1 .
Non-Malkus fields complicate our linear wave problem because regular singular points appear in the governing equations, as demonstrated in Section 2. The singular colatitudes θ = θ c result from the Alfvén resonance, at which the azimuthal phase speed |ω|/(m/R 0 sin θ c ) of a wave is equal to the local Alfvén speed |V Aϕ (θ c )| = (|B 0 B(cos θ c )|/ √ ρ 0 µ m ) sin θ c (e.g.Uberoi, 1972;Goedbloed and Poedts, 2004).For a given angular frequency λ, such latitudes exist if min(B 2 ) ≤ λ 2 /m 2 α 2 ≤ max(B 2 ).For example, such singular latitudes ordinarily exist for low-frequency neutral modes under a non-Malkus configuration expressed by a function B that has at least one zero in 0 < θ < π.Even if such points exist for a given angular frequency λ, however, its eigenmode can be constructed when its eigenfunction is permitted to have singular profiles only at the points, and 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(B 2 ) ≤ λ 2 /m 2 α 2 ≤ max(B 2 ).
An eigenmode belonging to this continuous spectrum is hereinafter referred to as a "continuous eigenmode" or "continuous mode".In contrast to the eigenmode of a discrete spectrum, focusing on a single continuous eigenmode is futile since its eigenfunction is singular at its corresponding critical latitudes.Nevertheless, the integration (or superposition) with respect to λ in which the integrand is the eigenfunctions weighted by a coefficient depending only on λ, which may be considered as the inverse Fourier transform, can constitute a physically relevant solution.It should be noted that it is not always true that the Fourier transform of a well-behaved function is not pathological.Solutions composed only of continuous eigenmodes do not exhibit behaviour typical of collective oscillations by normal discrete eigenvalues but represent the transient growth of initial disturbances (Farrell, 1982) and non-diffusive attenuation (e.g.Adam, 1986).The latter involves phase mixing and 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).Related instability problems may pertain to the dynamics of the solar tachocline underlying the convection zone.Because unstable modes can be easily identified even in the presence of continuous spectra, this circumstance is significantly different from that of the study of linear waves.Gilman andFox (1997, 1999a,b) and the subsequent research (Dikpati and Gilman, 1999;Gilman and Dikpati, 2000;Zaqarashvili et al., 2010a) investigated 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 these studies, 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ϕ = (B 0 + B 1 cos 2 θ) sin θ cos θ (where 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 the latitudinal transport of the angular momentum (Spiegel and Zahn, 1992) and place an upper limit on the strength of a toroidal field stored within the layer through the ω-effect (Arlt et al., 2007a,b).The nonlinear evolution of the modes has been developed by Cally (2001) and Cally et al. (2003Cally et al. ( , 2004)), who reported the novel "clamshell" and "tipping" patterns.Furthermore, 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 three-dimensional (3D) thin shell model or "MHD hydrostatic primitive" equations (Miesch and Gilman, 2004) were newly proposed for evaluating the impacts of subadiabatic stratification and the weak vertical displacement of fluid particles.Gilman and Dikpati (2002) and Dikpati et al. (2003) demonstrated that the combinations of the differential rotation and toroidal fields in the shallow model easily become unstable again and that the growing perturbations have non-zero kinetic helicities, which are related to the α-effect (Moffatt, 1978).Furthermore, unstable MR waves in the layer may have caused some of the periodicities in the solar activity (Zaqarashvili et al. 2010a(Zaqarashvili et al. ,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(Dikpati et al., , 2018a,b) ,b) 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 were performed not only for non-axisymmetric unstable modes (Cally, 2003;Gilman et al., 2007;Kitchatinov and Rüdiger, 2008), but also for axisymmetric ones (Cally et al., 2008;Dikpati et al., 2009) in extensive studies of the 3D thin shell model.These studies include nonlinear simulations of these growth modes (Miesch et al., 2007;Hollerbach and Cally, 2009) and linear stability analyses that took into consideration the vertical profiles of the differential rotation and the background toroidal field (Arlt et al., 2007a,b).
Critical lines (or levels, latitudes, etc.) and their concomitant continuous modes have important bearings even on unstable modes, which are beyond the scope of this study.Although 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 with the continuous spectrum.This results in 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 intend to detail the critical latitudes and neutral continuous eigenmodes in advance to understand the linear stability in similar problems (e.g.Wang et al., 2022b), although our current problem does not comprise such unstable modes.
In contrast to the background magnetic field that we focus on in this paper, the dominant field may be radial in the outermost Earth's core owing to the small conductivity in the mantle (e.g.Knezek and Buffett, 2018).Nevertheless, critical latitudes can remain important.This is because the latitudes often reside in equations of similar eigenvalue problems when a radial background field depending on θ passes through a thin layer (e.g. the equations in Buffett and Matsui, 2019).Although the critical latitudes are sometimes neglected in approximations that are apparently appropriate for each target of study (Zaqarashvili, 2018;Buffett and Matsui, 2019), it is unclear whether the simplifications that drop the Alfvén resonance are always valid or not.Thus, it is necessary to scrutinise the influences of the critical latitudes on eigenmodes.Although a more realistic model should have a combination of a toroidal and poloidal imposed field with a finite layer thickness (e.g.Schmitt, 2010;Gerick et al., 2021;Luo et al., 2022), our simplified model would help us to understand such more complex, realistic ones.
Non-ideality and nonlinearity are also intimately related to the critical points.A small viscosity and magnetic diffusion can 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 intend to determine are functions on the real axis in the complex plane, they become non-singular with the removal of the singular points from the real axis.This appears to indicate that it is necessary to recover these ignored damping terms.However, very weak diffusions only result in thin boundary layers around the critical points (and near 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.The ideal model can, therefore, help us understand a weakly diffusive system (Rincon and Rieutord, 2003;Reese et al., 2004).Although decaying normal modes due to measurable dissipations (Steinolfson, 1985;Gizon et al., 2020) and nonlinear boundary layers (e.g.Tung, 1979;Maslowe, 1986) may also be 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 for seeking eigenmodes numerically based on the associated Legendre polynomial expansion in Section 2.4 after presenting some general results in Sections 2.1, 2.2, and 2.3.Section 3 presents numerical solutions for the case where the background field is expressed as the simplest equatorially antisymmetric non-Malkus field B = cos θ.In Section 3.1 and Appendix D, the structures of the obtained eigenfunctions are examined outside and near critical latitudes, respectively.In Section 4, we also conduct numerical integrations of ray-tracing equations at large wavenumbers (e.g.Bardsley and Davidson, 2017;Teruya et al., 2022) to interpret our eigenvalue problem from a different perspective.This approach tracks the paths of wave packets migrating with their group velocities.It should be noted that this is different from the calculations of the Stokes drift by Dikpati et al. (2020), who computed trajectories of fluid particles advected by an oscillatory flow induced by MR waves.Finally, we present the conclusions of this study in Section 5.

Mathematical formulation
We now begin with a description of the governing equations for the 2D ideal incompressible MHD in a thin layer on a rotating sphere.Here, we use the spherical coordinate system (r, θ, ϕ).Since the vertical flow should be weak within the thin stratified layer, we approximate that the velocity field u = (0, u θ , u ϕ ) relative to the rotating frame of reference is horizontal for simplicity.Furthermore, the magnetic field B = (0, B θ , B ϕ ) is assumed to have only horizontal components as our purpose is primarily to lay the foundation for more realistic models.To maintain consistency with this simplified setting, we select the azimuthal field B 0ϕ = B 0 B sin θ as the imposed field, as mentioned previously in the introduction section.In contrast to the toroidal field, radial ones normally bring vertical derivatives into their governing equations, and the operator is incompatible with a 2D model, which is computationally undemanding.Our 2D model can also be straightforwardly extended into the "MHD shallow water" system (Zaqarashvili et al., 2007(Zaqarashvili et al., , 2009;;Heng and Spitkovsky, 2009;Zaqarashvili et al., 2011;Márquez-Artavia et al., 2017).Moreover, Hardy et al. (2020) reported that small vertical motions under a strong stratification can enhance toroidal fields because of the "Malkus constraint".This suggestion supports our choice as a reasonable main field.
In the horizontal 2D model, we only have to consider the horizontal components of the momentum and induction equations.The horizontal momentum and the horizontal induction equations (e.g.Gilman and Fox, 1997) are written as with the solenoidal conditions where Π is the reduced pressure, ∇ H ≡ (ê θ /R 0 )(∂/∂θ)+(ê ϕ /R 0 sin θ)(∂/∂ϕ) is the horizontal nabla operator, and the material derivative is expressed as (D/Dt) ≡ (∂/∂t) + u • ∇ H .The 2D vorticity and 2D uncurled induction equations, obtained from (2), on the spherical surface r = R 0 are given as In the above equations, the radial components of the vorticity ζ and electrical current J are defined as where we introduce the stream function ψ and the magnetic vector potential A = (A, 0, 0).Owing to the solenoidal conditions of the fields u and B, the scalars ψ and A are related to the two vectors as Using these expressions, the governing equations (3) are written in terms of ψ and A (e.g.Raphaldini and Raupp, 2020) as where the operator J (f, g) for any two scalar functions f and g is defined as In what follows, waves having a small amplitude are considered in this system.We introduce a small positive parameter ε (≪ 1) representing the amplitude of the waves to rewrite ψ and A as The basic state is assumed to be rigid body rotation ψ 0 ≡ 0 with a latitude-dependent toroidal field B 0ϕ (θ).
With the expression of B 0ϕ given in the previous section, we obtain (dA 0 /dθ) = −R 0 B 0 B sin θ.It should be noted that a smooth background field at the poles necessitates that the function A 0 satisfies A 0 = ∞ n=1 T n P n , where T n (n = 1, 2, . ..) are constants, and P n is the Legendre polynomial of degree n. (6) is then written as where is the background electrical current.The normal mode approach is a common procedure in the study of linear waves.For a given azimuthal wavenumber m and an angular frequency ω -determined subsequently using a dispersion relation -which becomes (1) for the Malkus field B = 1 and is sought numerically in Section 3 when B = cos θ, we postulate that ψ 1 (θ, ϕ, t) ≡ Re[ ψ(µ; m, ω)e iφ(ϕ,t) ] = ψe iφ /2+c.c. and a 1 ≡ Re[ãe iφ ], where µ = cos θ, and φ ≡ mϕ−ωt is the phase of the waves.We use c.c. for the complex conjugate of the preceding terms.With the nondimensional time τ ≡ 2Ω 0 t, one also has φ = mϕ − λτ .Upon substituting this ansatz into (8), we obtain where It should be noted that sgn(α)ã/ √ ρ 0 µ m has the same dimensions as ψ.Hereinafter, our interest is limited to the case where m ̸ = 0.The above equations are easily transformed into a single ordinary differential equation of the form where the factor Λ(µ) ≡ λ 2 − m 2 α 2 B 2 is crucial to our problem.If there exist real values µ that satisfy Λ(µ) = 0 within the interval −1 < µ < 1, which are hereinafter denoted as µ c , those points are interior poles depending on λ.The poles yield continuous spectra, as stated in the introduction section.On the other hand, the Malkus field B = 1 obviously produces no singular points, with the exception of the endpoints µ = ±1.
On dividing (10) by the factor Λ (̸ = 0), one can obtain its dispersion relation (1) without any hurdles because (10) is then reduced to the associated Legendre differential equation.It should be noted that, if λ and ψ are an eigenvalue and its corresponding eigenfunction of (10), respectively, the same holds true for their complex conjugates.

Existence of continuous spectra
The existence of continuous spectra can be justified when the function Λ has zeros in −1 < µ < 1.Since (10) is a second-order differential equation, it should have two linearly independent solutions.Let ψI be a nonsingular solution, and let the other be tentatively expressed in the heuristic form ΨII = ψI µ f (µ * )dµ * with a function f (µ).Then, one has a non-zero Wronskian ψI (d ΨII /dµ) − (d ψI /dµ) ΨII = ψ2 I f , guaranteeing that ψI and ΨII are independent unless f = 0. On substituting the form of ΨII into (10), we obtain from which we obtain Λ(1 − µ 2 )f ψ2 I = const. .Therefore the solution should have the form with both C I and C II being constants.Let us examine how the integral changes if the interval of integration in (12) passes over zeros of Λ (we assume that ψ2 I ̸ = 0 within the interval).With an arbitrary starting point µ 0 of the integration interval (in the following equation, we suppose that µ 0 < min(µ, µ c ) and that at most one singular point exists between µ 0 and µ), we can express the second solution of the linearly independent set as the improper integral We may then deduce from (13) that ( 12) is also written as where the integer i is the index of the singular latitudes, and H(µ) is the step function.In addition, we have introduced a new second solution where P denotes the Cauchy principal value.In the vicinity of µ = µ c , we know that and so on.In this paper, we narrow down a target only to the case where (B 2 c ) ′ does not vanish.Using this series expansion, we obtain the magnitude C III,i of the discontinuities at the latitudes in the form with the integral I i given by This integral can be an arbitrary number, depending on how we take the limit ∆ 1,i , ∆ 2,i → +0.In fact, Van Kampen (1955) suggested that, in the problem of plasma oscillations, the counterpart of I may be considered as an arbitrary parameter, which can be determined by a normalisation condition of the counterpart of ψ, or the distribution function of plasma.As observed in (14a), the existence of more than two linearly independent solutions despite (10) being a second-order differential equation can result in 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 (10) in the vicinity of the singular point.On integrating (10) with respect to µ over the narrow range sandwiched between µ c − ∆ 1 and provided that the condition lim is fulfilled.The fact that ψI is surely a solution for (10) and that ψI H(µ − µ c ) always satisfies (15a) shows that the third term in (14a) is a weak solution for (10), and C III should then be an undetermined parameter.

Frobenius method
An alternative survey of the structure of the solution at the critical latitudes is conducted using the Frobenius method, which is a standard approach for providing a power series solution about a regular singular point of an ordinary differential equation (e.g.Braun, 1975).Let us suppose that a power series of the form ) is a solution for (10) around a critical latitude µ c , where ϱ is a number to be determined by substituting the form into basic equations.The equation that determines ϱ is referred to as the indicial equation.On substituting this assumption for ψ in (10), we obtain the equation ϱ 2 = 0 from its leading order term.As a result, two linearly independent solutions near the latitude are obtained on the basis of the method for dealing with a repeated root.The method is based on the fact that, for the temporary expression with the sequence ȃk (a 0 , ϱ) (k ≥ 1) selected to satisfy (10) apart from the term proportional to (µ − µ c ) ϱ−1 , we can rewrite the linear differential equation ( 10) as where M is the linear differential operator on the left side of (10).On setting ϱ = 0 and a 0 = 1 in (16b) and its ϱ derivative, we obtain ψ(c) where we used the formula It should be noted that (17b) has a logarithmic singularity.Some of the expansion coefficients are determined from ȃk (a 0 , ϱ) with slightly tedious but standard manipulations such as

and
If we substitute the first Frobenius solution ψ(c) I for ψI into the integral expression (14b) to calculate the second solution up to its second-order terms, the resulting expression certainly agrees with the second Frobenius solution ψ(c) II up to the same order (strictly speaking, ψII ≃ −( ) ′ around the latitude).Accordingly, we can associate ψI with ψ(c) I , and ψII with ψ(c) II , in the vicinity of µ = µ c .The series solutions (17) also justify (15b), since the term that becomes the largest contributor to the integral is The comparison between linear combinations of these linearly independent solutions and our numerical solutions is presented in Appendix D.

Necessary conditions for instability
The necessary conditions for instability often receive interest from many hydrodynamicists.We prove an inequality that gives an eigenvalue bound, of the form of the semicircle theorem, in Appendix A. It is written as only if Im(λ) ̸ = 0.Moreover, we also find another bound for the case where Im(λ) ̸ = 0 in the form These relations indicate that if unstable modes exist, they must propagate in the retrograde direction, and the value max{2B[d(Bµ)/dµ]} regarding the gradient of an imposed magnetic field must be positive.Hughes and Tobias (2001); Mak et al. (2016), andWang et al. (2022a,b) derived similar theorems for the MHD or MHD shallow water systems with a background shear flow.The theorem (19a) indicates that the magnetic shear ascribed to the spherical geometry may have a destabilising effect.We were not, however, able to identify unstable modes that are likely to be physically meaningful when B = µ, as described in Section 3.

Numerical method
A numerical method for solving our eigenvalue problem is described here.In this article, we focus on the simplest equatorially antisymmetric non-Malkus field B = µ.This choice prompts us to use the associated Legendre polynomial expansion, as (10) exactly becomes the associated Legendre differential equation if B = 1, and the recurrence formulae of these polynomials are useful in the present situation.For a fixed m, the polynomials P m n (n ≥ m) constitute a basis for function expansion in the Galerkin discretisation on a spherical surface.Thus, we obtain where N t denotes the truncation degree, and the normalising factor is written as On assuming (20a) to be an approximate solution for (9) and using useful relations of the forms and the orthogonality relation ν with the Kronecker delta δ n,ν , we obtain the following simultaneous equations for all the integers n that satisfy m ≤ n ≤ N t : where the number sequence is introduced for convenience.
The system of linear equations ( 22) is equivalent to the eigenvalue problem for the corresponding 2(N t − m + 1) × 2(N t − m + 1) matrix, and the arrays of the expansion coefficients ψ[n] and ã[n] are its eigenvectors.We performed numerical calculations for solving this eigenvalue problem using our Python code, which is based on the numpy.linalg.eigfunction of the NumPy library.In these calculations, the truncation number N t was set as 2000.As can be observed from ( 22), the subset of ψ[n] with n being odd numbers pertains only to the subset of ã[n] with n being even, and the same is true of the relationship between ψ[n] with n being even and ã[n] with n being odd.Based on this dichotomy, we refer to eigenmodes for which ψ[n] 's are non-zero only when n − m is even (i.e., ã[n] 's do not vanish only when n − m is odd, and u θ and b ϕ are equatorially symmetric, and u ϕ and b θ are equatorially antisymmetric) as the sinuous modes (cf.Márquez-Artavia et al., 2017).Conversely, eigenmodes for which ψ[n] 's become non-zero only when n − m is odd (or u θ is antisymmetric, and u ϕ 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 realised by evaluating where ⌊x⌋ is the integer part of x.Only eigenmodes that pass the screening via this validation are studied and illustrated in the results section.
The normalisation of the amplitudes of eigenmodes is valuable for understanding their characteristics by comparing physical quantities such as energies.The normalisation was realised by letting the mean total energies MKE + MME of the perturbations be (ρ 0 /8R 2 0 )e 2Im(ω)t , where the mean kinetic and mean magnetic energies of an eigenmode are expressed as respectively, with u = εu 1 + O(ε 2 ) and B = B 0 + εb 1 + O(ε 2 ).The energy partitioning between MKE and MME allows us to examine the force balance of the eigenmodes and to classify the types of waves.For instance, Figure 3 presents the energy partitions for various eigenmodes under the Malkus field and shows that , where the nondimensional angular frequency λ is calculated using (1).As |α| → 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 2.
MKEs dominate the mean total energies in the case of the fast MR waves, MMEs are predominant over MKEs in the case of the slow MR waves, and Alfvén waves exhibit almost equipartition between the two for a large |α|.This normalisation of eigenvectors is applied to all the figures of the profiles of eigenmodes displayed in the following sections.The associated Legendre polynomials P m n employed to construct the eigenfunctions from their corresponding eigenvectors are given by the scipy.special.lpmvfunction of the SciPy library.

Numerical results
We start this section by presenting the dispersion relation for our current problem.Figure 4 shows the real parts of the eigenfrequencies obtained numerically when m = 1 as functions of |α| (we find that their imaginary parts vanish, that is Re(λ) = λ, for all the eigenmodes except for unreliable eigenmodes, which are discussed briefly later).The figure has four panels.The left and right columns present the retrograde and prograde modes, respectively, and the upper and lower rows present the sinuous and varicose modes, respectively.It should be noted that, in the sinuous modes, ψ is symmetric and ã is antisymmetric about the equator, while ψ is equatorially antisymmetric and ã is equatorially symmetric in the varicose modes, as defined in Section 2.4.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.It should be noted that we used a nonlinear colour scale that is created using the arctangent function to highlight whether an eigenmode is similar to the Alfvén wave (MKE ≈ 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 select the range MKE < 0.49 or 0.51 < MKE), we furthermore set their size to be larger than that of the markers representing the Alfvén wave.Based on the information presented in Figure 3, even in the present situation, we would be justified in thinking of an eigenmode having a reddish marker as a mode similar to the fast MR wave and a bluish one as a mode similar to the slow MR wave.In all the panels, we can observe the bands crowded with eigenmodes just below the lines |λ| = m|α|.The kinetic and magnetic energies of the majority of the 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 should encompass the expected range m|α| min(µ) = −m|α| ≤ λ ≤ m|α| = m|α| max(µ)  without any gaps as explained in the previous sections, our eigenmodes satisfying (23) do not have eigenvalues smaller than some levels in terms of absolute values.This is because a fine structure in its eigenfunction appears around the equator, as the critical latitudes µ c = ±λ/m|α| approach the equator as λ → 0 for a given m and a given |α|, in addition to the vertical axes of the panels presented on a logarithmic scale.It is necessary to perform the calculation with a higher truncation degree for (20a) to express such a fine structure.Moreover, small values of |α| (≪ 1) reduce the typical meridional wavelengths of perturbations, as shown in Figure 10.For the aforementioned reason, the bands are cut off below certain values of |α|.If numerical calculations were performed with an infinite degree, the obtained eigenvalues would cover the entire range below the line |λ| = m|α|.We also performed calculations with truncation degrees less than that shown in Figure 4, for example, N t = 1000 (not shown).The outlines of these dispersion diagrams appear 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 4, some branches of the varicose modes look like discrete eigenvalues that lie below the bands.For these eigenmodes, MMEs are dominant.In particular, the lowermost two eigenvalues, which are represented by overlapping curves (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 (19a) and (19b).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 Figure 5) reveal that these eigenvalues depend strongly upon N t as opposed to normal discrete modes (cf.Carpenter and Guha, 2019).These Alfvén continuous modes for which MMEs dominate over MKEs, which are found only in the m = 1 varicose modes (see also Figure 6), could possibly be related to the "clamshell instability" (e.g.Wang et al., 2022b, see also the last paragraph of Appendix C).
Discrete branches equivalent to slow MR waves found in Figure 1 disappear from Figure 4 as a result of the modification of the main field.Instead, the blue markers, for which the fractions of MMEs of their eigenmodes are close to unity as in the case of slow MR waves, are distributed within the Alfvén continuum.These markers form blue upward wedges having the approximate slope λ ∝ |α| 2 at |α| ≈ 10 −2 and λ ≈ 10 −4 in the right panels of Figure 4. Therefore, we suggest that the discrete modes of slow MR waves turn into continuous ones under a non-Malkus field.This is similar to the case of equatorial Rossby waves in the study of Taniguchi and Ishiwatari (2006), who investigated eigenmodes in a linear shear flow on an equatorial β-plane.
Outside the continuous spectrum, that is, above the lines |λ| = m|α| in the diagrams, fast MR waves remain discrete eigenmodes even in the non-Malkus field.Their semi-analytical solutions can be obtained from  10) is reduced when B = µ and m 2 α 2 /λ 2 is small (see Appendix B).Moreover, we find that the lowermost wavenumber 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 4).This is because the Lorentz force does not act on this eigenmode.This mode is explained in detail in Appendix C.
Figure 6 depicts the eigenfrequencies for m = 2.They roughly epitomise the diagrams when m ≥ 3 (not shown).Their outlines do not change much from those of Figure 4 with the exception of the absence of the branch of the fast MR waves that penetrates the continuous spectrum.As in the m = 1 case, more conspicuous blue upward wedges exist around 10 −2 ≤ |α| ≤ 10 −1 and 10 −4 ≤ λ ≤ 10 −2 in the right panels of Figure 6.

Eigenfunctions of the Alfvén continuous modes
We investigate the eigenfunctions of the Alfvén continuous modes in this subsection and Appendix D. Figure 7 presents the typical structures of the perturbations in the stream function ψ 1 and the magnetic vector potential a 1 for the continuous modes.The dependences of their amplitudes ψ and ã on the colatitude are illustrated in Figure 7(a) and used for preparing the contour maps of Figure 7(b).From these figures, we note that spiky singular structures appear at the critical latitudes of the eigenmode.As shown in Section 2.2, the eigenfunction   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 |α| is sufficiently smaller than unity displays this behaviour (see Figure 11).Furthermore, the comparison of the eigenmodes having the same value of λ/m|α| in Figures 9 and 10 shows that the smaller the value of |α| is, the smaller the typical north-south wavelengths of their amplitudes.In Section 4, we will therefore examine the behaviour of the wave packets possessing large wavenumbers at small values of |α|, which may be applied to the Earth's core conditions, based on the ray theory to better understand our numerical results.
For moderate or large values of |α|, less striking differences exist in the evanescent property between the retrograde and prograde modes which possess the same absolute value |λ| of their angular frequency.This statement is based on Figure 11, which shows the heatmaps for m = 1 and |α| = 1, and other experiments with   several values of m and |α| (not shown).Therefore, the contrast in the property between the retrograde and prograde modes as demonstrated in Figures 9 and 10 can be attributed to the planetary β effect, that is, the effect of rotation.The ray-tracing analysis and 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 λ = −1/2 is discernible in the upper panels of Figure 11 (see also Appendix C).We attempted to use these plots to discover buried discrete eigenmodes other than this fast MR mode, although no such eigenmodes were found.
The evanescent property can be determined using the function which appears in an alternative form of the differential equation ( 10) The Mercator projection transformation y = (1/2) ln[(1 + µ)/(1 − µ)] yields a differential equation of the harmonic oscillation, which shows that the sign of L 2 determines the evanescent property of ψ√ Λ at the latitude.We can also rewrite (25b) into the form which, similarly, shows that the value of L 2 indicates the evanescent property.The equivalent of this differential equation was derived by Gilman and Fox (1999a), although the form of 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 (∂/∂t) + Ω 0 (∂/∂ϕ) in the non-rotating frame.The left panel in Figure 12 presents contour plots of L 2 as a function of λ and the colatitude when |α| = 0.01 and B = µ.We observe that the area where L 2 > 0 (or L 2 > −1) in the panel certainly agrees with the wavy regions in Figure 10.We now consider the case where |α| is small.On noting that λ = O(|α|) and |Λ| = O(|α| 2 ) for the continuous modes, we obtain , the oscillatory condition L 2 > 0 on the equatorial side (µ 2 < λ 2 /m 2 α 2 ) of the critical latitudes therefore requires that λ < 0 (the retrograde continuous modes), and L 2 > 0 on the polar side (µ 2 > λ 2 /m 2 α 2 ) for the prograde modes (λ > 0).This explains the contrasting evanescent behaviour between the retrograde and prograde continuous modes for a small value of |α|.

Interpretation in terms of ray theory and discussion
To better understand the continuous modes and their eigenfunctions, we reduce the system investigated thus far to a more restricted situation where |α| ≪ 1, which is appropriate for the Earth's core.We apply the ray theory, in which an inhomogeneous background field varies with a spatial scale much larger than typical wavelengths.
In addition, we track the path of a wave packet migrating with its group velocity.It should be noted that we now use θ as a dependent variable instead of µ to define a north-south wavenumber in terms of θ for the ray theory and consider B as the function of θ rather than µ.We introduce the local coordinates (Θ, Φ), which suitably measure the spatial scale size of the typical wavelength of a wave train.The introduction of small parameters aids in incorporating such a setting into the governing equations.In Section 3.1, we found that the typical meridional wavelengths decrease as the value of |α| decreases.Thus, it would be reasonable to select |α| as the parameter, if the value of |α| is sufficiently smaller than unity.The local coordinates are then stretched in the forms The temporal scale of the wave period (λ −1 = O(|α| −1/2 )) is similarly far from that of the migration of the wave train.The new shrunk time T , which is useful for measuring the latter, is given by The values of the exponents of |α| of the above variables are determined in Appendix E.
We then introduce a locally defined wavenumber and angular frequency, which depend on the global coordinates (θ, ϕ) 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 regarding wave dynamics (e.g.Lighthill, 1978), and we present an explanation of the corresponding details in Appendix E.Here we summarise the results.The expression of the perturbations of the stream function postulated in Section 2 is rewritten here as ψ 1 ≡ Re[M (ϕ, θ, T )e iφ L (Φ,Θ,τ ) ], where M is the wave amplitude, and φ L is the phase of the wave packet.With this ansatz, the local wavenumber and local nondimensional angular frequency are expressed as It should be noted that φ L depends on the local coordinates (Φ, Θ, τ ), while M , k, l, and λ depend on the global ones (ϕ, θ, T ).The local dispersion relation for our present problem is obtained from the leading order terms in the governing equations ( 8) in the form where λ s ≡ |α| −1/2 λ is the scaled nondimensional angular frequency.On replacing sin θ and B in (30) with constants, we obtain an equivalent to the nondimensional dispersion relation on a middle latitude β plane (Zaqarashvili et al., 2007).It should be noted that we denote λ s = λ H (ϕ, θ, k, l, T ) as the solution of (30) for λ s .From (30), the components of the nondimensional local group velocity c g are given by We eventually derive the ray-tracing equations and where the material time derivative moving with the local group velocity is with 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 θ.Furthermore, (32) shows that k sin θ and λ s (or λ) are invariant along a ray trajectory, but l is not.
We conduct the numerical time integration of the ray-tracing equations ( 32) with (31) for the movement of a wave packet originating at a given initial position (θ, ϕ) with a given initial local wavenumber (k, l) and a local dimensionless angular frequency λ determined by the local dispersion relation (30) (e.g.Teruya et al., 2022).In our code, the initial local longitudinal wavenumber k init is calculated from the relation (30) with the scipy.optimize.fsolvefunction of the SciPy library after the initial values of l, λ s , ϕ, and θ are specified.The succeeding time integration of ( 32) is based on an explicit Runge-Kutta method of order eight (the DOP853 algorithm in the scipy.integrate.solve_ivpfunction of the SciPy library).This integration is conducted without explicitly using (30) on the way, and the numerical errors in our calculations are monitored based on the value of the function D of (30).Their results are also compared to those of Section 3.1 in terms of the evanescent property.
Before demonstrating the 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 λ s = λ H (ϕ, θ, k, l, T ) at any time.Figures 13 and 14 present contour plots of the scaled dimensionless angular frequency λ s as a function of k and l for three different latitudes, when B = cos θ and B = 1, respectively.For ease of understanding Figure 13, we first explain Figure 14.These diagrams are obtained based on the equation transformed from (30) in the form which corresponds to (1) for global modes when B = 1.The wavenumber |k| sin θ corresponds to m, and k 2 + l 2 corresponds to n(n + 1) (see also Appendix F).In Figures 13 and 14, we take the plus of the plusminus sign in (33a) such that λ s is 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 14 correspond to the relations describing the propagation properties of wave packets that belong to prograde (λ s /k > 0) and retrograde (λ s /k < 0) Alfvén waves λ s ≃ ±|k| sin θ.Since this dispersion relation is independent of l, the contours become vertical.Furthermore, it can be concluded that, in Figure 14, the circular contour lines that are tangent to the line k = 0 represent the dispersion relation for fast MR waves λ s ≃ −k sin θ/(k 2 + l 2 ), and that the slightly curved part of the nearly vertical contour lines near the line l = 0 on the half plane λ s /k > 0 explains how slow MR waves λ s ≃ |k| sin θ(k 2 + l 2 ) propagate.The similarities between Figure 14 and the left and middle panels of Figure 13 suggest that the same is true for the case where B = cos θ.However, for B = cos θ, no branches of Alfvén and slow MR waves exist at the equator (θ = 90 • ), as shown in the right panel of Figure 13, since the main field vanishes there.It should be noted that the direction of the gradient (∂λ H /∂k, ∂λ H /∂l) at a point (k, l) in these plots is the same as that of the group velocity of a wave packet having a dominant local wavenumber (k, l) at the colatitude θ.For a wave packet belonging to either Alfvén or slow MR waves, the sign of the azimuthal component c g,ϕ of its group velocity is identical to that of the azimuthal component |α|(λ s /k) of its nondimensional local phase velocity, while those of the meridional components (c g,−θ and |α|(λ s /l)) are the opposite for the retrograde Alfvén packet.
Figure 13 illustrates the remarkable feature that the north-south component c g,−θ of the group velocity vanishes at l = 0 and l = ±∞, as can also be observed from (31b).The wave train can then be refracted at or absorbed into the latitude, moving only in the ϕ direction there (e.g.Acheson, 1972;McKenzie, 1973;Eltayeb, 1977;Eltayeb and Mckenzie, 1977;Grimshaw, 1979).In particular, from (30), the latter situation l 2 → ∞ with a reasonable condition λ s k sin θ ̸ = 0 results in the limit λ 2 s − k 2 B 2 sin 2 θ → 0, which indicates that the latitude is a critical one (Λ = 0).It follows that the nearly vertical lines for the Alfvén waves in Figure 13 should be linked to the rays corresponding to the Alfvén continuous modes observed in the results in Section 3. It should be noted that, although the nearly vertical lines in Figure 14 are similar to those in Figure 13, the Malkus field  The local dispersion relation (30) also allows us to understand the evanescent property of waves from the sign of the squared local meridional wavenumber l 2 .This value can be calculated from Figure 15 presents contour plots of l 2 as a function of k and λ s for three different latitudes in the case where B = cos θ.The 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 λ s = ±k|B| sin θ in the left and middle panels correspond to the relations for prograde Alfvén and slow MR waves (λ s /k > 0) and retrograde Alfvén waves (λ s /k < 0).Again, as shown in the right panel, no areas for Alfvén and slow MR waves are found at the equator (θ = 90 • ), since the background field vanishes there.
The propagation properties of the fast MR waves are indicated by the domain near the line k = 0 on the half plane λ s /k < 0. It should be noted that the wave packets that approach the curved lines l = 0 and the lines λ s = ±k|B| sin θ, where l 2 → ∞, are refracted at or absorbed into the corresponding latitude.Figure 15 aids in the short-term prediction of the migration of a wave packet.We now consider the two cases: (i) when the wave packet moves toward its corresponding critical latitude, and (ii) when the packet proceeds in the direction away from the critical latitude.Since k sin θ and λ s 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, λ s ) in the horizontal direction of the panels in Figure 15 (unless the outlines of their contour plots change significantly depending on the latitude).1.When a wave train moves toward the equator with λ s /k and l 2 positive (then λ s /l < 0 in the northern hemisphere from (31b) or Figure 13), k decreases and the point (k, λ s ) approaches the line λ s = k|B| sin θ from the right on the plots of Figure 15.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, which results in evanescent waves.It should be noted that, although plots for the Malkus field B = 1, which are illustrated in Figure 16, are similar to those in Figure 15, except for its right panel (θ = 90 • ), the point (k, λ s ) never reaches near the line λ s = ±k|B| sin θ owing to the constancy of k sin θ and λ s (unless the initial condition has already approximately satisfied this equality).When λ s /k is negative (under B = cos θ), a wave train that travels poleward (λ s /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, λ s ) approaches the line λ s = −k|B| sin θ from the left.It follows that the train belonging to retrograde Alfvén waves approaches the critical latitude from the equatorial side.
2. A wave packet moving in the opposite direction from its corresponding critical latitude is realised by a local meridional phase velocity |α|(λ s /l) which has a sign opposite to that in case (i).In other words, we focus on a packet that travels poleward (λ s /l > 0 in the northern hemisphere) with λ s /k being positive and one that moves equatorward (λ s /l > 0 in the northern hemisphere) with λ s /k being negative.The point (k, λ s ) then 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 from the polar side.
Whether a wave packet is to be refracted at or absorbed into a latitude where l 2 → ∞, or λ 2 s = k 2 B 2 sin 2 θ, is considered now.For a general profile of B, near such a colatitude θ = θ c , (33b) becomes if (dB 2 /dθ)| θ=θc ̸ = 0.In particular, when B = cos θ, the oscillatory condition l 2 > 0 requires that (θ−θ c ) cos θ c has the sign opposite to that of λ s /k; the packet that belongs to either prograde Alfvén or slow MR waves (λ s /k > 0) approaches the critical colatitude θ c from the polar side ((θ − θ c ) cos θ c < 0), while the retrograde Alfvén one (λ s /k < 0) does from the equatorial side ((θ − θ c ) cos θ c > 0).This is just a mathematical paraphrasing of the consideration of the manner in which the packets approach the critical latitudes mentioned in the previous paragraph.Since the term λ s k sin θ in ( 30) -which has the same sign as the numerator of ( 34)represents the planetary β effect, the aforementioned distinction between the prograde and retrograde waves is caused by the β effect.Moreover, on using ( 31) and ( 34), we can obtain an asymptotic expression of the group velocity This expression provides the travel time of the packet from a given latitude θ in the vicinity of the critical latitude to the latter in the form 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 explained in a similar fashion in Appendix F. The latitudes are often referred to as "turning latitudes".
We finally demonstrate the ray trajectories obtained from the numerical time integration of the ray-tracing equations, although rough trajectories can be sketched from the above examinations.Figures 17 and 18 present 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, and its position evolves in accordance with (32a).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 ϕ of their positions at time T are recorded using the same colour scheme as their upper panels.The numerical errors for the results shown in Figures 17 and 18 are diagnosed through D in (30) as D ≲ 10 −3 throughout their numerical integrations.The trajectories agree with the above predictions in terms of the refraction at the turning latitudes, the absorption into the critical ones, and the incident directions to those.From the asymptotic expression (35a) of the group velocity near a critical latitude, we obtain the period for a wave packet to circle along the latitude line as sgn(Ω 0 )T = |2πk sin θ/λ s |, which approximates the periods read from the lower right panels of Figures 17 and 18.
The facts that the wave packets cannot cross their corresponding critical latitudes and that whether the packets approach there from the polar or equatorial sides depends on the sign of λ s /k are consistent with the features observed from the numerical results for the continuous modes (Figures 9 and 10) in the eigenvalue    problem for global modes when |α| is small.Although the spatial scale of the 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 where |α| ≪ 1, the ray trajectories for Alfvén waves facilitate the rough prediction of the behaviours of the continuous modes without actually solving the eigenvalue problem.This approach is applicable even if the imposed field possesses more than two singular latitudes unlike the case where B = cos θ.For instance, four of the trajectories for another equatorially antisymmetric field B = sin θ cos θ are displayed in Figures 19 and 20.For this profile, two critical latitudes exist in each hemisphere when 0 = min(B 2 ) < λ 2 s /k 2 sin 2 θ < max(B 2 ) = 1/4.It should be noted that this field profile does not satisfy the condition for a smooth field at the poles, as mentioned in Section 2. However, we have selected this profile for simplicity only for now.The starting points (the black asterisks) of the rays in these figures are selected such that l 2 > 0, and the rays must stay in regions where l 2 > 0. As discussed in the previous paragraphs, the regions where l 2 > 0 are bounded by critical latitudes (and turning latitudes).For prograde waves (Figure 19; λ s /k > 0), the starting points should lie between the two critical latitudes in each hemisphere.For retrograde waves (Figure 20; λ s /k < 0), the starting points should lie either on the polar side or the equatorial side in each hemisphere.This corresponds to the prograde and retrograde continuous modes being evanescent in the polar and equatorial regions and the mid-latitudes, respectively.The value of the function L 2 of (25a) also provides information regarding 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 E).However, its approximate formula (27) for a small |α| agrees with (34). Figure 12(b) depicts the values of L 2 when B = µ 1 − µ 2 (= sin θ cos θ), 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.
Moreover, we learn why discrete branches of slow MR waves disappear under non-Malkus fields, as shown in the numerical results in Section 3, from the ray-tracing approach.Figure 21 presents a trajectory for a wave packet that (at least initially) belongs to the slow MR wave, because the condition has been satisfied by its initial condition.This inequality is obtained from the comparison between the first two terms in the local dispersion relation (30) 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 17).Since the propagation properties for slow MR waves and those of prograde Alfvén waves are continuous, as shown in Figures 13  and 15, the 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 (37) as l 2 → ∞, thus the Alfvén balance λ 2 s ≈ k 2 B 2 sin 2 θ is reached.From (33b), the latitudinal variation of l 2 during the migration of a wave packet is written as and we find that (∂l 2 /∂θ) k sin θ,λs is always positive (negative) in the northern (southern) hemisphere when B = cos θ and λ s /k > 0. Therefore, the mode conversion from slow MR into prograde Alfvén waves should have occurred between the initial colatitude θ init and the critical colatitude θ c in the example of Figure 21, since 0 ≤ l 2 ≤ l 2 init within the interval between θ 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 type of critical latitude: the Alfvén resonance λ 2 s = k 2 B 2 sin 2 θ.On the other hand, wave trains belonging to fast MR waves, which have discrete branches when λ < −m|α| 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).
The invariants including the products of perturbations are useful for understanding waves and their associated phenomena.As derived in Appendix E, our approximation leads to a conservation law in the form It follows from the constancy of k sin θ and λ s that the equation in which (∂D/∂λ s )|M | 2 is replaced by (k sin θ/λ 2 s )(∂D/∂λ s )|M | 2 is also correct.In our subsequent paper, we intend to 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), owing to its long travel time (36).The fact that the packet carries the pseudomomentum in line with (39) implies that the mean flow would be accelerated there because the damping of waves can cause angular momentum exchange between the waves and a mean flow in accordance with the wave-mean flow interaction theory (e.g.Bühler, 2009).This could possibly 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 confirm the validity of the ray theory when l 2 → ∞.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 We now examine whether the condition ( 40) is valid near its corresponding critical latitudes, where l diverges and |M | may also diverge.To obtain an asymptotic expression of |M | near the critical latitude θ c , we take advantage of the conservation law (39) in the form where S g is a region that moves and deforms with the local group velocity.Let θ 1 (T ) and θ 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 c g of the packet depends only on the latitudinal position θ, we obtain Thus, the latitudinal length |θ 2 − θ 1 | of the packet satisfies Using these relations, we can estimate that (∂D/∂λ s )|M |2 (∆ϕ sin θ) ∝ c −1 g,−θ with ∆ϕ being the longitudinal length of the packet and obtain M = O(|θ − θ c | −1/4 ).It should be noted that an identical asymptotic expression is also obtained from the steady problem in the last two paragraphs of Appendix E. Since l = O(|θ − θ c | −1/2 ) (from ( 34)) and M = O(|θ − θ c | −1/4 ), we cannot appropriately discuss the behaviour of a wave packet if it approaches within the distance |θ − θ c | = O(|α|) from its corresponding critical latitudes.The case where l = 0 is addressed in Appendix F.

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ϕ = B 0 B(cos θ) sin θ.In Section 3, we solved the eigenvalue problem for the simplest equatorially antisymmetric non-Malkus field B = cos θ.The most important finding is that the Alfvén and slow MR discrete branches disappear and a continuous spectrum appears instead.The eigenfunctions of the continuous modes and their concomitant critical latitudes were investigated in unprecedented detail and compared with the Frobenius series solutions in Appendix D. In Section 4, we conducted a ray-theoretical study in an inhomogeneous magnetic field to clarify the nature of the continuous modes from a perspective different from the mode theory.Generally, the ray theory is the theory of wave packets with slowly varying frequencies, wavenumbers, and amplitudes.In our problem, this theory has the advantage that it can be used for small absolute values of the Lehnert number α, where the mode theory encounters numerical difficulty.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.
The advent of the continuous mode can be interpreted in terms of the ray theory as wave packets that move toward a critical latitude and that are ultimately absorbed there.In particular, the wave packet belonging to the slow MR wave turns into the packet that belongs to the Alfvén wave before it is absorbed at its corresponding critical latitude.This indicates that the slow MR wave cannot be distinguished from the Alfvén wave in the case of a non-Malkus field and corroborates the idea that slow MR waves transform into Alfvén continuous modes.
We used the function L 2 to show that the observed difference in the evanescent property between the prograde and retrograde continuous modes results from the planetary β effect.This function corresponds with the square of the local meridional wavenumber l 2 in the ray theory, which provides the same explanation for the evanescent property.Whether the packet approaches the critical latitude from the polar or equatorial sides depends on the sign of the azimuthal component |α|(λ s /k) of its nondimensional local phase velocity, and this packet-approaching side is the oscillatory one (l 2 > 0) of the corresponding mode.
We note some interesting topics below that should be investigated in the future studies.
This work was supported in part by JSPS KAKENHI Grant Number JP24K07177 and JP24K00694, and performed with the support and under the auspices of the NIFS Collaboration Research program (NIFS24KIIC001).
We can derive another bound (19b) by using (44b) when Im(λ) ̸ = 0.The property of the Rayleigh quotient of the associated Legendre operator in the form unless Im(λ) = 0.This inequality is equal to (19b) since Re(λ) ≤ 0.

B Fast MR waves when B = µ
This appendix and Appendix C cover the fast MR waves when the simplest equatorially antisymmetric non-Malkus field B = µ is imposed.In particular, we now narrow our interests down to discrete branches λ < −m|α| outside the continuous spectrum.If we focus on waves with frequencies |λ| ≫ m|α|, sufficiently higher than the continuous modes, (25a) is reduced to where c 2 ≡ (m 2 α 2 /λ 2 )(7 + m/λ), and then (25b) approximately becomes the same form as the differential equation for the angular prolate spheroidal wave function S mn (c, µ) (e.g.Abramowitz and Stegun 1964; see also Zaqarashvili et al. 2009).The power series expansion for the eigenvalues λ mn ≡ −m/λ + m 2 α 2 /λ 2 of this differential equation yields the approximate dispersion relation for the fast MR waves outside the continuous spectrum in the form Their approximate angular frequencies λ approx were calculated from the above relation with the scipy.optimize.fsolvefunction of the SciPy library.Their results for m = 1 are shown in Figure 22, overlaid on the same graphs as the left panels of Figure 4.When |λ| ≫ m|α|, their values surely approximate the angular frequencies obtained numerically from the full eigenvalue problem. Figure 23 displays two of their approximate eigenfunctions S mn / λ 2 approx − m 2 α 2 µ 2 (≃ ψ), where S mn were given by the scipy.special.proang1 function of the SciPy library, with their corresponding eigenfunctions calculated from the full problem.It should be noted that, if c 2 = (m 2 α 2 /λ 2 approx )(7 + m/λ approx ) < 0, the parameter c in S mn (c, µ) is replaced by −ic, and then S mn (−ic, µ) is regarded as the angular oblate spheroidal wave function.In this case, we used the scipy.special.oblang1 function of the SciPy library.
We also integrated the ray-tracing equations (see Section 4) numerically for wave packets that belong to fast MR waves.Figure 24 shows three of their trajectories when B = cos θ.They move back and forth between two If both sides of both (50a) and (50b) did not vanish, the above equations would be transformed into the same form as (10) and inevitably possess critical latitudes µ 2 c = 1/4α 2 .We, therefore, presume that a solution for which both sides of (50a) are identically equal to zero exists.Because ∇ 2 h P 1 1 = −2P 1 1 , ∇ 2 h P 1 2 = −6P 1 2 , and µP 1 1 = P 1 2 /3 are derived from ( 21), the presumption requires that its eigenfunction is written as 25).This implies that this discrete mode buried in the continuous spectrum is exceptional.On its branch, the equipartition between MKE and MME occurs only when |α| = 5/12 ≈ 0.6455, and we found that the value agrees with that which one can read from the colours of markers in the upper left panel of Figure 4. Furthermore, we notice that its angular frequency λ = −1/2 satisfies even the approximate dispersion relation (49) for fast MR waves "outside" the continuous spectrum, since one then has c 2 = 20α 2 and λ 11 = 2 + 4α 2 .Nevertheless, its (exact) eigenfunction ψ[1] N 1 1 P 1 1 is not identical with the approximate eigenfunctions S 11 / √ Λ for the fast modes outside the continuous spectrum (see Appendix B).
It should be added, at the end of this appendix, that Wang et al. (2022b) recently studied two distinctive modes in the same system as our problem.One is the above buried discrete mode (λ = −1/2), and the other is the stationary mode (λ = 0) whose eigenfunction is represented as ψ 1 = 0 and a 1 ∝ P 1 2 e iϕ in (50).Their paper showed that the latter is closely linked to "clamshell instability" (see Section 1) and referred to this mode as the "tilting mode".On the other hand, we consider that the former is also interesting, because this mode remains stationary in the non-rotating frame although this is a rotating fluid problem, as shown in Wang et al. (2022b).Those two modes are similar in terms of the tilt of a toroidal magnetic field.For the buried mode of the fast MR wave, the perturbation ψ 1 in the stream function, which is proportional to P 1 1 e iϕ , makes the imposed field B 0ϕ = B 0 sin θ cos θ tilt, hence the perturbation a 1 ∝ P 1 2 e iϕ in the magnetic vector potential.Such a total magnetic field does not obviously induce the Lorentz force acting on the perturbation.It should be noted that the right-hand side of (50a) vanishes when sgn(α)ã/ √ ρ 0 µ m ∝ P 1 2 .

D 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 (17).Let ψnum be one of the numerical stream functions, and we fit it into the form C I ψ(c) I + C II ψ(c) II by adjusting the coefficients C I and C II on each side of the critical latitudes.The fitting procedure is as follows.The colatitude (0 ≤ θ ≤ π) is divided into N θ points at even intervals.On each side of the nearest point of a singular latitude, the N data points closest to the point are chosen among the N θ points for the fitting.For these points θ i (1 ≤ i ≤ N data ) either on the equatorial or polar sides, we may write Now, ψ(c) I and ψ(c) II are approximated by the second-order Frobenius solutions with (18).We obtain a candidate value for each of C I and C II from each equation, through least squares fittings of (51) with the numpy.polyfitfunction of the NumPy library.Thereby we have four candidate values for each of C I and C II for one critical latitude, two from the equatorial side, and two from the polar side.The upper panels of Figure 26  The above procedure is also applied to all the continuous modes obtained from our numerical calculations.Figure 27 This outcome demonstrates that the values of I num appear to be compatible with our expectation stated in Section 2.1; the values are arbitrary numbers and adjusted to satisfy boundary conditions.In addition, we observe that the value of sgn(µ c )[C ] has extremums at the zeros of C II .If discrete eigenmodes without logarithmic and step function singularities are buried in the continuum, the two values must simultaneously vanish.In fact, we attempted to use these values as a means to reveal buried discrete eigenmodes, though no such eigenmodes have been found.

E Derivations of the ray-tracing equations and their related formulae
This appendix explains detailed derivations of the ray-tracing equations and their related formulae, which were skipped in Section 4. We first introduce the parameter given by s ≡ − sgn (ln |α|).Then, |α| s ≤ 1 is always satisfied regardless of the magnitude of |α|.If we only consider the situation where either |α| ≫ 1 or |α| ≪ 1, the parameter |α| s becomes sufficiently smaller than unity.Utilising this small parameter, we can express the separations of the spatial and the temporal scales, as explained in Section 4. In other words, with the local coordinates (Φ, Θ, τ ) and the global ones (ϕ, θ, T ), we can write where p, q > 0.
In such separations, the local coordinates are supposed to be suitable to measure the phase φ L of waves.The local wavenumber and local nondimensional angular frequency, which may depend only on the global coordinates, are accordingly defined as  where φ G (ϕ, θ, T ) (= φ L (Φ, Θ, τ )) denotes the phase function whose independent variables are the global coordinates.Now we recognise that k, l = O(|α| 0 ) by construction, and that the zonal wavenumber m = O(|α| −sp ) measured in the global coordinates is large.This implies that λ = O(|α| s(q−p) ), and we introduce the scaled nondimensional angular frequency λ s ≡ |α| s(p−q) λ.In addition, (54) yield the curl-free conditions If we consider (ϕ, −θ) to be generalised coordinates in analytical mechanics, (k sin θ, l) and |α| sp φ G , from (54a) and (54b), correspond to generalised momentums and a generating function, respectively.The first equation of ( 55) is also equivalent to the property {k sin θ, l} = 0 of the Poisson bracket.On the basis of these facts, we may raise the topic of the Hamilton-Jacobi equation where λ H is the "Hamiltonian" for our problem.We however will not pursue the solution |α| sp φ G of this equation.Upon comparing (54c) to ( 56) and treating λ s as an independent variable, we obtain a dispersion relation in the form The nondimensional local group velocity with respect to the dimensionless time τ and the coordinates (ϕ, −θ) is then written as where the material time derivative moving with the local group velocity is expressed as (d g /dT ) ≡ (∂/∂T ) + |α| −sq (c g • ∇ G ).Using the chain rule, one can also get (59b) and (59c) from ( 54), ( 55), (56), and (58).
The difference between our ray-tracing approach and standard analytical mechanics is that we do not initially know an explicit expression of the Hamiltonian.To obtain it, we need to substitute the ansatz ψ 1 ≡ Re[M (ϕ, θ, T )e iφ L (Φ,Θ,τ ) ] into the original partial differential equation the first term on the left-hand side to the third term on the right-hand side.A plausible balance among them requires that 2(q − 2p) = q − 2p = 2s − 4p ( < 2s − 3p < 2s − 2p ).
We then have s = 1 (|α| ≪ 1), p = 1/2, and q = 1.Because possible balances for s = −1 do not exist, our approximation is probably inappropriate for the case where |α| ≫ 1.It follows that the leading order terms O(|α| 0 ) in ( 60 It should be noted that this equation is complex at a glance, but conforms to the same rule as the structure that the counterpart in Bretherton (1966) has.At last, the novel conservation law (39) is obtained from ( 31), (55), and ( 63) with (57).
In the remainder of this appendix and Appendix F, we shall set about the steady problem in which wave trains with specified values of k sin θ and λ s uniformly continue to be injected on a line of latitude.The discussion of the problem will provide the behaviours of the steady trains near critical latitudes, where l 2 → ∞.Its findings are also helpful in learning the behaviours of the trains near turning latitudes, at which l = 0 (see Appendix F).First of all, consider a region sandwiched between two ray trajectories drawn by the trains injected at two distinct longitudes (ϕ 1 and ϕ 2 , say).This region may be called a "ray tube".It should be noted that, because of the symmetry about the axis of rotation, one of the two trajectories can overlap the other by the rotation |ϕ 1 − ϕ 2 | about the axis.The distance between the two trajectories along a line of a colatitude θ is, therefore, |ϕ 1 − ϕ 2 | sin θ.We then prepare two cross-sections of the tube along two lines of latitude.Let S be the region bounded by the two sections on the tube.For the steady problem with the uniform injection, the surface integral of ( 39 .This means that steady trains jam into such a latitude and that their amplitudes increase there.It should be noted that the same is true for an isolated packet, as shown in the last paragraph of Section 4. For constant values of k sin θ and λ s , we can get the same expression for |M | as the above in a different way.The solution is again redefined as ψ 1 ≡ Re[ ψ(θ)e iφ ] with φ = |α| −1/2 (k sin θ)ϕ − |α| 1/2 λ s τ , where we conjecture that ψ can be expressed by the WKBJ (Wentzel-Kramers-Brillouin-Jeffreys) form

F Behaviours of wave packets near a turning latitude
The last appendix is devoted to an investigation into the behaviours of wave packets near turning latitudes, where l = 0, when |α| ≪ 1 (see also Section 4 and Appendix E).We here write l 2 = (l 2 t ) ′ (θ − θ t ) + O(|θ − θ t | 2 ) in the vicinity of a turning colatitude θ = θ t , where (l 2 t ) ′ ≡ (dl 2 /dθ)| θ=θt (and we assume that (l 2 t ) ′ ̸ = 0).If a packet advances toward its corresponding turning latitude from a nearby latitude θ, it can then reach there in a finite time It should be noted that, in the above equation, either 1 < (−mλ)κ 2 < κ 2 or 0 < (−mλ)κ 2 < 1 < κ 2 holds, because we find that 1 + mλ > 0 during the derivation of (78).To evaluate the integral in (78), we consider the integral of a complex-valued function of z (≡ x + iy) in the form where a ≡ min(1, √ −mλκ), b ≡ max(1, √ −mλκ), γ ≡ 1 − (a/b) 2 , and χ ≡ γκ 2 /(κ 2 − a 2 ).In addition, the contours C 0 , C 1 and C 2 of integration are shown in Figure 28, and the substitution y = ab/ b 2 − (b 2 − a 2 )v 2 has been used in the last equal sign of (79).When the radius of the semicircle bounded by the contour C 0 goes to infinity, the left-hand side of (79) becomes equal to the integral in (78).Because the integrals in the rightmost side of (79) are the complete elliptic integrals of the first and third kinds for the elliptic modulus γ and the characteristic χ (e.g.Abramowitz and Stegun, 1964), its calculation requires a numerical method.Figure 29 shows the approximate angular frequencies obtained from (78), the calculation of which was performed by the scipy.optimize.fsolveand the scipy.integrate.quadfunctions of the SciPy library.The condition (75) will provide the approximate dispersion relation for discrete modes when any basic field is imposed and can be utilised for the check of the numerical calculation of the full eigenvalue problem for such a field.

Figure 1 :
Figure 1: Schematic illustrations of the profiles of the imposed toroidal magnetic fields B 0ϕ ∝ B sin θ.

Figure 3 :
Figure 3: Ratio of the mean kinetic energy MKE to the mean total energy MKE + MME against the absolute value |α| of the Lehnert number when the zonal wavenumber m = 1 and the Malkus field is imposed.This plot is obtained from the relation MKE/(MKE + MME) = λ 2 /(λ 2 + m 2 α 2 ), where the nondimensional angular frequency λ is calculated using (1).As |α| → 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 2.

Figure 4 :
Figure 4: Eigenfrequencies versus the Lehnert number when the zonal wavenumber m = 1 and the simplest equatorially antisymmetric non-Malkus field B = µ pervades the system.The ordinates of each panel represent the real part Re(λ) of the dimensionless angular frequency, and the abscissas present the absolute value |α| of the Lehnert number.Retrograde modes (Re(λ) < 0) and prograde modes (Re(λ) > 0) are presented in the left and right columns, respectively.The upper and lower rows illustrate the 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.

Figure 5 :
Figure 5: Dependence of the real parts Re(λ) of the eigenvalues on the truncation degree N t when the zonal wavenumber m = 1, the absolute value of the Lehnert number |α| = 1, and the background field is the simplest equatorially antisymmetric non-Malkus one (B = µ).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 where N t = 2000 and 1999, respectively.

Figure 6 :
Figure 6: Same as Figure 4, but for m = 2.The two white asterisks of the upper right panel correspond to the two eigenmodes depicted in Figure 8.

Figure 7 :
Figure 7: Eigenfunction of the sinuous mode with the nondimensional angular frequency λ ≈ 0.05006 when the zonal wavenumber m = 1, the absolute value of the Lehnert number |α| = 0.1, and the simplest equatorially antisymmetric non-Malkus field B = µ is imposed.The critical colatitudes θ c ≈ 59.96 • and 120.04 • .(a) Amplitudes of the stream function ψ (red line) and the scaled magnetic vector potential sgn(α)ã/ √ ρ 0 µ m (blue line) as functions of the colatitude.(b) Contour maps of the stream function ψ 1 (left panel) and the scaled magnetic vector potential sgn(α)a 1 / √ ρ 0 µ m (right panel) in the Mollweide projection.

Figure 8 :
Figure 8: Eigenfunctions of two magnetic-energy dominant modes for m = 2 and |α| = 0.013.The axes and colours are identical to those in Figure 7(a).

Figure 9 :
Figure 9: Amplitudes of the eigenfunctions of all the obtained continuous modes when the zonal wavenumber m = 1, the absolute value of the Lehnert number |α| = 0.1, and the basic field is the simplest equatorially antisymmetric non-Malkus one (B = µ).The four panels are divided into the stream function and magnetic vector potential in the left and right columns, respectively, and sinuous and varicose modes in 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 λ.The darker the colours, the higher the absolute values | ψ| and |ã| of the amplitudes of the stream function and the magnetic vector potential.

Figure 12 :
Figure 12: Dependence of the function L 2 given by (25a) on the nondimensional angular frequency λ and the colatitude when the zonal wavenumber m = 1, the absolute value of the Lehnert number |α| = 0.01, and the basic fields are (a) B = µ and (b) B = µ 1 − µ 2 .The solid and dashed curves present the contour lines corresponding to L 2 = 0 and −1, respectively.

Figure 15 :
Figure 15: Local dispersion relation written as (33b) when the simplest equatorially antisymmetric non-Malkus field B = cos θ 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 λ s = |α| −1/2 λ for the colatitudes θ = 30 • (left panel), 60 • (middle), and 90 • (right).

Figure 17 :
Figure 17: Ray trajectories in the Mollweide projection for wave packets belonging to prograde Alfvén waves (the scaled zonal wavenumber k sin θ ≈ 1.22338, and the critical colatitude θ c ≈ 35.17 • in the northern hemisphere) when the simplest equatorially antisymmetric non-Malkus field B = cos θ pervades the system (upper panels).The scaled dimensionless angular frequency λ s = |α| −1/2 λ = 1.The black asterisks correspond to their initial position (θ init , ϕ init ) = (30 • , 0 • ).The directions (k, l) of their local phase velocities are indicated by hue, the colour scales of which are presented in the left lower panels.The lower right panels present the longitudes ϕ of their positions against the scaled nondimensional time T = |α|τ with the colour denoting the directions of their local phase velocities.It should be noted that both ϕ and T are multiplied by sgn(Ω 0 ) such that the "time" sgn(Ω 0 )T is always positive even when the rotation rate Ω 0 of the sphere is negative.

Figure 22 :
Figure 22: Same as the left panels of Figure 4, but for a different plotting range (grey markers).The red circles represent the approximate dimensionless angular frequencies λ approx calculated by the approximate dispersion relation (49).

Figure 23 :
Figure 23: Same as Figure 7(a), but for fast MR waves outside the continuous spectrum.The black dashed curves illustrate their approximate eigenfunctions S mn / √ Λ.The eigenfunctions have been scaled such that their amplitudes become comparable to those of their corresponding eigenfunctions which are obtained numerically from the full eigenvalue problem.

Figure 25 :
Figure 25: Same as Figure 7(a), but for the fast MR wave buried in the Alfvén continuous spectrum (the zonal wavenumber m = 1 and the nondimensional angular frequency λ = −0.5).
show a result when we perform this procedure with N θ = 7201, and N 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 I typically has different values (C ) between the two sides of a critical latitude, whilst C II has the same value on both sides.Since C III in (14a) can also be written as sgn(µ c )[C our numerical eigenmodes are consistent with the general results as described in Sections 2.1 and 2.2.We accordingly adopt the mean values of the two candidate values for each of C their definite values, which are used in the graph comparing ψnum with C I ψ(c) I + C II ψ(c) II (the lower panel of Figure26).
depicts their values of sgn(µ c )[C (p) I − C (e) I ] (= C III , red circles) and C II (blue circles) for m = 1 and |α| = 0.1 in the left ordinates as functions of λ.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 (14d), which is written in the present instance as at the extremum points of C II , and that sgn(µ c )[C

Figure 26 :
Figure 26: Comparison between the numerical eigenfunction ψnum of the sinuous mode with the dimensionless angular frequency λ ≈ 0.05006 (the critical latitude θ c = 59.96 • in the northern hemisphere) and linear combinations C I ψ(c) I + C II ψ(c) II of the Frobenius series solutions (17) with (18) when the zonal wavenumber m = 1, the absolute value of the Lehnert number |α| = 0.1, and the simplest equatorially antisymmetric non-Malkus field B = µ permeates the system.The undetermined coefficients C I and C II are estimated from least squares fittings of ψnum / ψ(c) I = C II ( ψ(c) II / ψ(c) I ) + C I (upper left panel) and ψnum / ψ(c) II = C I ( ψ(c) I / ψ(c) II ) + C II (upper right panel) with N 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 θ = 7201.The lower panel shows ψnum (black curve) and C I ψ(c) I + C II ψ(c) II with C I and C II determined by the fittings (red and blue curves) as functions of the colatitude.The vertical grey shaded area of this panel contains the 2N data = 400 points used in the fittings.

Figure 27 :
Figure 27: Values of sgn(µ c )[C (p) I − C (e) I ] (red circles and the left vertical axes), C II (blue circles and the left vertical axes), and I num defined as (52) (black circles and the right vertical axes), relevant to the coefficients of the Frobenius series solutions against the nondimensional angular frequency λ.The case where the zonal wavenumber m = 1, the absolute value of the Lehnert number |α| = 0.1, and the simplest equatorially antisymmetric non-Malkus field B = µ is imposed is shown in the left and the right panels for sinuous and varicose modes, respectively.N θ and N data are set to 7201 and 200, respectively, in the same way as Figure 26.

Figure 28 :
Figure28: Path of integration for the integral (79).The contours C 0 , C 1 and C 2 are represented as thick black curves.The black and red circles are poles and branch points, respectively, and the red lines correspond to branch cuts.The integral in (78) must be evaluated along the path represented as the cyan line.

Figure 29 :
Figure 29: Same as Figure 22, but for the approximate dimensionless angular frequencies λ approx calculated by the approximate dispersion relation (78) instead of (49).