Solitary magnetostrophic Rossby waves in spherical shells

Finite-amplitude hydromagnetic Rossby waves in the magnetostrophic regime are studied. We consider the slow mode, which travels in the opposite direction to the hydrodynamic or fast mode, in the presence of a toroidal magnetic field and zonal flow by means of quasi-geostrophic models for thick spherical shells. The weakly-nonlinear, long waves are derived asymptotically using a reductive perturbation method. The problem at the first order is found to obey a second-order ODE, leading to a hypergeometric equation for a Malkus field and a confluent Heun equation for an electrical-wire field, and is nonsingular when the wave speed approaches the mean flow. Investigating its neutral, nonsingular eigensolutions for different basic states, we find the evolution is described by the Korteweg-de Vries equation. This implies that the nonlinear slow wave forms solitons and solitary waves. These may take the form of a coherent eddy, such as a single anticyclone. We speculate on the relation of the anti-cyclone to the asymmetric gyre seen in Earth's fluid core, and in state-of-the-art dynamo DNS.


Introduction
Linear waves in an inviscid, perfectly-conducting fluid permeated by a uniform magnetic field B 0 in a frame rotating with rate Ω satisfy the dispersion relation (Lehnert 1954) where ω is the frequency, k is the wavenumber vector, ρ is the density, and µ 0 the magnetic permeability. This yields a wide variety of magnetic Coriolis (MC) waves, including fast (modified inertial) and slow (magnetostrophic) waves; the latter being unique to rotating magnetohydrodynamics (MHD). In this manuscript we consider magnetostrophic waves for which (Ω · k) 2 /|k| 2 ≫ (B 0 · k) 2 /(ρµ 0 ): In particular, one class which has the relation ω ≈ − (B 0 · k) 2 |k| 2 ρµ 0 βk . (1.2) Here β denotes the beta parameter, k is the azimuthal wavenumber, and the minus sign indicates waves travel opposite to the hydrodynamic Rossby wave, ω = βk/|k| 2 . This class is sometimes referred to as slow hydromagnetic-planetary or magnetic-Rossby (MR) waves (Hide 1966). Relation (1.2) indicates they are dispersive, and depend on the background field and the wavelength; these waves have been suggested to be important in Earth's fluid core and for the geomagnetic westward drift (e.g. Hide 1966;Malkus 1967;Canet et al. 2014;Hori et al. 2015;Nilsson et al. 2020).
Other classes of MC waves include torsional Alfvén waves for which Ω · k ≈ 0 and (Ω · k) 2 /|k| 2 ≪ (B 0 ·k) 2 /(ρµ 0 ) (Braginskiy 1970;Roberts & Aurnou 2012;Gillet et al. 2015). More recently inertial-Alfvén waves (Bardsley & Davidson 2016) have been claimed to account for the geomagnetic jerks (Aubert & Finlay 2019). Laboratory experiments have identified several types of magnetostrophic waves in spherical Couette flows with a dipolar magnetic field being applied (Schmitt et al. 2008). We note the wave dynamics relies on both the direction and the morphology of the background magnetic field, as illustrated in the simple planar model (1.2). Here we focus on the problem with a purely azimuthal basic field; for this case (1.2) reduces to ω ∝ k|k| 2 , indicating its linear and cubic relationship to the azimuthal wavenumber.
The linear theory for MC waves in stably-stratified, thin layers is well-studied (e.g. Braginskiy 1967;Gilman 2000;Zaqarashvili et al. 2007;Márquez-Artavia et al. 2017) as observational exploration of the geomagnetic field and the solar corona has developed to reveal periodic patterns (Chulliat et al. 2015;McIntosh et al. 2017). Stratification in general introduces a correction term to the dispersion relations of MC waves, whilst in a thin layer the direction of travel is usually reversed; however, this is not always true in spherical geometries. The unstratified thick shell problem considered here is sufficient to provide some fundamental understanding of the nonlinear problem.
Theoretical investigation is expanding to consider their nonlinear properties such as turbulence (Tobias et al. 2007) and triadic resonances (Raphaldini & Raupp 2015). London (2017) found a couple of cases in which nonlinear equatorial waves in the shallow water MHD should be governed by Korteweg-de Vries (KdV) equations and so behave like solitary waves. They were mostly fast MR modes, recovering the equatorial Rossby wave soliton (Boyd 1980) in the nonmagnetic limit, but he reported one case in which the wave would slowly travel in the opposite azimuthal direction. Hori (2019) investigated magnetostrophic MR waves in a Cartesian quasi-geostrophic (QG) model. The slow, weakly-nonlinear waves led to evolution obeying the KdV equation unless the basic state -all the magnetic field, topography, and zonal flow -is uniform. Slow MR waves have been seen in spherical dynamo DNS travelling with crests/troughs that were isolated and sharp, unlike the continuous wave trains that might be expected (Hori et al. 2015(Hori et al. , 2018. Hydrodynamic Rossby wave solitons have been extensively studied, motivated by atmosphere and ocean dynamics (e.g. Clarke 1971;Redekopp 1977;Boyd 1980). In the long wave limit it has been demonstrated that the QG soliton relies on the presence of a shear in the basic flow or topography. Redekopp (1977) further analysed nonlinear critical layers arising from singularities as the wave speed approaches the basic flow speed, and discussed their relevance for the persistence of Jupiter's Great Red Spot.
The present manuscript demonstrates that weakly nonlinear slow MR waves in spherical containers yield soliton solutions. We adopt simple QG MHD models and asymptotically derive the evolution equation for the long wave when the basic magnetic field and flow are both azimuthal. We demonstrate: (i) the amplitude at the first order is described by the KdV equation for the chosen basic states, (ii) the problem is dictated by an ODE, which has no singularities as the wave speed approaches the basic flow speed, and (iii) the single soliton (solitary wave) solution to the KdV equation implies an isolated eddy that progresses in a stable permanent form on magnetostrophic timescales.

Theoretical foundations
We consider an inviscid, incompressible, ideal quasi-geostrophic (QG) model of electrically conducting fluid within a rapidly rotating shell, bounded by inner and outer spheres of radii r i and r o , respectively (e.g. Busse 1970;Gillet & Jones 2006). We use polar coordinates (s, ϕ, z) with rotation Ωẑ.
For rapid rotation, the incompressible horizontal QG fluid motion can be expressed as u ≈ ∇ × ψ(s, ϕ)ẑ with ψ a streamfunction, so it is independent of z. When the magnetic field is not too strong to violate the QG approximation, we further assume the magnetic field may be written as B ≈ ∇ × g(s, ϕ)ẑ with g being the potential (e.g. Busse 1976;Abdulrahman et al. 2000;Tobias et al. 2007;Canet et al. 2014). No penetration on the spherical boundaries at z = ±H = ± r 2 o − s 2 enables us to represent the Coriolis term of the axial vorticity equation in terms of the topography-induced beta parameter. The equations for the z-components of the vorticity and the magnetic potential in dimensionless form are then: where ∆ H = (1/s)∂/∂s(s∂/∂s) + (1/s 2 )∂ 2 /∂ϕ 2 , and J [f 1 , f 2 ] = (∂f 1 /∂s ∂f 2 /∂ϕ − ∂f 2 /∂s ∂f 1 /∂ϕ)/s for any functions f 1 and f 2 . Here the length, the magnetic field, and the velocity are, respectively, scaled by the radius of the outer sphere r o , the mean field strength B 0 and the MC wave speed where the aspect ratio η = r i /r o . As β → ∞ at s = 1, the governing equations are singular there; these boundary conditions ensure that the regular solution is selected.
Of particular interest is the regime when Le −1 is large. Taking the limit leads to a balance between the vortex stretching and the Lorentz term in the vorticity equation: whilst (2.2) retains its same form. The nonlinear problems have two source terms acting on the magnetostrophic wave: below we asymptotically solve the weakly nonlinear cases.
To seek solitary long-wave solutions we introduce slow variables with a small parameter ǫ (≪ 1) and a real constant c: Note that this assumes a long spatial scale in the azimuthal direction compared with the radial direction. This is reasonable for small m. We then expand variables with ǫ as for the basic state satisfying where D = d/ds. At zeroth order the equations of vorticity (2.4) and of electric potential (2.2), and the boundary condition (2.3) are all trivial.
At O(ǫ), (2.4) and (2.2) become respectively. Substituting (2.8) into (2.9) gives a homogeneous PDE with respect to g 1 : where L represents the linear differential operator comprising of s, ∂/∂s or D, B, β, U , and c. Inserting the boundary conditions (2.3) at this order into (2.9) yields We then seek a solution in the form of g 1 = Φ(s)G(ζ, τ ), so that Now the linear operator L is the ordinary differential operator with the partial derivatives with respect to s replaced by D. Given a basic state, the ODE (2.12) together with the boundary conditions is an eigenvalue problem to determine the eigenfunction Φ with eigenvalue c; it can have many eigensolutions. We note that the second-order ODE (2.12) remains non-singular as U /s → c, but not as B 2 /β → 0 unless s = 0. Below we concentrate on cases in which (2.12) has no internal singularities, i.e. there is a discrete spectrum. We consider cases where the z-averaged toroidal magnetic fields do not pass through zero (e.g. figure 3 of Schaeffer et al. (2017); figures 1-2 of Hori et al. (2018)). We proceed to the next order to obtain the amplitude function. Eqs. (2.4) and (2.2) at O(ǫ 2 ) yield (2.14) Eliminating ψ 2 using (2.13) and ψ 1 using (2.8), (2.14) becomes the inhomogeneous PDE where D 2 = (1/s)DsD. The boundary conditions here are ∂g 2 ∂ζ = 0 at s = η, 1 . (2.16) The adjoint linear problem corresponding to (2.10) is The adjoint boundary conditions are Note that the substitution B 2 Φ † /sβ = Φ reduces the adjoint problem to the ordinary linear problem (2.10) so, provided B 2 Φ † /sβ is non-zero in the sphere, the adjoint eigenfunction Φ † can simply be found by dividing the solution of (2.10) by B 2 /sβ.
The solvability condition to (2.15) is thus given by (2.20) Eq. (2.19) is the Korteweg-de Vries equation if the coefficients, α and γ, are both nonzero.
In the following section we examine the coefficients for different choices of the basic state. We note that the presence of U does not directly impact either α or γ. It however dictates Φ and Φ † through the linear problems at O(ǫ) and then may contribute to the terms at O(ǫ 2 ). This is in contrast with the hydrodynamic case (e.g. Redekopp 1977), where the basic flow enters the nonlinear term at O(ǫ 2 ) too. The mean-flow effect on the magnetostrophic wave arises from the equation for the magnetic potential (2.2).
Solutions to (2.19) may take the form of solitary (single or multiple soliton), cnoidal, similarity, and rational waves (e.g. Whitham 1974;Drazin & Johnson 1989). For instance, for a single soliton the asymptotic solution up to O(ǫ) is This is an eddy that has the solitary characteristics in azimuth, riding on the basic state with the linear wave speed. The finite-amplitude effect α accelerates the retrograde propagation if γ < 0, but decelerates it when γ > 0. The characteristic waveform is clearly visible in the magnetic potential.

Illustrative examples
We solve the eigenvalue problem (2.12) and the adjoint problem ( to as a Malkus field hereafter), the second B is inversely proportional to s (an electricalwire field), and the third one is (3/2) cos {π(3/2 − 50s/19)} + 2, which was adoped by Canet et al. (2014) to model a profile of the radial magnetic field B s within Earth's core (a CFF field). For the Malkus and wire fields the terms J in (2.12), (2.17) and (2.20) all vanish, whereas this is not the case for the CFF field. The Malkus field case has been extensively studied in the literature (e.g. Malkus 1967;Roberts & Loper 1979;Zhang et al. 2003;Márquez-Artavia et al. 2017). We also consider the inclusion of a basic zonal flow U that is prograde with either a linear or quadratic dependence on s. Table 1 summarises the results, listing the eigenvalue λ = |c|/2 (see below) and c for the n-th mode, the coefficients α, γ, and δ 0 as calculated from the eigenfunction Φ, the adjoint eigensolution Φ † and (2.20), and whether/at which s the wave speed c approaches the basic angular velocity U /s. Here the n-th mode has (n − 1) zeros within the explored interval. Negative values of c indicate retrograde waves. More notably, in the all cases we obtain nonzero α and γ for all n examined and so the KdV equations are appropriate. The fraction |α/γ| and their signs characterise the solitons.
For the Malkus field (B = s) and no mean flow U , we let x = 1 − s 2 and Φ(x) = xy(x) to rewrite the ODE (2.12) as where λ 2 = −c/4. This is a hypergeometric equation, which has a solution where F denotes the hypergeometric function (e.g. Abramowitz & Stegun 1965). The eigenvalue λ is determined by the condition Φ = 0 at s = η. The adjoint solution is related to the axial electrical current generated at this order as −D implying the current is nonzero at s = 1. Figure 1 shows the solutions in the Malkus case. Figure 1(a) shows profiles of B(s), the topography β, the eigenfunctions Φ for n = 1 and 2, and their adjoint eigenfunctions Φ † (3.2). This yields α ≈ −12.85 and γ ≈ 0.87 for n = 1; the nonlinear effect is more significant than the dispersive one. Figure 1(b) illustrates a single soliton solution (2.22) of ψ for n = 1. If the amplitude ǫ is too large, neglected higher order terms will be significant; if ǫ is too small the azimuthal scale of the solitary wave is too large to fit in, so we choose ǫ = 0.1 as a reasonable compromise. The streamfunction ψ is negative, indicating a clockwise solitary eddy. The retrogradely propagating vortex ψ 1 is slightly more concentrated at the outer shell than the magnetic potential g 1 (not shown). As c < 0 and γ > 0, the dispersion term reduces the retrograde propagation speed. We note that a clockwise vortex is observed in Earth's core (Pais & Jault 2008) and geodynamo simulations (Schaeffer et al. 2017): its implications are discussed in the final section.
The same basic states admit high-n modes with more isolated structure to have the KdV equations with nonzero α and γ ( Table 1). The speed |c| increases with n, confirming the dispersivity of the wave. The eigenfunction Φ for n = 2 is negative at small s, and then turns positive when s 0.787 (dashed-dotted curve in figure 1a), so the eddy is clockwise in the outer region and anticlockwise in the inner region ( figure 1c).
We next consider the basic field given by the wire field, B = 1/s, whilst U = 0. By using Φ(x) = xe λx y(x), (2.12) may be reduced to a confluent Heun equation The solution regular at s = 1 corresponding to the eigenvalue λ is where H c represents the confluent Heun function with the accessory parameter q c = λ 2 +2λ−1 and exponent parameters α c = λ 2 +3λ, γ c = 2, δ c = 1 and ǫ c = 2λ (Olver et al. 2010). This case admits a simple form of the coefficients (2.20) such that To evaluate the function we use the algorithm of Motygin (2018) below. Figure 2(a) gives profiles of the basic state and eigenfunctions. The figure shows that Φ for n = 1 has a peak nearer the outer boundary, compared with that for the Malkus field; it is still propagating retrogradely and is dispersive. This case yields α ≈ −36.9 and γ ≈ 1.25 for n = 1 and with ǫ = 0.1 the soliton is a more compact, clockwise eddy (figure 2b). Analysis of the individual terms of the coefficient α 0 in (2.20) implies that the presence of high order derivatives is favourable for nonlinear effects. For n = 2, dispersive effects are enhanced compared to nonlinear ones. The solitary eddy is clockwise in the outer region when s 0.894 and anticlockwise in the inner region (figure 2c).
To explore more general cases we implement the Matlab routine bvp4c to solve the eigenvalue problems. We retain the boundary condition Φ = 0 at s = η = 0.35, but use the modified condition Φ + (1 − s)DΦ = 0 close to the outer boundary s = 0.99999 to avoid the numerical issue arising from singularities when s → 1. We also impose a normalising condition DΦ at the inner boundary: the values for the Malkus field and the CFF field are given by (3.2), whereas the one for the wire field is by (3.4). The number of gridpoints in s is 500 in all cases. Given the obtained c, the same routine is adopted to solve the boundary value problems for Φ † . For consistency with the earlier cases we set Φ † = 1 at the outer boundary. The codes are benchmarked with the exact solutions. With modified boundary condition, our computational results match the expected eigenvalues λ = |c|/2 and eigenfunctions Φ for 1 n 3 with errors less than 0.01 % and 0.2 %, respectively. Now the third basic field, B = (3/2) cos {π(3/2 − 50s/19)}+2, is examined. Figure 3(a) depicts the basic state, the eigenfunctions for n = 1, and additionally J (represented by the red dotted curve). It is nonzero except at s ≈ 0.40 and 0.78 and is negatively peaked at s ≈ 0.59. The eigenvalues c do not differ from those in the Malkus case very much (Table 1). For n = 1, Φ has a peak at s ≈ 0.61 (blue dashed curve), as so does the basic field. This case gives α ≈ −11.5 and γ ≈ 2.85. Indeed the term including J dominates over the ODE (2.12) and also over α 0 (2.20); if the term Φ 2 D(J /β) were absent, α would become ≈ 1.68. Figure 3   Therefore it affects the speed c of propagation of the mode, whilst leaving its other properties unchanged (Table 1). For a more realistic flow, U = 4s(1 − s), with the Malkus field, the structures of Φ and Φ † are not drastically altered (leading to δ 0 ≈ 0.049). The dominance of the nonlinearity over the dispersion, |α/γ|, is however weakened. The presence of the same basic flow in the wire field case also exhibits this property. Finally, we comment on the behaviour of solutions in the vicinity of the point s at which U /s equals c, the location of a critical layer for the hydrodynamic Rossby wave soliton (e.g. Redekopp 1977). We impose a fast mean zonal flow, U = 80s(1 − s), in the Malkus field case; figure 4(a) shows the basic state and additionally the deviation from the wave speed, U /s− c (blue dotted curve). The curve shows this case has such a critical point at s ≈ 0.838. Nevertheless the impact is hardly seen in the eigenfunctions Φ and Φ † : there are no discontinuities in the derivative DΦ (figure 4b) and hence in the solitary wave solutions (figure 4c). This remains true for the wire field case with U = 320s(1 − s).

Concluding remarks
In this paper we have performed a weakly nonlinear analysis of magnetostrophic waves in QG spherical models with azimuthal magnetic fields and flows. The model we considered is an annulus model (Busse 1976;Canet et al. 2014) of the form utilised by Hide (1966) for linear magnetic Rossby (MR) waves. We found that the evolution of the long-wavelength, slow-MR waves in the spherical shells obeyed the KdV equation, whether the toroidal magnetic field and/or the zonal flow were sheared or not. The model we consider here is formally valid for cases where the azimuthal lengthscale is much longer than that in radius; the most obvious application of which is for thin spherical shells. For thicker spherical shells like those representative of Earth's fluid outer core, the ratio of these lengthscales is of the order ten. For thinner shells relevant to other astrophysical objects one might expect the asymptotic procedure to give a better approximation to the true behaviour. We find that solutions may take the form of a single soliton solution (for n = 1) which is a clockwise, solitary eddy when basic state magnetic field is any of a Malkus field (B ∝ s), a magnetic wire field (B ∝ 1/s), and a CFF field (comprising of a trigonometric function). In addition to these steadily progressing single-solitons we also find N -soliton solutions; as these satisfy the KdV equation we know that these may have peculiar interactions including a phase shift after a collision and FPU recurrence (e.g. Drazin & Johnson 1989).
We conclude by noting that inversion of the geomagnetic secular variation appears to detect an anticyclonic gyre in Earth's core (Pais & Jault 2008;Gillet et al. 2015;Barrois et al. 2018); it is off-centred with respect to the rotation axis and is believed to have existed for more than a hundred years. Moreover, DNS of dynamos driven by convection in rapidly-rotating spherical shells have exhibited the emergence of a large vortex which circulated clockwise and modulated very slowly (Schaeffer et al. 2017); in these simulations the averaged toroidal magnetic field tended to strengthen beneath the outer boundary. Our solution tentatively supports the idea that such an isolated single eddy should persist, while drifting on MC timescales of O(10 2-4 ) years. The long wave can be initiated through instabilities due to differentially rotating flows (Schmitt et al. 2008), due to thermally insulating boundaries (Hori et al. 2014), and due to the magnetic diffusivity (Roberts & Loper 1979;Zhang et al. 2003). The steadily drifting feature of the solitons should of course be altered during the long-term evolution when dissipation plays a role in the dynamics. The presence of dissipation may also alter the eigenfunction (Canet et al. 2014) and thus the detailed morphology of the soliton too.
We note an alternative to account for the eccentric gyre is a flow induced by, for example, the coupling with the rocky mantle and the solid inner core, as DNS by Aubert et al. (2013) had demonstrated. The issue ends up in a debate which has lasted for decades: does the geomagnetic westward drift represent the advection due to a large scale fluid motion (Bullard et al. 1950) or hydromagnetic wave motion (Hide 1966). We shall investigate these issues further, as well as the role of critical layers, by solving initial value problems in a future study.