The dynamics of magnetic Rossby waves in spherical dynamo simulations: A signature of strong-field dynamos?

We investigate slow magnetic Rossby waves in convection-driven dynamos in rotating spherical shells. Quasi-geostrophic waves riding on a mean zonal flow may account for some of the geomagnetic westward drifts and have the potential to allow the toroidal field strength within the planetary fluid core to be estimated. We extend the work of Hori et al. (2015) to include a wider range of models, and perform a detailed analysis of the results. We find that a predicted dispersion relation matches well with the longitudinal drifts observed in our strong-field dynamos. We discuss the validity of our linear theory, since we also find that the nonlinear Lorentz terms influence the observed waveforms. These wave motions are excited by convective instability, which determines the preferred azimuthal wavenumbers. Studies of linear rotating magnetoconvection have suggested that slow magnetic Rossby modes emerge in the magnetostrophic regime, in which the Lorentz and Coriolis forces are in balance in the vorticity equation. We confirm this to be predominant balance for the slow waves we have detected in nonlinear dynamo systems. We also show that a completely different wave regime emerges if the magnetic field is not present. Finally we report the corresponding radial magnetic field variations observed at the surface of the shell in our simulations and discuss the detectability of these waves in the geomagnetic secular variation.


Introduction
Observations of waves can provide us with information on many aspects of geophysical and astrophysical flows. An example is found in the study of Earth's atmosphere and ocean. The rotation of the planet gives rise to different wave modes including inertial, Rossby, and Kelvin waves (e.g. Pedlosky, 1979). They often appear in stably stratified environments, leading to a mixture with internal gravity waves. Tropical meteorology succeeded in distinguishing each wave mode in cloudiness data by performing space-time analysis in comparison with the linear theory of equatorial shallow water models (see Kiladis et al. (2009) for a review). This advances the knowledge of the individual wave modes and their roles in, for example, transferring energy and momentum.
It is hence quite natural to seek wave motions within the interior of the planet. The low-viscosity electrically-conducting fluid in the outer core is believed to host dynamo action that generates the global magnetic field. This generated field, combined with the rapid planetary rotation, can substantially influence the dynamics of waves in the core. The study of rotating magnetohydrodynamic (MHD) waves therefore offers another approach to planetary dynamo theory. The primary effect of the magnetic field is to split hydrodynamic modes into fast and slow modes. This provides a wide range of timescales -from days to thousands of years -on which waves in the fluid core can operate.
Geomagnetic secular variation and the core flow models deduced from it give evidence of wave motions in the core (see a recent review by Jault and Finlay (2015) and references therein). Axisymmetric modes have been seen in the core. The excitation of torsional oscillations (TOs) has become evident and is a plausible candidate for 6 year variations that are observed in core flow models and length-of-day (LOD) fluctuations (Gillet et al., 2010). This finding is used to infer the radial profile of the poloidal magnetic field within the core and to suggest a z-mean rms strength of approximately 3 mT. Buffett et al. (2016)  oscillations excited within a thin, stably stratified layer (e.g. Braginsky, 1999). A comparison with the predicted frequency would imply the thickness of the layer is approximately 130-140 km.
However, axisymmetric modes cannot reveal the azimuthal component of Earth's magnetic field, which may be considerably stronger than the poloidal component, so these attempts will be naturally extended to non-axisymmetric modes. A prominent feature of the geomagnetic variation is the westward drift on timescales of 300 years, which is clearly observed in the Atlantic hemisphere (e.g. Finlay et al., 2010). Recent geodynamo modelling successfully reproduced the spatial structure of the secular variation (Aubert et al., 2013). Revisiting a hypothesis of Hide (1966), Hori et al. (2015) (hereafter referred to as HJT15) demonstrated in dynamo simulations that these longitudinal drifts could be produced by the propagation of slow magnetic Rossby (MR) waves riding on mean flow advection. The advantage of their approach is that it did not specify the configuration of the background magnetic field, but computed it from a dynamo model. This enabled them to estimate a z-mean strength of the internal toroidal field of about 10 mT at a depth of 0:8 r core , where r core stands for the core radius. There is a rich literature on non-axisymmetric modes (e.g. Malkus, 1967;Zhang et al., 2003;Canet et al., 2014), but it mainly uses simple imposed fields chosen for mathematical convenience rather than geophysical relevance. Chulliat et al. (2015) analysed the geomagnetic secular acceleration in updated global models, such as CHAOS-5, including Swarm satellite data, and reported a 6-8 years westward drift of the equatorially anti-symmetric component. They attributed this to a fast MR wave excited in the thin stable layer. Since current satellite missions are increasing both the temporal and spatial coverage of data, a solid theory and methodology will be fruitful.
A related, but more theoretical, issue is what types of waves are found in strong-field dynamos, and how do they differ from the waves that occur in weak-field dynamos. We distinguish between strong-field dynamos, in which the inertial and viscous forces are small compared to the magnetostrophic forces, namely Coriolis, pressure, Lorentz and buoyancy forces, and weak-field dynamos, in which viscous or inertial forces play a significant role (e.g. Roberts and King, 2013). When the magnetostrophic forces are in balance, it is expected that Taylor's condition (Taylor, 1963) will constrain the configuration of the magnetic field generated by the induction process and then diagnostically determine the fluid motions (see Hollerbach (1996) for a review). Some key parameters are the Elsasser number K, quantifying the relative strength of the Lorentz force to the Coriolis force, and the Ekman number E, measuring the strength of the viscous force to the Coriolis force. In the limit of rapid rotation, E ! 0, the presence of magnetic field with K increasing to Oð1Þ could destabilise rotating convection, thicken the convective rolls, and lower their frequency (e.g. Chandrasekhar, 1961;Fearn, 1979;Jones et al., 2003); the appearance of these effects was found to be highly dependent on, for example, the background magnetic fields and boundary conditions (e.g. Zhang and Fearn (1993), Zhang and Schubert (2000), Jones (2015) and references therein). This led to a scenario of strong-field and weak-field dynamos. In strong-field dynamos the convection is influenced by the magnetic field, but the flow may nevertheless be quite columnar.
Convection-driven dynamo simulations, retaining inertia and viscosity, have provided some evidence to suggest that we are approaching strong-field regimes, as well as quasi-Taylor states.
Plane layer models for E 6 Oð10 À5 Þ have attained such regimes (Rotvig and Jones, 2002;Stellmach and Hansen, 2004;Hughes and Cattaneo, 2016). Spherical simulations for E ¼ Oð10 À4 Þ have reported the possible approach to a Taylor state (Aubert, 2005) but a rather minor impact on non-axisymmetric convective structures (Soderlund et al., 2012). The effect of the field on the flow seems to be model-dependent, as simulations with different boundary conditions and driving have increasingly demonstrated the influence of magnetic field on convective length scales (Sakuraba and Roberts, 2009;Hori et al., 2010;Hori et al., 2012) and subcritical behaviour (Sreenivasan and Jones, 2011;Hori and Wicht, 2013). This model dependence is known in linear magnetoconvection studies (see above). However, a clearer approach to the strong-field regime has been demonstrated recently by Yadav et al. (2016) as the Ekman number is reduced. Dormy (2016) shows that there is a relationship between the magnetic Prandtl number, Pm, and the Ekman number that must be respected to remain on the strong-field branch. Here Pm is the ratio of the kinematic viscosity to the magnetic diffusivity. Even at modest E $ 10 À4 , strong-field dynamos may be obtained if Pm is large enough, but if E is reduced, Pm can also be gradually reduced, so as shown by Yadav et al. (2016) even Pm $ 0:5 is large enough provided E ¼ 10 À6 . Teed et al. (2014) found that torsional waves were most clearly seen in regimes with strong-field dynamos. We therefore explore here whether slow MR waves are also a signature of a strong-field dynamo.
Slow MR waves are symmetric about the equator, and are quasi-geostrophic (QG) modes with a long wavelength in the zdirection (parallel to the rotation axis) and a short wavelength in the transverse direction (e.g. Malkus, 1967;Zhang et al., 2003). Consequently they are faster than non-QG rotating MHD (or MC) waves and hence there is a greater likelihood for their detection in geomagnetic data. Also, this class of waves emerges associated with rotating spherical convection. The magnetic mode on which we are focusing is preferred at the onset of magnetoconvection when magnetic diffusion is weaker than thermal diffusion (Hori et al., 2014). Slow MR waves can be also excited in a thin stable layer, in which they generally travel eastward: we refer to Márquez-Artavia et al. (2017) for a comprehensive classification of linear waves in shallow water models. This paper extends the investigation of MR waves in spherical dynamo simulations, in which magnetic fields are selfconsistently generated. The aims are threefold. (i) Guided by the previous study (HJT15), we present more cases in which we were able, or unable, to identify the wave modes by performing spacetime analysis of the output data. The longitudinal drifts observed in the radial velocity match very well with the predicted wave speeds. (ii) Of particular interest are the dynamics of these waves: whether the identification could indeed represent a predominant magnetostrophic balance, and to what extent assumptions required for the wave theory could be appropriate. (iii) In the light of the analysis of the internal dynamics, we examine whether these wave motions could be detected in data of the magnetic field that is inferred at the top of the core. In Section 2, we present the mathematical formulation for our numerical models and the wave mode. The results are detailed in Section 3. Section 4 summarises our findings and we discuss implications for the study of planetary dynamos.

Formulation
We numerically model convection and magnetic field generation in a rotating spherical shell filled with an electrically conducting fluid. For applications to Earth's core, we adopt the Boussinesq approximation for an incompressible fluid. The details of our models are described in Teed et al. (2014) (hereafter referred to as TJT14), and we give only brief details here. The governing equations for temperature, T, velocity, u, magnetic field B, and pressure, p, are solved in a dimensionless form: @T @t ð4a; bÞ withê z andê r being the unit vectors in the z-and r-directions, respectively. The equations are scaled by taking the shell thickness D ¼ r o À r i for length, the magnetic diffusion time, D 2 =g, for time, ðql 0 gXÞ 1=2 for magnetic field, and D 2 =g for temperature. Here r i and r o are the inner and outer core radii, respectively, g is the magnetic diffusivity, m is the kinematic viscosity, q is the density, l 0 is the permeability of free space, X is the rotational angular velocity, and is the internal sink rate. We assume a volumetric sink term in the temperature equation for modelling compositional convection, as well as its boundary conditions of zero heat-flux at r ¼ r o and a prescribed heat-flux at r ¼ r i such that the energy contained in the fluid region is conserved. Other boundary conditions are assumed to be no-slip, electrically insulating, and co-rotating. The fundamental parameters are the Ekman number, E, the Prandtl number, Pr, the magnetic Prandtl number, Pm, and the Rayleigh number, Ra, which are defined as respectively. Here j is the thermal diffusivity, a is the thermal expansivity, and g o is the reference gravity at the outer boundary. We assume that gravity increases linearly with radius.

Theory
Rossby waves, whether they are hydrodynamic or MHD, are derived from the equation of vorticity n ¼ r Â u.
where r 2 h f ¼ 1 s @ @s s @f @s þ 1 s 2 @ 2 f @/ 2 and r h Á A ¼ 1 s @ @s sA s þ 1 s @ @/ A / for any vector field, A. The integral in N C is performed by using the sloping boundary conditions, u z ¼ Çu s s=H at z ¼ AEH. We used r Á n ¼ r Á J ¼ 0 as well as the solenoidal conditions (4a,b).
To seek perturbations about a background state, we split the velocity and magnetic fields into their mean and fluctuating parts. Furthermore, to focus on the background state given by the axisymmetric component, we further separate the mean parts into axisymmetric and non-axisymmetric parts, such that u ¼ U The averaging operators and fluctuating parts appearing here are defined by f ðt; s; zÞ ¼ 1 2p Substituting (10) into the Lorentz term, N L , we find its individual terms: Up to this point, everything is exact and no assumptions about the relative magnitudes of the different components of the flow and field, or the length scales on which they vary. However, the equations are very complicated, and to get a system which we can understand we must make assumptions about the relative sizes of the various terms. We start by linearising the fluctuating parts, i.e. consider only terms of first order in the primed quantities. We assume that the zero order quantities describe a slowly evolving flow and field state, and that the first order terms describe relatively fast wave motions perturbing that quasi-steady state. We also ignore terms which are second order in the fluctuating primed quantities, though as we see later, in actual simulations nonlinear effects are visible. Next, we assume the azimuthal length scale of our disturbances is short compared with the variation in the s and z directions, and short compared with variations of the mean quasi-steady flow and field. Of the terms in (13), the second is of second order, and the fourth and fifth are small under our length scale assumptions. We eliminate the third term by assuming that the axisymmetric part of the mean azimuthal field is bigger than the non-axisymmetric part. We are left with the first term on the right-hand-side, h e B Á rj 0 z i, representing the restoring part for MHD waves with respect to the background field e B. For modes with reasonably large azimuthal wavenumber, say m P 5, these assumptions are approximately true. We view the theory based on them as an 'ideal' theory to serve as a starting point, and we can then explore how the actual simulations depart from this idealised model. Similarly, the Reynolds term N R can be expanded as Here the first term, on the right-hand side, h e U Á rn 0 z i, describes the advection effect due to the background mean flow e U , which under our assumptions is the dominant term.
Separating the restoring and advective terms from the remaining terms, we rewrite the vorticity equation as z i denote the residual parts of the Reynolds and Lorentz terms, respectively. The Coriolis term still involves the mean and fluctuating parts; as the mean component is negligible in simulations, we omit this component hereafter. In the same manner, we rewrite the induction Eq.
(2) as where I S ¼ B Á ru À e B Á ru 0 and I A ¼ u Á rB À e U Á rb 0 .
The terms on the left-hand sides of (15) and (16) give the basic equations for MR waves (HJT15). The equations are linear with respect to fluctuating variables, n 0 z and j 0 z , but are coupled to each other, and exclude nonlinear terms including b 0 Á rj 0 z in the Lorentz force and u 0 Á rn 0 z in the Reynolds force. To elucidate the fundamentals of the wave modes, we first concentrate on the linear aspects, bearing in mind that we are assuming the background field and flow have axisymmetric azimuthal components which are at least comparable to the other components, and that the azimuthal wavelengths of our modes are short compared to the radial and axial components. Note that this is not obvious in every case and other components possibly become significant, as we shall discuss later. However, these assumptions do surprisingly well because the convection driving the modes in the models consists mainly of tall thin columns. We then operate d dt ¼ @ @t þ h e U / i s @ @/ over the left-hand side of (15) to obtain Substitution of the left-hand side of the induction equation (16) into this and using our length scale assumptions gives For some simple fields this equation can be solved analytically (see Canet et al. (2014) for detailed analysis). We instead suppose that u 0 s is approximately geostrophic, so that u 0 s ðHÞþ u 0 s ðÀHÞ % 2hu 0 s i, and that the radial gradient of the axial vorticity, n 0 z , is smaller than the azimuthal gradient, consistent with our previous assumptions, i.e. n 0 z % À 1 s @ @/ hu 0 s i. This is valid only for reasonably large m components, but it considerably simplifies the problem leaving Here we seek solutions with a form of hu 0 s i $ e ıðm/ÀxtÞ at given s and obtain the dispersion relations of the fast and slow modes as where the Rossby, Alfvén, and advection frequencies arê respectively. We see that the wave frequency is the sum of the dynamical wave frequency plus an advective term due to the mean flow.
The fast modesx þ essentially recall the hydrodynamic Rossby waves, which travel prograde with the frequencyx R about the advection part. They arise from a balance between the first two terms of (19), dhn z i=dt and N C . By contrast, the slow modesx À are a unique solution of rotating MHD, sometimes called MR waves or MC-Rossby waves. Their properties become evident when taking the limitx 2 M =x 2 R ( 1 on the slow mode,x À , to obtain (using the binomial approximation) and the observed frequency will be the sum ofx MR and the advection frequency x adv . This implies a much lower frequency and a retrograde propagation unless the advective flow is large and eastward. The corresponding phase speed is given V MR ¼x MR =m, and similarly for the Rossby and Alfvén phase speeds. The magnetic Rossby speed goes up as the wavenumber m increases or the radius s decreases. A balance between the last two terms, N C and N L , is vital for this mode, indicating that the time variations arise from the induction equation while the momentum equation is almost in balance. These slow waves will be distinguished from Alfvén or Rossby (fast MR) modes in terms of dispersion relations x ¼ xðmÞ, phase velocity x=m, and vorticity balances.
At fixed s and hence h f B 2 / i, all dispersion relations (20) are comprised of MR branches at lower wavenumber m and Alfvén branches at higher m. The transition will occur when x 2 We did not observe signals of Alfvén branches in our simulations, but it could be possible if faster or smaller-scale disturbances are provided, for instance, by more vigorous convection. Studies of equatorial atmospheric dynamics demonstrate an impressive ability to distinguish several wave modes through space-time spectra and theoretical dispersion relations (e.g. Kiladis et al., 2009).
Our assumption of a short azimuthal length scale means terms involving f B / dominate over the terms involving the poloidal field, f B s and f B z . We speculate that if these terms do become significant, the dispersion relation would become almost proportional to m. However, solving the linear equations in this case becomes difficult. Applying the assumption n 0 z % À 1 s @ @/ hu 0 s i helps to simplify our equation considerably. To pursue analytical solutions when all the field components are relevant, we are required to make further assumptions, such as uniform e B s and constant H; this would not give expressions, (20) and (19), for an azimuthal field.
Whereas the consideration of the restoring forces predicts the eigenmodes, an excitation mechanism determines what frequencies and/or wavenumbers indeed set in. Any terms appearing on the right-hand-side of (15) can initiate disturbances leading to wave motions. In our simulations, excitation is mostly created by the instability driven by the buoyancy N B . This is supplied everywhere at the inner boundary r ¼ r i . Topographic Rossby waves naturally occur, associated with convection in rotating spheres and spherical shells (e.g. Busse, 1970): thermal Rossby waves,x TR ¼x R =ð1 þ PrÞ, are preferred in the hydrodynamic case.

Numerical models
The models explored in this study and their global properties are listed in Table 1. The control parameters range over 1 6 Pm 6 5 and 5 Â 10 À6 6 E 6 10 À4 , while Pr ¼ 1. In five of the runs, the Rayleigh number is fixed at 8:32Ra c where Ra c denotes the critical Rayleigh number for the onset of nonmagnetic convection. These are the runs selected from the previous study by TJT14 and partly analysed by HJT15; in this paper we shall present a detailed analysis of the models. We also add two new runs for Ra ¼ 16:6Ra c to investigate the effects of higher Ra. Unlike axisymmetric TOs, the nonaxisymmetric waves are closely linked to the thermal instability and hence can be affected by the convective vigour. At the Pm regime explored here, the slow MR modes propagating retrograde are favoured at the onset of magnetoconvection, whereas other diffusive Rossby modes prevail at lower Pm (Hori et al., 2014).
Monitoring the time evolution of the kinetic and magnetic energies, we confirmed that each model reached a quasi-steady state and then we chose a short time interval, s, to analyse its time variation. The intervals s are 0.01 magnetic diffusion times for most models and are taken longer, s 6 0:02, for the large-E models 4R2 and 4R5. By equating B s at the CMB from our simulations to its known value from the geomagnetic data the time intervals can be translated to the dimensional time, s E (TJT14): all our analyses presented below correspond to s E less than 83 years (see HJT15).
In our simulations magnetic fields are self-consistently generated. They are overall dominated by axial dipoles that do not reverse during the time intervals. These stable large-scale fields act as the background field (such as e B) for the disturbances (u 0 and b 0 ) discussed below. The morphology of the background fields will be presented in Section 3.1. To characterise each run for the magnetostrophic wave motions, we pay attention to the following output parameters, defined in our scaling as The Elsasser number, K, measures the relative strength of the Lorentz to Coriolis forces. The smallness of the Taylorization parameter, T , indicates to what extent the system resembles a pure Taylor state. This parameter increases with s, as reported by Wicht and Christensen (2010), and thus suggests a better Taylorization nearer the rotation axis. The parameters U C and U 0 C quantify the geostrophy of fluctuating zonal flows with respect to the total flows and the fluctuating parts only, respectively. The latter indicates the dominance of axisymmetric TOs on short timescales. Investigating extensive magnetoconvection runs, Teed et al. (2015) found that TOs were identified when the parameter U 0 C J 0:4. Additionally, since the non-axisymmetric motions of the cylindrically radial velocity are also of interest, an equivalent geostrophy parameter U s C defined with the radial component, u 0 s , is also to be checked. The values of these quantities for each run in our suite of simulations appear in Table 1. Other output parameters for lower Ra models are found in Table 1 of TJT14. The magnetic Reynolds number of all our runs, Rm ¼ ffiffiffiffiffiffiffiffi ju 2 j p , ranges from 100 to 450. The Rossby number, Ro ¼ RmE=Pm, hence remains no greater than 0.001. The small Ro is consistent with the observations of the stable dipolar field solutions (Christensen and Aubert, 2006).
We recall the classification made by TJT14 and find strong-field solutions for all the presented models except 4R2. These strongfield dynamos show K greater than unity, T less than 0:2, and relevant TOs detected (Table 1). For the non-axisymmetric dynamics, we define measures for the length scales in the kinetic power spectrum: the mean harmonic degree and the peak harmonic order, m peak , i.e. the value of m for which ju m Á u m j=ju Á uj; is greatest, summing over all possible '-values for that particular m. These values are expected to remain large for weak field regimes but be smaller for strong-field regimes (Section 1). The mean value ' has often used in analysis of recent dynamo simulations. The dependence of ' on E appears to retain the nonmagnetic scaling of E À1=3 [not shown; e.g. Roberts and King (2013)]. This may not be a good measure when the spectrum has several peaks indicating more than one distinct scale. Also, spectra with respect to the harmonic order m, rather than ', better represents the convective structure in rapidly rotating spheres and spherical shells. The peaks m peak hence indicate the enlarging effect, depending on the field strength.
The influence of the generated magnetic fields on the flows becomes evident when comparing results with the corresponding nonmagnetic simulation and evaluating the force balance (as shown below). These magnetic effects are hardly found in model 4R2, so we term this model a weak field solution. We used the Leeds spherical dynamo code to solve the full equations, (1)-(4a,b); see Willis et al. (2007) and Jones et al. (2011) for a detailed description. In the code, a predictor-corrector method is adopted for choosing timestep sizes, the longitudinal and latitudinal grids are expanded in the spherical harmonics, and the radial grid uses a finite difference method. The number of grid points in the r; h, and / directions was N r ¼ 160; N h ¼ 288, and N / ¼ 576 for most runs, respectively, but needed to be increased up to N r ¼ 192 for low-E or high-Ra models. Here the h and / resolutions were given by the maximum spherical harmonic Table 1 Control parameters and global properties of our dynamo models. Prandtl number Pr ¼ 1 throughout. K; T ; UC ; U 0 C and U s C are defined in Eq. (23), and ', m peak in Eq. (24) and below. The column m peak presents the peak modes as well as the strongest secondary modes in order. Column TO denotes whether torsional oscillations were found or not (Yes/No). Results for nonmagnetic convection are given in parentheses.
The output data were transformed to cylindrical polar coordinates for comparisons with the QG theory. The resolution was reduced for post-processing, for which grid points in the s and z were typically fixed at N s ¼ 128 and N z ¼ 96, respectively.

Predicted wave and advection speeds
In Fig. 1a, we compare phase speeds of slow MR, V MR ¼ jx MR =mj, Alfvén, V M ¼ jx M =mj, and nonmagnetic Rossby waves, V R ¼ jx R =mj, as a function of normalised radius s=r o for model 4R5, using formulae (21) and (22). The magnetic modes, V MR and V M , were calculated using the z-mean toroidal field, iðsÞ, from the simulation. The nonmagnetic speed V R , and the nonrotating one V M , are much greater than V MR so they have been rescaled down. This time-scale separation indeed helps to distinguish each wave mode. The Alfvén speed, V M , plotted with black solid curves, is a proxy for the profile of the background toroidal field: the structure of f B / is presented in Fig. 2a. The field component is strengthened beneath the CMB at the equator, as commonly seen in spherical dynamo simulations (e.g. Christensen and Wicht, 2015). This implies the V M profile is fastest at s % 0:9r o . The blue solid and dotted curves plot V MR and V R for a chosen wavenumber of m ¼ 5, respectively. This wavenumber is chosen because it gives the clearest image of the MR waves in the simulations, see below. The MR speed, V MR , becomes slower with increasing s, as expected from (22), with large values near the TC. The speed of the waves is similar to that defined in the linear analysis of Zhang (1995). Waves are quite geostrophic, travel westward, and their frequency increases with wavenumber m.
These wave motions will be observed, riding on the geostrophic mean zonal flow, with f ¼ h f Fig. 1b shows the profile f ¼ fðsÞ for the same model, including the sign. The zonal flow is prograde near the inner boundary, when s K 0:45, but retrograde at an outer radius. So at the middle of the shell (s=r o ¼ 0:5) the background flow becomes extremely slow. Comparing this with thex MR profile allows us to ascertain at which radius wave propagation will dominate over mean flow advection. To explore the wave dynamics, we choose the mid-depth, s ¼ 0:5r o , for analyses presented in the following subsections. Similarly, Figs. 1c-d demonstrate wave and mean flow speeds for a low-E model, 6.5R2. The profiles are similar to those in the earlier run; there are differences in the details, such as the maximum speeds and the radius at which f changes sign. Fig. 2b illus-trates that the field f B / outside the TC is concentrated at low latitude. Analogous plots are found in all the models except 4R2, so we avoid presenting other plots.
Figs. 1e-f and 2c depict the exceptional case, 4R2. In this case the magnetic modes, V M and V MR , are orders of magnitude slower than those in other runs, indicating that the background magnetic field is rather weak here. This makes the bifurcation from the Alfvén to MR waves less drastic (see the following subsection). The s-profiles indicate that the morphology of the background field is more complex than others: the field is found to hold two wavenumbers in s outside the TC (Fig. 2c). The flow profile, f, is also remarkably distinct; it is retrograde for all s. The distinction in the f structure is reminiscent of the work by Aubert (2005), who discussed the influence of magnetic fields on axisymmetric zonal flows. Indeed the generated field hardly affects the zonal flows in this particular model: an analogous profile in a nonmagnetic model is shown later in Fig. 7.

Space-time analysis of internal radial velocity
The wave equation (19) gives a description for the z-averaged radial velocity, hu s i, which is the variable to be analysed in this subsection. Fig. 3a displays a snapshot of the spatial structure of hu s i, in the view from the northern pole, for model 5R5. The presence of the strong field with K % 22 fattens the convective struc-tures here (cf. the nonmagnetic case in Section 3.4). The raw (i.e. not averaged) radial velocity, u s ðzÞ, sliced at the equatorial plane is very similar to hu s i. This is confirmed by checking meridional slices of u s ðzÞ (Fig. 3b); columnar (i.e. z-independent) structures are found even for the very strong magnetic field. The geostrophy parameter U s C for the cylindrically radial velocity amounts to 0.35 and ensures the dominance of the geostrophic component in the whole flow. Also, the equatorial plots show that the azimuthal gradient therein is steeper than the radial one. Therefore a key assumption for the theoryj @ @s hu 0 / ij ( j 1 s @ @/ hu 0 s ij -leading to n 0 z % À 1 s @ @/ hu 0 s i, is found to be appropriate. Fig. 3c-d, for model 6.5R2Ra, demonstrates similar slices to confirm the high m approximation on n 0 z and the two-dimensionality of the flow. The moderately strong field, K % 6, enlarges azimuthal scales,  compared to the corresponding nonmagnetic case, but they remain rather small for this lower E model. Fig. 4 shows time-azimuthal sections of hu 0 s i at s ¼ 0:5r o for runs 4R5, 5R5, 6.5R2, 6.5R2Ra, and 4R2. The left column displays plots in the physical domain (i.e. t-/ space), and plots in the spectral domain (i.e. f-m space, where f ¼ x=2p) are shown on the right.
To calculate the spectra, we performed two-dimensional FFTs of hu 0 s i at the chosen s. These spectra are important for comparing with the predicted dispersion relations, but also for determining the dominant m or f components. With the wavenumbers being determined, we calculate the respective phase speedsx MR and compare them with observed longitudinal drifts. The chosen m for each model is presented in the figure. In the physical domain, white dashed lines draw the advection speed f (at the chosen radius) and black solid lines indicate the combined phase speed, f þx MR =m, for the selected m. Since the background flows at the mid-depth are very slow, the white lines appear to be almost vertical in all cases. No black lines are shown for model 4R2, for which we do not find MR waves. Analogously in the spectral domain we continue to use white dashed lines for the advective dispersion relations, f ¼ x adv =2p ¼ fm=2p, and black solid curves for the total one, x ¼ ðx adv þx AE Þ=2p. In the strong field models the fast modes, x adv þx þ , are far off the frequency window; this branch appears for f > 0 only in the exceptional model, 4R2. In the same spectra, we also indicate the Alfvén modes, f ¼ ðx adv AEx M Þ=2p, by white solid lines; these are linear in m.
Model 4R5 -for large E and very large K (% 18) -illustrates wave identification most clearly ( Fig. 4a-b). From the spectra we find m ¼ 5 and 3 modes excited significantly. Migrations in t-/ space almost perfectly fit with the calculated total phase speeds for m ¼ 5. As the convective rolls in the model spread throughout the radius (not shown), these wavenumbers are dominant at any s. Here recall that the theory assumes that the azimuthal scales are smaller than the radial ones; we hence exclude the lower wavenumber mode m ¼ 3 for the identification. A lower E model, 5R5, displayed another successful identification for m ¼ 5 and 8 ( Fig. 4c-d). The crests and troughs observed there were narrower and sharper than those in the larger E model. Models for identical Ra=Ra c but smaller Pm, 5R2 and 6.5R2, yield weaker generated fields of K % 2 and larger m are significant. Fig. 4e-f demonstrates the plots for model 6.5R2. For the weak background field the dispersion relation, x adv þx À , predicts a slower wave speed. The spectral analysis shows a strong signal of ðm; f Þ % ð9; À300Þ; however the frequency is higher than that of slow MR waves and too low for the Alfvén waves. There are some features that travel at the m ¼ 9 MR phase speed (see Fig. 4e), but there also features travelling at different speeds. The signals for ðm; f Þ % ð9; À100Þ; m ¼ 8 and 12, may be interfering with the m ¼ 9 mode to give a more confused picture than in Figs. 4a and c. The migrations are very slow, but even though the phase speed is not so well-defined, the sharp wave forms are found to be persistent.
In the higher-Ra moderate-K models, 5R2Ra and 6.5R2Ra, we also see more complex drift patterns. Fig. 4g for model 6.5R2Ra shows that the duration of the migrating crests and troughs becomes shorter. Vigorous convection gives rise to more chaotic motions and hence interrupts wave propagation more frequently. Nevertheless, we are able to find signals distributed over the predicted dispersion relation x adv þx À (Fig. 4h). Note that the advective velocity f at s ¼ 0:5r o is positive for this run. This may explain the prograde drifts seen in real space (white dashed lines). The total phase velocity for the preferred m ¼ 9 mode remains retrograde and gives a correct speed that matches slow retrograde drifts. However the significant signal of m ¼ 7 and f > 0 cannot be met with any of the dispersion relations shown in the spectral domain. This indicates the limitation of the present theory; we may here be seeing diffusive MR waves that can propagate prograde (Hori et al., 2014). In larger-E model 5R2Ra, larger azimuthal scales (m ¼ 4; 7, and 9) are selected and prograde migration is less clear.
Finally, the weak field model, 4R2, demonstrates a failed case ( Fig. 4i-j). At these parameter regimes, the present setting for fixed heat-flux boundary conditions can cause a mixture of very wide convective rolls, such as m ¼ 1, and rotationally-constrained thinner ones (Hori et al. (2012); also see later in Section 3.4). This results in a hemispherical structure seen here in the physical domain. The spectral analysis in Fig. 4j shows no relevant signals along the MR dispersion relations except at m ¼ 1 and 2. They are instead better aligned with the Alfvén modes; however we exclude this because the generated field is weak here. To host Alfvén waves, a requirement is for the system to satisfy the very strong-field limit x 2 M =x 2 R ) 1. The force balance, as presented in the next subsection, shows a minor role for the Lorentz force in this run. This is consistent with the fact that axisymmetric TOs were also not identified in this dynamo model (TJT14). Table 2 summarises some properties of the non-axisymmetric motions of all the runs, all taken at s ¼ 0:5r o . Column MR indicates whether magnetic Rossby waves were detected: only run 4R2 failed to show any. The value of m is determined by finding the largest peak sitting around the dispersion relation in the wavenumber-frequency power spectrum (right hand panels of Fig. 4) and V MR is the corresponding phase speed, which can be compared with the advective phase speed f. In all cases V MR is larger, showing that at this s-value migration is mainly due to wave motion rather than advection by a mean flow. V rel M is the relative strength of the internal azimuthal field to the radial field as measured by the ratio of Alfvén waves at s ¼ 0:5r o . Note that in these dynamo models the radial field is stronger than the azimuthal field. We don't currently know whether this holds for the actual field in the core.
A striking feature of the observed MR waves are their waveforms because they do not show wave packets, but rather feature isolated crests and troughs. This is surprising as the highly dispersive waves (22) may be expected to form wave trains comprised of several m components. To show more details, Fig. 5 depicts the evolution of the amplitude hu 0 s i for model 5R5 (as shown in Fig. 4c). Given a disturbance, it grows to a crest or trough, whilst travelling retrogradely. Meanwhile, waveforms steepen and shift to the positive side: for instance, a crest peaked when t ¼ 0:005 between p=2 < / < 2p=3. These are reminiscent of steepening, particularly of cnoidal or solitary waves, which are typically known in the weakly-nonlinear dynamics of inviscid, dispersive waves (e.g. Whitham, 1974). The theory suggests that the effects will be more relevant as the system becomes inviscid: this agrees with our observation that lowering E produced sharper waveforms. This indicates that the nonlinear terms, which we omitted for the theory, are important in creating the observed wave patterns, while the linear part is fundamental to determine the wave speeds; we shall address this in Section 3.5.

Vorticity balance
To elucidate the nature of the MR waves, we evaluate the individual terms of the z-averaged vorticity equation (7) in terms of the Table 2 Properties characterizing the nonaxisymmetric motions at radius s ¼ 0:5ro for our models. For each run a preferred wavenumber, m, its MR-speed, VMR, and the advection speed, f, are presented. The relative strength of the internal azimuthal field to the radial field is measured by the ratio V rel M of Alfvén waves at s ¼ 0:5ro.  Fig. 6 depicts time-azimuthal sections of those terms at the same radius for the model 5R5. In every plot we retain the white and black lines from Fig. 4 to mark the predicted phase speeds. We also use identical colour contour steps for every plot with the maximum of the individual terms listed in Table 3. Fig. 6, for model 5R5, illustrates that the vorticity equation is dominated by the Lorentz force, N L , and Coriolis force, N C , terms.
Other terms such as the inertia, @hn 0 z i=@t, Reynolds force, N R , and viscous force, N V , can become relevant locally and temporarily.
Their amplitudes remain smaller than those of the two dominant terms (Table 3), indicating their minor roles throughout the time evolution. The significance of N C and N L agrees with the fact that this model nicely demonstrated propagation of the slow magnetostrophic waves. We recall that the analysis here is made for the geostrophic component. This reveals a predominant dynamical balance between the Coriolis and Lorentz forces within the QG approximation. This confirms former findings in linear rotating magnetoconvection (Zhang, 1995): it is now seen in nonlinear dynamo systems.
The buoyancy term, N B , at this radius is weaker than the other contributions, as can be seen from the values in Table 3, as well as the amplitude in Fig. 6. This term is most significant at the inner boundary, at which the buoyancy source is set. The disturbances arising from the bottom spread towards an outer shell and induce the longitudinal wave motions at a given s. Therefore, in spite of its small magnitude at mid-depth, the time-azimuthal patterns are found to almost perfectly correlate with those of N C and N L . The buoyancy term is therefore crucial for driving the observed wave motions.
Primary roles of N C and N L are found in all models except 4R2. Whereas the magnetically dominated run 5R5 yields a sizable N L , the two terms were almost in balance for moderate-field models, such as 5R2Ra. However Table 3 shows that other terms can become significant locally. Despite the clear wave identification, the large-E model, 4R5, has significant contributions from @hn 0 z i=@t; N R , and N V . Lowering E helps to suppress these terms; this is crucial for steepening waveforms. Table 3 also shows that a higher Ra seemingly increases the significance of N R . For some runs, N R can occasionally become comparable with N C , but it is extremely localized when it does so. This indicates that magnetostrophic balance remains dominant most of the time.
This balance does not hold for model 4R2, in which MR waves were not identified. The Lorentz term, N L , is weaker by an order of magnitude, and instead N R ; N V , and @hn 0 z i=@t are stronger (Table 3). We thus confirm only a minor role for the magnetic field in this model, and exclude the excitation of Alfvén waves. One may expect that a weaker field could host fast MR modes, or nonmagnetic Rossby waves. However, we do not find any direct evidence of such waves. For the fast wave motions, a predominant balance between @hn 0 z i=@t and N C is mandatory. The significant magnitude of N R and N V in the model suggests that this is not the case.

Hydrodynamic model
To make clear the impact of magnetic fields, we explore the corresponding nonmagnetic models where the induction equation is not solved and hence magnetic field generation is switched off. Fig. 7 displays a snapshot of the non-axisymmetric structure of hu s i for a run with E ¼ 10 À5 , termed NM_5R5. In the absence of the magnetic field, convective rolls overall get thinner in azimuth

Table 3
The maximum of each term of the vorticity equation, (7) and (8), where the two most significant terms for each model are indicated in bold. At radius s ¼ 0:5ro. Results for nonmagnetic convection are given in parentheses.
In Fig. 9 we evaluate each term of the vorticity equation for this nonmagnetic model. We see that N R is as significant as N C and @hn 0 z i=@t, so that although the wave speed is primarily the thermal Rossby wave speed the nonlinear Reynolds stress is affecting the waveforms.

Restoring force and nonlinearity
We have seen that the formula for toroidal field, given by (20), is able to account for some of the observed longitudinal drifts.
Meanwhile, the poloidal component, f B s and f B z , possibly acts as a restoring force. To quantify this, we measure the ratio, V rel M , of Alfvén waves, V M , for the azimuthal component to those for the radial component, Here U A is equivalent to the propagation speed of TOs. Table 2 lists the relative speeds at the mid-radius and shows that the radial field components are stronger for all the models except 4R2. Indeed, in standard dynamo simulations, the axisymmetric poloidal field is found to be equal to or stronger than the toroidal one (e.g. Christensen and Wicht, 2015). Note that the relative strength in the Earth's core is unknown; some estimation has been made [Zhang and Fearn (1993), Shimizu et al. (1998)  . Spatial structures of (a) z-averaged radial velocity husi (cf. Fig. 3a) and (b) angular velocity f (cf. Fig. 1b of HTJ15) for nonmagnetic model NM_5R5.   (Fig. 4c) and in the Coriolis part of the restoring force (Fig. 6c). We therefore display a filtered Fig. 10a in Fig. 10d, and see that the pattern visible in Figs. 4c and 6c has reappeared. This suggests that for wave motions of the preferred wavenumber mode, m ¼ 5, the toroidal field has primary importance. It is, however, quite possible that the radial background field can have some influence over the wave speed, particularly for lower values of m.
The observation of wave steepening, the surprisingly thin wave fronts seen in the left panels of Fig. 4, implies a considerable nonlinear effect on the amplitude (Section 3.2). In the linear theory we omitted two types of nonlinear terms in the vorticity equation: Lorentz, hb 0 Á r H j 0 z i À hj 0 Á r H b 0 z i, and Reynolds, hu 0 Á r H n 0 z i À hn 0 Á r H u 0 z i, terms. Evaluation of these two terms shows that the maximum of the nonlinear Lorentz term is orders of magnitude greater than that of the Reynolds term at any chosen time (not shown). Indeed, the nonlinear Lorentz term is equivalent in magnitude to the restoring part. An interesting question here is whether only a limited number of terms from N L can model the pattern of N C , or hu 0 s i. In Fig. 11 we test this by taking a sum of The selected terms reproduce some features including sharper crests and troughs. We hence speculate that, although the linear theory is essential for explaining its wave speeds, the nonlinear Lorentz term is important for creating the observed waveforms. This will help us to study the fundamentals of the nonlinear dynamics, for example, by adopting reduced models.

Space-time analysis of surface magnetic field
We now address the question whether MR waves could be detectable in geomagnetic data. The westward drift is analysed using the radial component of the geomagnetic field, which is inferred at the top of the core (e.g. Finlay et al., 2010). The QG theory, when no boundary layers are taken into account, suggests that the internal wave motions at given s can be seen at the top at latitude % arc cosðs=r o Þ in each hemisphere. Therefore one may expect identification of MR waves in the secular variation if the flow is sufficiently two-dimensional. Note that the geostrophy varies with the Ekman number E and the background magnetic field, which can be quantified by the Elsasser number K. Fig. 12 depicts plots for space-time analyses of the radial magnetic field B r observed at the outer boundary r ¼ r o in model 6.5R2, in which low E and K % 2 give a well-defined geostrophy. These are analogous to the plots shown of the internal fluid motions discussed in Section 3.2. To focus on the secular variation, we remove The spectrum at 60 N is dominated by signals of m % 9 and 12 and f < 0; prograde modes of f > 0 also look significant. The predicted wave speed for m ¼ 9 can fit some magnetic drifts observed in the physical domain. At lower latitude 39 N drift patterns seemingly get noisier. As jfj goes up and V MR does down as s increases to 0:77r o (see Fig. 1c-d), so flow advection becomes more relevant here. A higher m of 15 increases the contribution due to wave propagation, and this can be distinguished from the contribution due to advection. However, the spherical harmonic components of the geomagnetic field with m > 12 are hard to detect due to crustal field contamination, so these higher wavenumbers will not be easy to identify. In Figs. 12e and f, we further test this detectability by excluding all the wavenumber modes when m > 12 from the magnetic data at each latitude. The filtered plot at 60 N retains the wave patterns identified in the whole data in Fig. 12a. Some drifts at 39 N remain visible when filtering, but they run almost parallel to the advection speed here.
Figs. 12g and h display t-/ sections at 60 S an 39 S in the southern hemisphere. When the flow is quasi-geostrophic we expect the B 0 r signal in the southern hemisphere to be the same as in the northern hemisphere, but with a sign change. In this model, we see an excellent correspondence between 60 N and 60 S as well as between 39 N and 39 S, as guided by the black and white lines; some very small differences can be seen. The QG internal dynamics, regardless of predicted boundary layers and flux expulsion effects, is indeed visible in the magnetic data observed outside the dynamo region.
We examined the B 0 r signal in other models as well. We were able to identify some wave signals in every dynamo case, but the clarity of the signal strongly depends on the case examined. Fig. 13 compares t-/ and f-m plots of B 0 r at r ¼ r o for models 4R5, 5R5, and 6.5R2Ra at latitude 60 N. The model 4R5 for a strong field K % 18 demonstrates that the wave patterns seen in the surface field become less evident than the equivalent hu 0 s i plot of Fig. 4a. The frequency spectrum (Fig. 13b) shows some eastward moving features, which were hardly visible in Fig. 4b. This becomes more obvious in model 5R5: despite the excellent identification in hu 0 s i (Fig. 4c), it is difficult to find the corresponding patterns of the surface field. Nonetheless, the spectrum still retains the signals, although weaker, sitting around the wave dispersion relation. Model 6.5R2Ra, which demonstrated an eastward drift of hu 0 s i, illustrates magnetic eastward drifts even more clearly; the calcu- lated wave speeds (black lines) help to identify westward drifts corresponding to the internal wave motions. All this shows that detecting MR waves in the magnetic field at the top of the core will not be straightforward, compared to that in the QG flow models. Our simulations indicate that the background magnetic field for K no larger than 5 provides a reasonable observation in the surface field. It is not entirely clear yet what determines the detectability of the B 0 r signal, but it may be that it is more strongly affected by nonlinearity than the hu 0 s i signal. Nonlinear interactions between the waves and the underlying quasi-steady state may be responsible for the appearance of eastward propagating features in the frequency spectrum, but further exploration is needed.

Discussion and concluding remarks
We have presented further evidence of magnetic Rossby (MR) waves operating within rotating spherical dynamos, which are used for simulating planetary dynamos in fluid cores. The rotating MHD wave motions are non-axisymmetric but equatorially symmetric, representing a QG mode in a rotating thick shell problem. Linear theory shows that these waves will propagate retrogradely in azimuth on magnetostrophic timescales, which are given by the toroidal component of the background magnetic field with respect to the rotational rate. It therefore has the potential to infer the 'invisible' toroidal magnetic field deep down in the core.
Adopting the methodology introduced by HJT15, we performed space-time analyses of an extended range of simulation data and reported successful cases as well as a failed one. In the models explored in this study, we were able to detect MR waves if axisymmetric torsional Alfvén waves (TOs) were excited. Torsional waves are most strongly excited in dynamos with larger values of K, i.e. strong field dynamos (TJT14). We found that slow MR waves were also found at larger values of K, so that TOs and slow MR waves are seen together or not at all in the analysed simulations. As noted by Dormy (2016), strong field dynamos can be found at moderate E $ 10 À4 if Pm is large enough (e.g. run 4R5), but if E is lower, Pm does not need to be quite so big (e.g. run 6.5R2). As noted in the introduction, the existence of MR waves in magnetostrophic balance does not a priori imply the dynamo itself is magnetostrophic. Nevertheless, our numerical experience suggests that slow MR waves are seen when Dormy (2016)'s criteria for a strong field dynamo are approximately satisfied, and are not seen when they are violated. We therefore make the conjecture that the existence of slow MR waves is a signature of a strong field dynamo.
Dynamo models with strong magnetic fields are most easily found when Pm is greater than 1. Both fattened convective rolls and slower wave propagation were found even though the convection is approximately geostrophic. Generally, the form of the waves is consistent with that in linear analyses (Zhang, 1995). Pm > 1 is when the linear theory of convection predicts MR waves at onset (Hori et al., 2014). The geodynamo operates at small Pm, but at much lower E than we can reach numerically. At small Pm convective onset occurs in the form of eastward propagating diffusive modes. However, we argue that as the magnetic Reynolds number is large in the core, westward propagating non-diffusive MR waves are possibly found in the geophysical regime. The disturbances we discuss in this paper are all associated with convection, since the supercriticality has been kept close to the onset value (Ra=Ra c 16). Exploring more vigorous convective regimes would be useful, as computing resources improve. Alternative approaches such as magnetoconvection simulations and experiments King & Aurnou, 2015) may also help in this regard.
To examine the argument given in HJT15 that the waves found in the simulations are indeed MR waves, we evaluated the individual components of the z-vorticity equation. We found that the Coriolis and Lorentz terms are indeed the dominant terms, supporting the view that the waves are magnetic Rossby waves. The buoyancy term is weak in magnitude, but plays a crucial role in exciting the non-axisymmetric waves. The importance of the other terms, such as the inertia, Reynolds force, and viscosity, varies with the model parameters. They are suppressed for an Ekman number E 6 10 À5 , in the presence of a strong magnetic field with an Elsasser number K P Oð1Þ.
We also performed some simulations with the magnetic field switched off. As expected, a very different picture emerges, with eastward propagating thermal Rossby waves becoming visible. The vorticity balance is also completely changed, with Coriolis, inertia and Reynolds force now being the dominant players. We speculate that nests of convection (Brown et al., 2008), preferred longitudes of strong convective activity, may be connected with energy being transported by thermal Rossby waves into these convectively active regions.
Of geophysical importance, we examined how the waves affect the radial component of the magnetic field at the CMB, as this is what is seen in the geomagnetic secular variation. Our results showed up a possible difficulty, as although the waves can be seen in the B 0 r signal at the top of the core (see Fig. 12), when the field is very strong this signal is less clear-cut than the hu 0 s i signal. This could be because a very strong field with K > Oð1Þ tends to make the flow less geostrophic, so the signal at the CMB is not directly related to the core flow in the interior. It could also be due to the importance of nonlinear terms in the induction equation. If the perturbed field is small compared to the mean field, then we expect a simple linear relation between the perturbed field and the perturbed flow, but if the perturbation fields are comparable to the mean field the relationship is less simple. This suggests that the internal core field should be K ¼ Oð1Þ to host detectable MR wave motions; if the field is too weak no MR waves occur, if the field is too strong, nonlinearity and ageostrophy make it difficult to see evidence of the linear dispersion relation in the observed signals.
An interesting finding from this work is that nonlinearity can indeed influence the waveforms. It is known that the nonlinear dynamics of dispersive waves is distinct from that of nondispersive -sometimes called hyperbolic-waves (e.g. Whitham, 1974): dispersive modes in the inviscid, weakly nonlinear regime appear to form cnoidal or solitary waves. This may explain our observations of narrow wave crests and troughs in the low E simulations. In our simulations, however, the finite amplitude effect of the Lorentz force did not seem to impact on the wave speeds very greatly, as the linear theory gave surprisingly good results; nonlinearity has the potential to alter wavespeeds. It is also possible that nonlinearity is important in the induction equation. Since nonlinear theories on rotating MHD waves are in their infancy, this line of research could bring a new physical insight.