Magnetoconvection in a rotating spherical shell in the presence of a uniform axial magnetic field

ABSTRACT We report simulations of thermal convection and magnetic-field generation in a rapidly-rotating spherical shell, in the presence of a uniform axial magnetic field of variable strength. We consider the effect of the imposed field on the critical parameters (Rayleigh number, azimuthal wavenumber and propagation frequency) for the onset of convection, and on the relative importance of Coriolis, buoyancy and Lorentz forces in the resulting solutions. The imposed field strength must be of order one (corresponding to an Elsasser number of unity) to observe significant modifications of the flow; in this case, all the critical parameters are reduced, an effect that is more pronounced at small Ekman numbers. Beyond onset, we study the variations of the structure and properties of the magnetically-modified convective flows with increasing Rayleigh numbers. In particular, we note the weak relative kinetic helicity, the rapid breakdown of the columnarity, and the enhanced heat transport efficiency of the flows obtained for imposed field strengths of order one. Furthermore, magnetic and thermal winds drive a significant zonal flow in this case, which is not present with no imposed field or with stronger imposed fields. The mechanisms for magnetic field generation (particularly the lengthscales involved in the axisymmetric field production) vary with the strength of the imposed field, with three distinct regimes being observed for weak, order one, and stronger imposed fields. In the last two cases, the induced magnetic field reinforces the imposed field, even exceeding its strength for large Rayleigh numbers, which suggests that magnetically-modified flows might be able to produce large-scale self-sustained magnetic field. These magnetoconvection calculations are relevant to planets orbiting magnetically active hosts, and also help to elucidate the mechanisms for field generation in a strong-field regime.

influence of the Coriolis force due to the rapid rotation of the planet. Magnetic fields can also significantly affect the convective flows via a magnetic feedback force, the Lorentz force. The relative strengths of the Lorentz and Coriolis forces can be estimated with the Elsasser number (which we will define more precisely in section 2). An Elsasser number of order unity is often interpreted as an indication that both forces play an equally important role on the dynamics. In the Earth's core, the Elsasser number is estimated to be of the order of 10 (Gillet et al. 2010), so the geodynamo (and planetary dynamos in general) operate under the major influence of the Lorentz force, in addition to the Coriolis force.
Numerical simulations of convective dynamos aim to provide a better understanding of the core dynamics and the generation of planetary magnetic fields (e.g. Christensen and Wicht 2015). Although the governing equations solved by these numerical models are good approximations of the modelled system, the controlling parameters used in the models can vastly differ from realistic values. A well-know example is the problematic Ekman number Ek, which measures the ratio of the rotation period to the viscous timescale at the system size: Ek is O(10 −15 ) in the Earth's liquid core, but is typically set to O(10 −5 ) in numerical models. This artificial increase of the fluid viscosity is necessary to eliminate the small scales that cannot be resolved on the numerical grid. As a result, viscous forces play a more important role in the large-scale dynamics in numerical models than is realistic. Notwithstanding these unavoidable computational limitations, geodynamo simulations are able to reproduce some of the key characteristics of the geomagnetic field such as a mainly dipolar field. But despite this generation of a relatively strong dipolar magnetic field, the detailed role played by the magnetic field on the dynamics often remains unclear and debated (e.g. Soderlund et al. 2012, King and Buffett 2013, Oruba and Dormy 2014, Cheng and Aurnou 2016. Recent geodynamo simulations have focussed on exploring the parameter space where the sustained magnetic fields visibly affect the dynamics; these are sometimes termed "strong-field" dynamos, although this term is ambiguous. Different approaches to produce strong-field dynamos have been followed, such as using large values of the magnetic Prandtl number (the ratio of viscosity to magnetic diffusion) at the relatively high Ekman numbers considered by most studies (Dormy 2016), or performing computationally intensive simulations with Ekman numbers as low as currently feasible (Ek ≈ 10 −7 ) (Yadav et al. 2016a, Aubert et al. 2017, Schaeffer et al. 2017. Changes of the convective flow due to the Lorentz force in these dynamos include a visible increase of the typical flow lengthscale, and the loss of the regular columnar structure of the flow, which is a feature of rotationally-dominated convection.
In this work, we adopt a different approach to examine the influence of magnetic fields on convection by studying magnetoconvection, where a magnetic field is externally imposed on the convective system. The great advantage of magnetoconvection is the ability to control the strength and configuration of the imposed magnetic field precisely. The basis of our knowledge of the effect of large-scale magnetic fields on rotating convection comes from early studies of linear magnetoconvection in a planar domain (Chandrasekhar 1961); for instance, a well-known result of planar rotating magnetoconvection is that a uniform vertical magnetic field favours convection for Elsasser numbers of the order of unity, with the emergence of a preferred system-size mode. There have been numerous studies of magnetoconvection in spherical geometries using imposed azimuthal fields (e.g. Cardin and Olson 1995, Olson and Glatzmaier 1995, Gillet et al. 2007), but few using an imposed uniform axial field (Fearn 1985, Sarson et al. 1997, Sakuraba and Kono 2000, Sakuraba 2002, 2007, Gómez-Pérez and Wicht 2010, Teed et al. 2015. Here we choose the uniform axial field configuration, firstly because most planetary magnetic fields have a dominant dipolar field at the planet surface, and numerical dynamo simulations show that magnetic fields generated by rotating spherical convection are naturally organised in dipole-dominated configuration (Olson et al. 1999, Sreenivasan andJones 2011). Secondly, this configuration is particularly relevant for satellite bodies surrounded by the ambient magnetic fields of their host planets. For example, Io is thought to be able to produce its own field due to the presence of Jupiter's magnetic field providing a strong ambient field (Sarson et al. 1997). Thirdly, uniform axial fields are stable (i.e. force-free). Azimuthal magnetic fields are subject to magnetic instabilities that can onset before the convective instabilities for strong imposed field (Fearn and Weiglhofer 1991a, 1991b, 1992a, 1992b, Zhang and Fearn 1993, Proctor 1994, Zhang 1995, a situation that we wish to avoid for simplicity. Linear magnetoconvection in a rotating spherical shell with an imposed uniform axial field was investigated in a numerical study by Sakuraba (2002). He showed the existence of an overall minimum of the critical Rayleigh number (which measures the strength of the buoyancy driving against diffusive effects) for Elsasser number of the order unity, similarly to the case of magnetoconvection in the presence of a uniform toroidal magnetic field (Fearn 1979). Sakuraba found that a variety of convective modes can be unstable, depending on the strength of the imposed magnetic field. As the Elsasser number is increased, geostrophic, magneto-geostrophic and magnetostrophic modes are successively preferred at the onset of convection, where these modes are distinguished by the increasing role played by the Lorentz force in the primary force balance. In the nonlinear model of Sarson et al. (1997Sarson et al. ( , 1999, which uses a two azimuthal mode approximation, a weak-field solution and a strong-field solution were observed, depending on the strength of the buoyancy forcing. In the latter case, the system essentially behaves like a self-sustained dynamo without being strongly affected by the imposed field. In a study of the fully nonlinear system, Sakuraba and Kono (2000) found that the transition between the two types of solutions occurs abruptly for an Elsasser number of order unity. The transition is characterised by an increase in the kinetic and induced magnetic energies, while the convective cells retain the same azimuthal lengthscale. On the other hand, Gómez-Pérez and Wicht (2010) found that the flow lengthscale visibly increases in the presence of a strong imposed field. Interestingly, Sarson et al. (1997Sarson et al. ( , 1999 and Sakuraba and Kono (2000) found evidence of hysteretic behaviour in this system, where the final state depends on the initial condition. The present paper follows on from these earlier studies and aims to investigate the modification of the convective flows in the presence of a strong magnetic field and, in particular, the morphology of the magnetic fields that these modified flows induce.
Dynamos driven by planar rotating convection can produce strong large-scale magnetic fields that disturb the small-scale convective flows and lead to the emergence of box-size flows (Stellmach andHansen 2004, Cooper et al. 2020). These box-size flows are not able to sustain large-scale magnetic fields themselves, leading to intermittent dynamo behaviour, where the system fluctuates between strong-field and weak-field states. This intermittent behaviour has not been reported so far in spherical dynamo simulations with strong fields (e.g. Dormy 2016, Schaeffer et al. 2017. By studying whether magnetically-modified flows can induce large-scale magnetic fields, we can speculate whether an intermittent behaviour might be expected in spherical dynamos. Using spherical magnetoconvection results, Gubbins (2000a, 2000b) previously suggested that this behaviour might indeed be relevant for the geodynamo.
The linear magnetoconvection study of Sakuraba (2002) shows that the relatively large Ekman number adopted in previous nonlinear studies (Ek ≈ 10 −4 ) misses some of the complexity of the linear convective modes observed at smaller Ek (Ek ≈ 10 −5 ). Furthermore, fully nonlinear studies were carried out in the moderately supercritical regime (with Rayleigh numbers less than twice the critical value at onset). We will therefore model magnetoconvection for values of Ek that allow us both to study the rich small-Ek dynamics and to explore widely the supercritical regime (with Rayleigh numbers up to 50 times the critical value).
The layout of the paper is as follows. The mathematical and numerical formulation of the model is described in section 2. In section 3, we present the effects of the imposed magnetic field on the onset of convection. We discuss the evolution of the flow properties in section 4, including changes in the efficiency of the convective heat transfer and the formation of zonal flows. The characteristics of the induced magnetic field and the generation of the axisymmetric field are presented in section 5. Finally, we discuss the significance of our results in section 6.

Governing equations
We study thermal Boussinesq convection in a rotating spherical shell in the presence of an externally imposed magnetic field. The convection is driven by imposing a temperature difference T between the inner sphere at radius r i and the outer sphere at radius r o . The aspect ratio is set to χ = r i /r o = 0.35. The system rotates at the rotation rate about the z-axis. Gravity varies linearly in radius, g = −g o r/r o e r . The imposed magnetic field is uniform and axial, B 0 = B 0 e z . This is implemented via a boundary condition, as described in section 2.2. We use D = r o − r i as the unit of length. The computational domain is therefore located between r i /D ≈ 0.54 and r o /D ≈ 1.54. Times are scaled by D 2 /ν, the velocity u by ν/D, the temperature T by T, the magnetic field B by (ρμηΩ) 1/2 , and the pressure P by ρνΩ, where ρ is the density, ν the kinematic viscosity, η the magnetic diffusivity, and μ the magnetic permeability. The dimensionless governing equations are where B is the total magnetic field comprised of both the imposed and induced parts, and j = ∇ × B is the electric current density. The dimensionless parameters are the Ekman the Rayleigh number the Prandtl number and the magnetic Prandtl number where κ is the thermal diffusivity and α the thermal expansion coefficient. In all our simulations, Pr and Pm are kept fixed to 1. Note that we use the definition of the Rayleigh number from the dynamo benchmark (Christensen et al. 2001), which is related to the traditional definition of the Rayleigh number, R = αg o TD 3 /(νκ), by Ra = R Ek/Pr. The boundary conditions are fixed temperature, electrically insulating, no-slip and impenetrable at r = r i and r = r o .

Numerical method
All our simulations are performed using the code PARODY, which was originally written by Dormy et al. (1998), and subsequently modified, parallelised and optimised by J. Aubert andE. Dormy (e.g. Aubert et al. 2008, Dormy 2016). PARODY is a Fortran pseudo-spectral code that solves the magnetohydrodynamic equations (1)-(5) for a Boussinesq fluid in a 3D spherical geometry. The code was previously benchmarked against five independent numerical codes used in the geodynamo community (Christensen et al. 2001). The velocity and magnetic fields are decomposed into poloidal and toroidal scalars, which are then expanded in spherical harmonics Y m l in the angular coordinates, with l the degree and m the order. A second-order finite difference scheme is used on an irregular radial grid, which is finer near the boundaries and uses a geometrical progression for the radial increment. A Crank-Nicolson scheme is implemented for the time integration of the diffusion terms and an Adams-Bashforth procedure is used for the other terms. The poloidal-toroidal decomposition and the spherical geometry allow a relatively simple implementation of the magnetic boundary conditions. We use the SHTns library for the spherical harmonic transforms (Schaeffer 2013). The numerical resolution used for each simulation is given in table C1 in appendix C.
We have modified the code to include an externally imposed magnetic field. The imposed field is axial and uniform, B 0 = B 0 e z , and so, it only has a poloidal component that projects on the spherical harmonic (l, m) = (1, 0). The field is imposed by modifying the boundary condition at r = r o for this specific harmonic of the poloidal magnetic field, giving dp 0 where p 0 1 (r) gives the spectral coefficient of the poloidal magnetic field of degree l = 1 and order m = 0 at radius r. The derivation of this boundary condition is given in appendix A.
For the initial condition, a small perturbation is added to all the spectral coefficients of the temperature field. All the other fields are set to zero, except in the case of no imposed magnetic field (B 0 = 0), where the magnetic field is initialised by an axial dipole and a toroidal field of harmonic (l, m) = (2, 0), both of amplitude of order unity in the dimensionless unit.

Definitions of the output quantities
The Elsasser number is often used in dynamo studies to estimate the ratio of the Lorentz force to the Coriolis force. It is defined as where B is a r.m.s. measure of the magnetic intensity in dimensional form. In our units, Λ = B 2 * , where the dimensionless r.m.s. magnetic intensity is B * = √ 2E m t , with E m the volumetric (dimensionless) magnetic energy density in the fluid shell. To quantify the strength of the imposed field, we will quote values of B 0 . Previous studies (e.g. Sakuraba 2002) use instead the Elsasser number based on B 0 , which is simply Λ 0 = B 2 0 . To measure the ratio of the kinetic to magnetic energies, we use the squared Alfvén number (Aubert et al. 2017) where U is the r.m.s. velocity in dimensional form. (This can also be identified as the square of the ratio of the fluid velocity to the Alfvén velocity.) In our units, A 2 = PmEkU 2 * /B 2 * , where the dimensionless r.m.s. velocity is U * = √ 2E k t , with E k the volumetric kinetic energy.
We define the Reynolds number of the convective flow as where E p k is the volumetric (dimensionless) kinetic energy density based on the poloidal velocity, and · t is a time average taken over at least one third of a viscous timescale.
The magnetic Reynolds number is defined as Rm = DU/η = PmU * . The efficiency of the convective heat transfer is usually measured by the Nusselt number Nu, which is the ratio of the total heat flux to the conducted heat flux in the absence of convection. In our system, Nu is calculated at the outer boundary and is where · S denotes an average over a spherical surface, χ = r i /r o is the conducted heat flux at r o , and T is the total temperature, which includes the static background temperature. (The conducted heat flux is simply the negative of the radial derivative of the static background temperature.) Selected output values from our simulations are given in table C1 in appendix C. A complete list and the numerical resolution is given in Data.xlsx to be found in the supplementary material.

Critical parameters
We first examine the effect of the imposed axial field on the three critical parameters at the onset of convection: the critical Rayleigh number Ra c , the marginally stable azimuthal wavenumber m c , and the frequency of the marginally stable mode ω c . The onset of convection is determined numerically using PARODY by observing the growth rate of the kinetic energy of each azimuthal wavenumber at given Ra. The critical parameters are given in tables 1-2 for B 0 ∈ [0, 10] and Ek ∈ {10 −4 , 10 −5 }. In table 1 for B 0 = 0, we compare Ra c (which we denote Ra 0 c hereafter) and m c with the theoretical values obtained from the asymptotic analysis of Dormy et al. (2004) in the case of non-magnetic rotating convection Note: The last two columns for B 0 = 0 indicate the theoretical values Ra a and m a obtained in the asymptotic study of non-magnetic rotating spherical convection of Dormy et al. (2004). * The onset of the modes m = 2 and 3 is very close for B 0 = 3 and Ek = 10 −5 , both having a critical Rayleigh number of 64.0.  with differential heating in a spherical shell with χ = 0.35 and no-slip boundaries. In the asymptotic theory, Ra c and m c both scale as Ek −1/3 in the leading order asymptotics. Note that we have used the relevant scaling factor to account for the different definitions of the Ekman and Rayleigh numbers (recall in particular that Ra = REk/Pr in our definition). The agreement between numerical and theoretical values is good and reduces when Ek decreases, as expected. Figure 1(a) shows Ra c as a function of B 0 . We find that B 0 must be greater than 1 to significantly influence the onset of convection for both Ekman numbers. As B 0 increases above unity, Ra c first decreases and reaches a minimum at B 0 = 3 for Ek = 10 −4 and at B 0 = 4 for Ek = 10 −5 . The decrease in Ra c is significant, reaching a minimum of 0.65Ra 0 c for Ek = 10 −4 and 0.56Ra 0 c for Ek = 10 −5 . The decrease of Ra c is accompanied by a decrease of m c , which is most noticeable at the smallest Ek, as seen in figure 1(b). At larger B 0 , Ra c increases, reaching a value that is approximately independent of Ek and larger than Ra 0 c for both Ek. This increase is accompanied by an increase in m c , which is again most noticeable at the smallest Ek.
The evolution of Ra c and m c with B 0 is similar to that observed in the linear magnetoconvection study of Sakuraba (2002), which mainly focussed on the case of internal heating with no inner core and zero flux boundary condition at r = r o . The main difference between their results and ours is that Sakuraba observed that the preferred mode at the onset is an axisymmetric mode for B 0 ∈ [1, 2]. This axisymmetric "polar" mode consists of a large-scale meridional circulation that crosses the equatorial plane and transports heat between the two hemispheres. The absence of this polar mode at the onset in our simulations is presumably due to our choice of inner core size and fixed temperature boundary conditions, as Sakuraba shows that the stability curve of the polar mode increases towards larger Ra in these conditions.
The qualitative behaviour of Ra c is also similar to that expected from linear magnetoconvection in rotating planar geometry (Chandrasekhar 1961). In the case where rotation vector and imposed magnetic field are vertical, Ra c starts to decrease when B 0 = O(Ek 1/6 ), while the lengthscale of the convection switches to a system-size scale (Roberts and King 2013). This change corresponds to a change in the primary force balance from a viscously-dominated regime for B 0 < O(Ek 1/6 ), where the viscous force is needed to break the Proudman-Taylor constraint, to a magnetically-dominated regime, where the Lorentz force becomes sufficiently large to assume this role. For the modest values of Ek considered here, the regime change is expected to occur around B 0 ≈ Ek 1/6 ≈ 0.1. In our configuration, we find that the regime change occurs for larger B 0 , B 0 = O(1), with no apparent dependence on Ek. However, our simulations are relatively coarse-grained in B 0 , while Ek only varies by a decade. We would therefore not be able to pick up an Ek 1/6 trend. In planar magnetoconvection, the minimum of Ra c occurs for B 0 = O(1) and, in the limit B 0 → ∞, Ra c scales as B 2 0 and the critical horizontal wavenumber as B 0 (Roberts and King 2013). In our spherical system, we also find that Ra c becomes independent of Ek at large B 0 , although we have too few data points to check the scaling with B 0 . However, the critical wavenumber retains an Ek-dependence with no clear dependence on B 0 . Table 2 gives the frequency ω of the waves at the dominant azimuthal wavenumber m n for simulations with Ra just above Ra c . By dominant wavenumber, we mean the wavenumber corresponding to the peak in the power spectrum of the kinetic energy. The frequencies were obtained by calculating the azimuthal drift of the dominant mode in Hovmöller maps (longitude-time maps at a fixed radius in the equatorial plane) of the radial velocity. Note that any azimuthal drift due to the zonal velocity was subtracted to extract the frequency of the waves. The frequency can be compared with the expected frequency of thermal Rossby waves of azimuthal wavenumber m n at the onset of non-magnetic convection, which we estimate by using the theoretical values provided by the asymptotic study of Dormy et al. (2004) for our configuration (see appendix B). Outside the tangent cylinder (where the height of the outer sphere decreases outwards), thermal Rossby waves always propagate in the prograde direction and their frequency decreases for increasing azimuthal wavenumber. For B 0 = 0, ω is close to the theoretical frequency of the thermal Rossby waves. For B 0 > 0, ω is much slower than the frequency of the thermal Rossby waves, indicating that the Lorentz force plays a significant role in the propagation mechanism. Near the minimum of Ra c (B 0 = 3 and 4 for Ek = 10 −4 and 10 −5 respectively), the increase of the frequency of the dominant mode m n = 4 when Ek decreases is relatively weak compared with the increase expected for the thermal Rossby wave (scaling as Ek −2/3 in the leading order asymptotics, see appendix B). For B 0 = 10 at Ek = 10 −5 , the wave propagates in the retrograde direction, revealing the dominance of the Lorentz force (Finlay 2008). The force balance in these waves will be analysed in section 3.3 below. Figure 2 shows the 3D isosurfaces of the radial velocity u r for selected B 0 with Ra just above the onset of convection for Ek = 10 −4 and 10 −5 . For each Ek, the three cases were chosen to represent the viscously-dominated regime (B 0 = 0) and the field-dominated regimes at the minimum of the Ra c stability curves (B 0 = 3 or 4) and at large B 0 (B 0 = 10). In all cases, the convective flow is mainly columnar. For B 0 = 0, we observe the well-known structure of thermal Rossby waves with a tilt of the column in the prograde direction (Zhang 1992). For B 0 ∈ [3, 4], the flow adopts a banana shape with a visible z-dependence, especially for Ek = 10 −5 . For B 0 = 10, the columns are straight with no tilt. The preservation of the columnar shape of the convection cells for B 0 > 0 can explain why the axial magnetic field starts to influence convection for B 0 = O(1) in our simulations, while studies of spherical magnetoconvection with an imposed azimuthal field find that a smaller field amplitude  (Fearn 1979). In the axial field problem, the magnetic perturbations are created by the distortions of the imposed field due to the z-gradients of the velocity (via the (B 0 · ∇)u part of the induction term), as opposed to the azimuthal field problem, where they are mainly due to the φ-gradients. For columnar flows, the z-gradients of the velocity are much smaller than the φ-gradients, so the effect is only comparable for larger B 0 .

3D structure of the flow
To study the axial structure of the flow, we consider the modified Proudman-Taylor constraint following Sakuraba (2002). Taking the curl of the Navier-Stokes equation (1) and assuming that inertia, buoyancy and viscous forces are negligible gives where the nonlinear contributions to the Lorentz force have also been neglected. Since the imposed field is uniform, j is only produced by the induced magnetic field. We define the modified velocity as so that the modified Proudman-Taylor constraint is ∂ u/∂z ≈ 0. Figure 3 shows the s-components of u, B 0 /(2Pm)j and u on the axial cylindrical surface of radius s = 0.7 for the three cases of figure 2 at Ek = 10 −5 . The onset of the dynamo for B 0 = 0 occurs at larger Ra (Ra ≈ 500) so this case does not have a magnetic field. For B 0 = 4, the contours of u s are bent along the z-axis, which gives the banana shape observed in 3D in figure 2. In this case, B 0 /(2Pm)j s opposes u s near the equatorial region, which produces a modified velocity that is fairly z-invariant. The system therefore closely follows the modified Proudman-Taylor constraint (15). For B 0 = 10, both u s and u s are largely z-invariant, although the deviation from z-invariance is greater near the outer boundary for the velocity. In this case, B 0 /(2Pm)j s is greater near the outer boundary, where it is correlated with u s . Interestingly, the axial structure of u s for B 0 = 4 and 10 is similar to that of u s for B 0 = 0, with the maxima of the absolute value located near the outer boundary. This is perhaps logical, given that the increased significance of viscosity in the boundary layers acts to reduce the applicability of the (modified) Proudman-Taylor constraint there, allowing for local variations. Moreover j s has a contribution from ∂B φ /∂z, and since the magnetic boundary condition requires the toroidal scalar to go to zero at the outer radius, the maxima of j s at the boundary is also unsurprising. The fact that the maxima of |u s | are encountered around the equatorial region rather than near the outer boundaries (for B 0 > 0) has an effect on the convective heat transport, as shown in section 4.2.

Force balance
To conclude this section, we investigate the dynamical force balance near the onset of convection. To do so, we look at the vorticity equation to eliminate the pressure term. Since the flow is mainly columnar, we consider the equation of the z-component of the vorticity ω z , where Figure 4 shows ω z and the terms in equations (18)-(22) on the axial cylindrical surface of radius s = 0.7 for the three cases of figure 2 at Ek = 10 −5 . The terms are plotted using individual snapshots that are representative of the flow over time. For B 0 = 0, the Coriolis and buoyancy terms are the most significant in the force balance; the viscous term is non-negligible, but is subdominant. Positive buoyancy terms are largely balanced by negative Coriolis terms, both being largest between anticyclones and cyclones (in the prograde direction); with opposite signs, there is a similar balance between cyclones and anticyclones. For B 0 = 4, the Lorentz term is one of the dominant forces, alongside the Coriolis and buoyancy terms (a MAC balance). The Coriolis and Lorentz terms balance each other at higher latitudes, whereas the contribution of the buoyancy term is more relevant near the equatorial region. Here, "sinks" in the buoyancy term are opposed by the Lorentz force, whilst buoyancy "sources" are balanced by the Coriolis force. For B 0 = 10, the force balance is primarily between the buoyancy term and the Lorentz term in the equatorial region (a MA balance), but towards the outer surface the balance is again between Coriolis and Lorentz contributions. The sum of the r.h.s. terms is relatively large for B 0 = 0, consistent with the greater frequency of propagation seen for that case in table 2. In other words, one effect of the Lorentz force for B 0 > 0 is to reduce this prograde propagation, as noted earlier for the frequencies at onset.

Evolution of the flow with increasing Rayleigh number
In the rest of paper, we will present results obtained at Ek = 10 −5 for varying Ra. Similar observations can be drawn from the simulations at Ek = 10 −4 ; exceptions are discussed in the text. In this section, we study how the presence of the ambient magnetic field affects the different components of the flow. Figure 5 shows the evolution of the Reynolds number based on the poloidal velocity Re p (defined in equation (13)) as a function of Ra for B 0 = 0, 4 and 10. Cases with B 0 = 4 produce stronger Re p close to onset than cases with B 0 = 0 and 10. However, for Ra away for the onset, there are no particularly notable differences in Re p for different B 0 , although B 0 = 4 tends to produce slightly smaller The magnetic Prandtl number is set to 1 in all our simulations, so the magnetic Reynolds number Rm is equal to the Reynolds number. In self-sustained spherical convective dynamos, the magnetic Reynolds numbers (based on the total velocity) needs to be greater than approximately 40 for dynamo action (Christensen and Aubert 2006). This is verified in our simulations with B 0 = 0, where the dynamo onsets at Ra = 500 with a Reynolds number based on the total velocity of 73 and Re p = 32. Since the Reynolds number is essentially similar for all B 0 at a given Ra, we might consider that cases with Ra > 500 for B 0 > 0 have magnetic Reynolds numbers that are large enough for dynamo action, so these solutions might behave essentially as dynamo solutions. This point will be explored in more detail in section 5. To study the evolution of the columnar structure of the flow with Ra, we estimate the variations of the axial vorticity along z using the measure of the columnarity C ωz proposed by Soderlund et al. (2012). C ωz is based on the non-axisymmetric vorticity ω and defined as

Convective flow
where · z denotes the average in the axial direction. The integral is evaluated as a sum on a discrete finite grid, where the summation is made over cylindrical radius s and φ outside the tangent cylinder. This measure is calculated from taking the time average of C ωz from various data snapshots. Figure 6(a) shows C ωz as a function of Ra for the same parameters as figure 5. The values of C ωz near the onset of convection are consistent with our observations in section 3, i.e. the flow at B 0 = 4 (B 0 = 10) is less (more) z-invariant than for B 0 = 0. Soderlund et al. considered that the columnar structure breaks down when C ωz < 0.5. By this definition, the columnar flows first breaks down for B 0 = 4 for Ra ≈ 350. The degradation of the columnarity occurs later for B 0 = 0 and 10. Examples of the 3D structure of the flow are shown in figure 7 for Ra = 350. As expected, the columnarity has visibly broken down for B 0 = 4, but the flow remains columnar to a good degree for B 0 = 0 and B 0 = 10. In this latter case, the flow displays quasi-2D sheetlike structures. Convection in this case is heterogeneous with quiescent regions located between the pair of intense radial jets. This sheet-like structure is still noticeable for higher values of Ra but becomes less uniform between Ra = 500 and Ra = 750.
To further characterise the flow structure, we measure the kinetic helicity, which describes the spatial correlation between the components of the velocity and vorticity and is defined as The helicity is often thought to be essential to the generation of large-scale magnetic fields (Olson et al. 1999, Sreenivasan and Jones 2011, Moffatt and Dormy 2019 and is therefore an interesting measure to assess the dynamo capability of the flow. Since we study the evolution of the flow with varying Ra, we use the relative helicity defined as (Olson et al. 1999) where · h corresponds to the average in the northern or southern hemispheres. Similarly to C ωz , the measure is calculated from taking the time average from different data snapshots. H rel is negative (positive) in the northern (southern) hemispheres (Olson et al. 1999), so we average the absolute value of the two hemispheres. Figure 6(b) shows H rel as a function of Ra. For B 0 = 0, the relative helicity is approximately constant for nonmagnetic solutions (Ra < 500). There is a small increase in the helicity after the onset of dynamo action, while the saturated field remains weak (for 500 ≤ Ra ≤ 750, where Λ < 1; see section 5); but for higher Ra the helicity decreases significantly. After the columnarity breaks down, the time-averaged helicity in each of the two hemispheres remains approximately equal in absolute value. Around 50% of the relative helicity is produced by the correlation of the z-components, u z ω z h . This percentage decreases slightly at the highest values of Ra, coinciding with the decrease in columnarity. Although the correlation of the z-components is always the largest contributor to the relative helicity, the dominance of the z-components is not as great as might naively be expected. In their asymptotic (small Ek) linear magnetoconvection analysis (with an imposed azimuthal field), Sreenivasan and Jones (2011) note that the contributions from u z ω z and u s ω s are equal. Nevertheless some authors (e.g. Soderlund et al. 2012) concentrate on the "axial helicity", u z ω z ; and indeed, the variations in this quantity are informative. In the presence of an imposed axial field (for B 0 = 4 and 10), the correlation between u z and ω z is poor, and the values of the relative helicity are always small. One possible consequence is that these flows would be inefficient for the generation of large-scale magnetic fields (in the absence of the imposed field). This is in contrast to the conclusions of Sreenivasan and Jones (2011), who found enhanced correlation of u z and ω z (and thus enhanced helicity), from calculations with no imposed field, and also from linear magnetoconvection analysis with an imposed azimuthal field. (Most of their calculations used stress-free boundary conditions, where this effect is stronger; but the effect is also present with no-slip conditions, as in the calculations presented here.) They note the effect as stronger for solutions of dipolar symmetry (i.e. B r antisymmetric with respect to reflection in the equator), and argue that this effect may explain the dominance of this symmetry for planetary dynamos. Our imposed vertical field has this equatorial symmetry, but the Lorentz force produced by magnetoconvection in its presence clearly does not enhance the vertical flow within convective columns (and hence the helicity) in the same way. In the study of Sreenivasan and Jones, the Lorentz force was dominated by B φ , which is zero near the equator so cannot balance the buoyancy force there. The Coriolis force must therefore balance the buoyancy force near the equator, meaning that the z derivative of u z cannot be zero there, giving rise to coherent axial flows along the columns. In our case the Lorentz force is dominated by B z , so the buoyancy and Lorentz forces are able to balance near the equator. Figure 8 shows the spectra of the kinetic energy with respect to the spherical harmonic degree and order for selected Rayleigh numbers. For B 0 = 0 and 10, the peak of the spectra in m always remains close to their respective critical mode at the onset of convection. For B 0 = 0, there is no drastic change in the spectra for the non-magnetic (Ra = 350) and dynamo (Ra = 750) cases, except for a slight shift of the dominant m towards smaller values. For B 0 = 4, the peak is broader, even close to onset, and the energy is increasingly transferred to smaller azimuthal wavenumbers for increasing Ra.
Finally, we conclude this section by looking at the activity inside the tangent cylinder. Convection inside the tangent cylinder onsets when 750 < Ra < 1200 for B 0 = 0. The convection begins at the edge of the tangent cylinder, encroaching from the convection flows in the rest of the shell, before spreading to all latitudes as Ra increases. The delay of the convection onset inside the tangent cylinder is well-documented in studies of nonmagnetic convection (Jones 2015) and is caused by the hindering effect of rotation on the linear onset, which is more pronounced towards the poles, where the direction of gravity and the rotation axis are aligned. For B 0 = 4, the convection onsets inside the tangent cylinder earlier than for B 0 = 0 (when 350 < Ra < 500); in this case, the presence of both rotation and an imposed field with Elsasser number of order unity is favourable to convection, similarly to the situation outside the tangent cylinder. For B 0 = 10, convection has not occurred inside the tangent cylinder for the highest values of Ra that we have considered (Ra = 2800). This delay is also expected from results of linear magnetoconvection in plane layer geometry with no rotation or high Elsasser number (Chandrasekhar 1961): analogously to the rotating case, the hindering effect of the magnetic field on the convection onset is more pronounced towards the poles when the direction of the imposed field is aligned with the direction of gravity.

Efficiency of the heat transfer
We now study whether the changes to the 3D structure of the radial velocity for varying B 0 affect the efficiency of the heat transfer. We calculate the convective heat flux as where the overbar denotes an azimuthal average. Contour plots of F are shown in figure 9 for Ra near the onset of convection. Note that these cases are not exactly located at the same Ra/Ra c , so we compare the shape of the contours rather than the amplitude of F, and check that our observations remain valid for varying Ra around these values. F occupies  To measure the improved convective heat transfer in the presence of the imposed field, we look at the evolution of the Nusselt number (defined in equation (14)) with Ra in figure 10. Nu − 1 is approximately 2 orders of magnitude larger near the onset for B 0 = 4 than for B 0 = 0 and 10, so the efficiency of the convection is greatly improved near onset. As Ra increases to values greater than 1000, the three cases tend to similar values of Nu − 1. (Note that this is despite the slightly smaller Re p values noted for B 0 = 4 at high Ra; i.e. notwithstanding the relatively lower Re p , the B 0 = 4 flow remains efficient at transferring heat.) This is somewhat surprising given the differences observed in the kinetic energy spectra at larger Rayleigh numbers, and in particular, that the flow is dominated by different lengthscales (this remains true at our largest Rayleigh number). This shows that global diagnostics such as Re p and Nu do not fully capture the significant differences observed in the flow at large Ra. In all cases, the measure of columnarity indicates that the columnar structure has nearly broken down for these Rayleigh numbers for all three B 0 , which could allow equally vigorous convection in all 3 cases. There has been considerable discussion of the variation of the heat flow in simulations (and experiments), and on its scaling with varying Ra, Pr and Ek (see Jones 2015, for a review and for references).

Zonal flow
Dynamo simulations tend to have relatively weak zonal (i.e. axisymmetric azimuthal) flows compared with non-magnetic convection simulations (Aubert 2005), unless heterogeneous boundary conditions are used (Aubert et al. 2013). This is caused by the action of the Lorentz force to prevent the shearing of the poloidal field. We assess how the presence of an imposed axial magnetic field affects this behaviour. Figure 11 shows the ratio of the kinetic energy of the zonal flows to the total kinetic energy as a function of Ra. For B 0 = 0, the zonal flow amplitude is visibly reduced when the dynamo starts at Ra = 500, as expected from previous dynamo studies. This ratio is much stronger for B 0 = 4 than for the other cases, reaching up to 30% for Ra = 200. The strong axial field for B 0 = 10 largely suppresses zonal flows, whose energy remains around 1% only of the total kinetic energy. Figure 12 shows a meridional slice of the zonal flow for the case B 0 = 4 and Ra = 200. The zonal flow is retrograde near the inner core and prograde near the outer boundary, which is a similar to the pattern observed in non-magnetic rotating spherical convection at moderate Ra (Christensen 2002). The zonal flow has a strong degree of axial dependence in the region of the convective columns, with the largest amplitude near the equatorial region. Ferraro's law of isorotation (Ferraro 1937) states that, in the limit of high electrical conductivity (high Rm), the field lines of the time-averaged axisymmetric poloidal field lines and the isocontours of the time-averaged angular velocity, u φ /s, adjust to become aligned, thereby minimising the production of an azimuthal magnetic field via an omega effect (see e.g. Roberts 2015). To check the validity of Ferraro's law in our simulations, we superpose the axisymmetric poloidal magnetic field lines in figure 12. While the poloidal field remains largely directed along z, the z-dependence of the zonal flow means that the isocontours of the zonal flow cross the poloidal field lines around the mid-latitudes. This is at odds with Ferraro's law, but is not unexpected for a fluid of finite electrical conductivity.
The presence of strong zonal flows for B 0 = 4 is surprising because the convective columns near the onset of convection do not have a tilt as for B 0 = 0. In the non-magnetic case, the column tilt produces the correlation between radial and azimuthal velocities, leading to significant Reynolds stresses that generate the zonal flows (Busse and Hood 1982). We therefore anticipate that the Reynolds stresses are not the primary source for the zonal Figure 11. Ratio of the kinetic energy of the zonal flows to the total kinetic energy as a function of Ra for Ek = 10 −5 . (Colour online) flows for B 0 = 4. To determine the sources and sinks of the zonal kinetic energy for B 0 = 4, we examine the power budget for the zonal flow, following Aubert (2005). To do so, we multiply the azimuthal average of the φ-component of the Navier-Stokes equation (1) by u φ to obtain 1 2 where Figure 13 shows the time average of each of these terms in a meridional slice for the case B 0 = 4 and Ra = 200. The power budget is dominated by the contributions from the Coriolis and Lorentz terms, which nearly balance each other, with the Lorentz term acting mainly as a sink near the equatorial region and a source towards mid-latitude. The contribution from the Reynolds stresses is small and acts as a sink around the tangent cylinder and a source at larger cylindrical radius, producing a mostly invariant pattern along z (and opposing the residual of P L + P C ). The viscous term is a sink term in most of the domain,  the z-gradient of u φ shown in plot (a). The two plots are nearly identical, which indicates that the viscous and nonlinear inertial terms only play a minor role in the balance. The thermal wind is stronger near the equatorial plane, while the magnetic wind is stronger close to the outer boundary. The magnetic wind is mainly due to the interaction of the magnetic perturbations with the imposed field and opposes the shearing of the poloidal magnetic field lines by the z-gradients of the zonal flow. The thermal wind is strong in this case because the convection is efficient at transporting heat as described in section 4.2, which creates strong latitudinal variations of the temperature. The enhanced zonal flows observed for B 0 = 4 (shown in figure 12) are therefore indirectly produced by the enhanced convective transport. Sakuraba (2007) studied a case of magnetoconvection with imposed axial field for Ek ≈ 10 −5 , Ra ≈ 500 and B 0 ≈ 1.4 (where we have rescaled their parameters to match our definitions). They also observe the formation of thermal-wind driven zonal flows that are retrograde near the inner core. This is despite the different choice of thermal boundary conditions, with Sakuraba using mixed boundaries (fixed temperature at the inner core r = r i and zero flux at r = r o ). The shape of the isotherms are very sensitive to the choice of thermal boundary conditions, but the same choice of inner boundary condition seems to result in similar behaviour here, regardless of the differing outer condition.

Evolution of the magnetic energy with Ra
In this section, we study the characteristics of the magnetic field induced by the flow for varying Ra and B 0 . We first study the evolution of the strength of the magnetic field. Figure 15(a) shows the Elsasser number Λ (see definition in section 2.3) as a function of Ra. For B 0 = 0, the Elsasser number increases monotonically up to values greater than 1 for Ra > 1000. For Ek = 10 −5 , the magnetic field always remains dominated by the dipole in the range of Rayleigh numbers investigated here, which goes up to Ra = 2800. Higher values of Ra are required to enter the multipolar regime (Christensen and Aubert 2006). For Ek = 10 −4 , the multipolar regime is entered at smaller values of Ra (Ra = 2100), at which point the Elsasser number decreases below unity. When B 0 > 0, Λ drops towards Λ 0 (= B 2 0 ) near the onset of convection, so the solution branch is continuous at Ra = Ra c . For B 0 = 4, reaches a local maximum for Ra = 100, before increasing monotonically for Ra > 150. This local maximum corresponds to the kink in the plot of Re p observed in section 4.1, and is the result of a change of branch between Ra = 100 and 150. For B 0 = 10, increases steadily after the onset of convection. Figure 15(b) shows the ratio of Λ/Λ 0 as a function of Ra. The induced field exceeds the imposed field (Λ/Λ 0 > 2) for B 0 = 4 at Ra > 750. For B 0 = 10, this situation has not yet occurred, although the amplitude of the induced field gets increasingly close to the imposed field: Λ/Λ 0 ≈ 1.75 for our largest Rayleigh number. For both B 0 = 4 and 10, there are no obvious changes in the evolution of Λ around Ra = 500, where the magnetic Reynolds number reaches the value of the dynamo threshold in the absence of imposed field. Figure 15(c) shows the evolution of the squared Alfvén number, A 2 , the ratio of the kinetic to magnetic energies as defined in section 2.3 (equation (12)). For reference, A 2 is thought to be of the order of 10 −4 in the Earth's core (Aubert et al. 2017). For B 0 = 0, A 2 decreases with increasing Ra, but always remains fairly large, never falling below 0.1. For B 0 = 4 and 10, A 2 is much smaller, of the order of 10 −5 and 10 −7 close to the onset respectively. However, in both cases, A 2 increases with Ra up to approximately 10 −2 and 5 × 10 −3 . The kinetic energy therefore increases more rapidly than the magnetic energy in the magnetoconvection case, contrary to the dynamo case.
Note that we did not observe notably large temporal fluctuations of the kinetic and magnetic energies in any of our simulations. For instance, at Ra/Ra c ≈ 12, the kinetic energy fluctuates by approximately 10% around the mean value, while the fluctuations of the magnetic energy tend to be smaller at larger B 0 (around 10% of the mean value for B 0 = 0 and 1% for B 0 = 10). Figure 16(a) shows the evolution of the toroidal ratio i.e. the ratio of the toroidal magnetic energy to the total magnetic energy. For B 0 = 0, the toroidal ratio is approximately 0.5 for all Ra. For B 0 > 0, the toroidal ratio is much smaller: simulations at B 0 = 4 slowly increases above 0.01 up to around 0.3, while for B 0 = 10, the toroidal ratio is less than 10 −7 near the onset of convection and increases up to 10 −2 for the largest Ra. The poloidal component therefore largely dominates the magnetic field in magnetoconvection with an imposed axial field, while the magnetic energy is equally distributed between poloidal and toroidal components in our dynamo simulations. This observation is consistent with the results of Sakuraba and Kono (2000). Figure 16(b) shows the evolution of the ratio of the axisymmetric toroidal to the toroidal magnetic energies. For B 0 = 0, this ratio remains constant around 0.15 for the highest Ra. For B 0 = 4, approximately 50% of the toroidal magnetic energy is contained in the axisymmetric mode around Ra = 200. This observation is consistent with the formation of strong zonal flows in this case, as the zonal flows distorts the poloidal magnetic field to produce the axisymmetric toroidal field via an omega effect. The zonal flows are relatively weaker at larger Ra, so the axisymmetric toroidal field weakens compared to the total toroidal field then. For B 0 = 10, the axisymmetric toroidal field is weak, especially at small Ra, which is consistent with the weak zonal flows observed in this case. To conclude this section, we study the partition of the poloidal magnetic energy between dipolar and multipolar components. Figure 17 shows the spectra of the time-averaged poloidal magnetic energy for fixed Rayleigh numbers. The dominant degree is l = 1 even for the highest values of Ra, indicating the dipolar nature of the field. As mentioned earlier, for Ek = 10 −4 at B 0 = 0, we found that the field becomes multipolar for Ra = 2100, but no multipolar regime was found for Ek = 10 −5 for Rayleigh number up to Ra = 2800. Christensen and Aubert (2006) found that convective dynamos produce multipolar fields when the local Rossby number (Ro l = U/(ΩL u ) with L u is a typical flow lengthscale) is greater than 0.1. In our simulations, Ro l ≈ 0.03 for Ra = 2800 and Ek = 10 −5 , so we expect that higher values of Ra are required to enter the multipolar regime in this case, despite C ωz being less than 0.5. The cases with B 0 > 0 show a gradual degradation of the relative dipole strength when Ra increases, which is more pronounced for B 0 = 4. However, the field remains overwhelmingly dipolar even at large Ra where the induced field is of the same order of magnitude as the imposed field. In these cases too, Ro l is still less than 0.1.

Morphology of the magnetic field
We now analyse the morphology of the magnetic field, starting from the axisymmetric magnetic field. Figure 18 shows meridional slices of the axisymmetric azimuthal magnetic field and the poloidal magnetic field lines. In the magnetoconvective cases, the poloidal field lines remain largely z-invariant, indicating that the imposed field is relatively unperturbed by the induced poloidal field, with the most significant variations observed at larger Ra. For B 0 = 4 and Ra = 1200, the induced field has a larger amplitude than the imposed field (see figure 15(b)) but, even then, the magnetic field largely retains the shape of the imposed magnetic field. For B 0 = 4, the azimuthal field consists of a single ring located outside the tangent cylinder and of opposite sign in each hemisphere. For Ra = 200, the azimuthal field can be related to the zonal flow shown in figure 12 via an omega effect. Only weak azimuthal field is present inside the tangent cylinder, which is consistent with the absence of zonal flow there (and despite the presence of convection for Ra ≥ 500). By contrast, for B 0 = 0, a strong azimuthal field is produced on the tangent cylinder and diffuses inwards (there is no convection inside the tangent cylinder for Ra = 500 in this case). For B 0 = 10, the azimuthal field is very weak, even for higher values of Ra, and multiple patches are seen in both hemispheres.
To determine whether the induced magnetic field reinforces or diminishes the imposed magnetic field, we plot in figure 19 the induced component of the axisymmetric axial magnetic field, B z − B 0 . Only the northern hemisphere is shown as B z is largely symmetric about the equator. For B 0 > 0, the axial magnetic field is strongest at higher latitudes, and of the sense to reinforce the imposed field. (The sense for B 0 = 0 is coincidental, as the opposite sense is an equally valid solution.) In the equatorial regions near the outer boundary (and in the exterior), the induced dipole field results in locally negative induced B z values; these are always weaker than the imposed field, however, so the total axisymmetric B z is always positive. For B 0 = 4, there is a significant octopole component of induced field close to onset, resulting in opposite polarity fields near the equator. At Ra = 500, the induced field inside the tangent cylinder is a significant fraction of the imposed field. Overall, for both B 0 = 4 and 10, the induced axial dipole therefore reinforces the imposed field. Figure 20 shows equatorial slices of the induced axial magnetic field (i.e. B z − B 0 ) in the equatorial plane, superimposed on the axial vorticity. For B 0 = 0, we show the nearest case to the dynamo onset, Ra = 500 ≈ 4.7Ra 0 c . In this case, the maximum of B z is located within the anticyclones, which is a familiar feature of convective dynamos (e.g. Jones 2011). In the magnetoconvective cases for Rayleigh numbers close to onset, the axial magnetic field is concentrated between cyclones and anticyclones, with regions of positive induced axial fields (which reinforce the imposed field) correlated with cold (inward) radial flow. At larger Rayleigh numbers (Ra = 350) however, the axial magnetic field is concentrated in the anticyclones in the equatorial plane, similarly to the dynamo case. This observation is in agreement with the study of Sakuraba and Kono (2000) for Λ 0 > 2, which emphasises that the magnetic flux concentration inside anticyclones leads to an asymmetry between cyclones and anticyclones with enlarged anticyclones. They argue that the magnetic confinement is due to convergent flows towards the core of the anticyclones in the equatorial plane that compensate the Ekman pumping at the outer boundary.

Generation of the axisymmetric poloidal field
In order to determine which lengthscales play a major role in the generation of the axisymmetric poloidal magnetic field, we study the contributions to the electromotive force (e.m.f. ). The equation for the poloidal magnetic component B p is obtained by taking the dot product of the magnetic induction equation with the vector r, where we have used that rB r = L 2 B p and the L 2 operator isolates the angular part of the Laplacian, Applied to a spherical harmonic Y m l , it gives L 2 Y m l = l(l + 1)Y m l . The e.m.f. is E = u × B. From equation (33), the axisymmetric poloidal component B p is driven by the where the superscript indicates a component of the velocity or magnetic field corresponding to a given m. Figure 21 shows the r.m.s. value of E m φ as a function of m for m ∈ [0, 20] for Ra close to onset and at Ra = 500. Figure 22 shows the corresponding meridional slices of E m φ for the dominant value(s) of m and also the sum m=20 m=0 E m φ . E φ is mainly symmetric about the equator in all cases, so we only show the northern hemisphere. For comparison, we also plot the azimuthal component of the axisymmetric electric current (which is related to the mean e.m.f. via Ohm's law) For B 0 = 0, the largest contribution to the mean e.m.f. is produced by the modes m ≥ 8, all contributing at a similar level until m = 20. The large-scale modes 1 ≤ m ≤ 7 have relatively weak kinetic energy (see the kinetic energy spectrum in figure 8(b)), so it is perhaps unsurprising that their contribution to the mean e.m.f. is small. The sum of the mean e.m.f. produced by the modes 0 ≤ m ≤ 20 is very similar to j φ , so the small scales corresponding to m > 20 only play a minor role in the generation of the axisymmetric poloidal field. In fact, this is true for all the cases shown in figure 22.
For B 0 = 4, the main contribution comes from m = 0 and m = 4 at onset, with m = 4 generating a mean e.m.f. that opposes that of m = 0 (which is largely due to the interaction of the axisymmetric meridional circulation with the imposed field), but cannot entirely match it. The net mean e.m.f. is thus similar to that of m = 0 with a reduced amplitude. It produces a ring of negative j φ in the equatorial region, topped by a ring of positive j φ at higher latitudes. Positive (negative) induced B z is produced inside (outside) the ring of positive j φ and vice-versa (see figure 19). For Ra = 500, modes up to m = 8 contribute significantly to E φ . Here again the non-axisymmetric modes generate a mean e.m.f. in opposition to that of m = 0. The net e.m.f. now produces mainly a ring of positive j φ whose maximum is located at mid-latitudes near the outer boundary. As a result, the induced B z is positive in most of the domain and only negative near the outer boundary at low latitudes. For Ra away from the onset, the large-scale non-axisymmetric modes that dominate the kinetic energy spectrum ( figure 8(b)) therefore largely contribute to the generation of an induced field that reinforces the imposed field.
At B 0 = 10 near onset, the dominant contribution comes from m = 9 (which dominates the kinetic energy spectrum). The resulting j φ is positive and contained in a narrow band outside the tangent cylinder, which corresponds to the location of the columnar flow (see figures 2 and 20). Accordingly, the induced B z is mainly positive inside the tangent cylinder and negative outside. At higher Ra, the dominant non-axisymmetric mode is now m = 11 (also the dominant non-axisymmetric mode of the kinetic energy spectrum), with m = 0 contributing at a similar level. These two modes produce e.m.f. that are in opposition in most of the domain. As the flow structure extends further in cylindrical radius when Ra increases (see figure 20), the resulting positive j φ now fills most of the domain outside the tangent cylinder, producing a positive induced B z throughout most of the domain. Interestingly, the net mean e.m.f. is fairly similar for B 0 = 4 and B 0 = 10 at Ra = 500 (hence similar induced B z in figure 19), although the contributions from individual modes are very different. This highlights a notable difference in field generation mechanisms between B 0 = 4 and B 0 = 10, although the net outcome in both cases is the production of an axisymmetric poloidal field reinforcing the imposed field in the interior.

Discussion
We have studied the effects of an imposed axial magnetic field B 0 = B 0 e z on rotating spherical convection at onset and for supercritical values of the Rayleigh number. Our parameter survey included both Ek = 10 −4 and Ek = 10 −5 , for a range of values of Ra/Ra c varying from 1 to 50, with fixed Pr = Pm = 1, and with fixed temperature, no-slip, and electrically insulating boundary conditions. We found significant changes in the onset of convection when B 0 O(1) (where B 0 = 1 corresponds to an Elsasser number of unity), with a decrease in the critical Rayleigh number Ra c and azimuthal wavenumber m c that is more pronounced at small Ek. The critical Rayleigh number increases again for larger values of B 0 ∼ O(10), becoming independent of the choice of Ek. The structure of the flow at onset is mainly columnar, regardless of the choice of B 0 . The primary force balance gradually changes as B 0 increases, and the Lorentz force becomes one of the dominant forces balancing the Coriolis and buoyancy forces. Interestingly, we find that the primary force balance varies locally: for B 0 = 10, the buoyancy and Lorentz forces balance in the equatorial region, but the balance is mainly between the Coriolis and Lorentz forces near the outer surface. The change in the force balance is evident in the wave propagation, which is prograde and matches the frequency of thermal Rossby waves for B 0 < O(1), and subsequently slows down as B 0 increases, until reversing direction at B 0 = 10 for small Ek.
For higher values of Ra, the columnar flows break down, with the earliest collapse seen for B 0 = 4. In this case, more efficient transport of heat is seen in the equatorial region, as the Lorentz force offsets the Coriolis force here. The improved convective heat transfer leads to strong latitudinal variations of the temperature resulting in a strong thermal wind. This in turn has led to zonal flows in our magnetoconvective simulations, despite the equivalent dynamo simulations normally having weak zonal flows. These zonal flows can distort the poloidal component of the magnetic field, to produce a stronger axisymmetric toroidal field. For larger values of B 0 , zonal flows are suppressed by the stronger axial field. Yadav et al. (2016b) found that the presence of dominant zonal flows (and hence large shear) in non-magnetic convection simulations reduces the efficiency of the convective heat transfer compared with dynamo cases, where the zonal flows are drastically reduced. This occurs when using stress-free boundary conditions, which favour the development of large zonal flows. Here, on the contrary, we have a strong thermal wind combined with efficient heat transfer. However the zonal flows are never dominant over the convective flows in our simulations because the zonal flow amplitude is limited by the use of no-slip boundaries. Additionally, Yadav et al. (2016b) found that no-slip dynamo simulations have significantly improved heat transfer compared with the non-magnetic cases when the Elsasser number of the self-sustained magnetic field is order unity, which is consistent with our results.
Near onset, the induced magnetic field generated by magnetoconvection is much smaller than the imposed field. The existence of subcritical behaviour, where a convective solution that is obtained by switching off the imposed field exists at the same parameters as the trivial conducting state, would require that the induced magnetic field is sufficiently strong to modify the flow effectively. Here we found that the field needs to be order unity to significantly influence the flow, so we expect that subcritical solutions would be difficult to maintain near onset. Indeed, we carried out multiple attempts to maintain convective dynamos below the onset of convection by gradually reducing the imposed field, but all were unsuccessful. However, for B 0 = 4, the induced field exceeds the imposed field for Rayleigh numbers approximately ten times critical. The induced field reinforces the imposed field in most of the domain, and the magnetic field largely retains the z-invariance of the imposed field, even at the largest Rayleigh numbers. The magnetically-modified flows are therefore able to generate a large-scale field with an intensity and morphology comparable to the imposed field. As a result, the intermittent behaviour observed in the planar dynamo simulations of Stellmach and Hansen (2004), Cooper et al. (2020) and suggested by the magnetoconvection results of Gubbins (2000a, 2000b) might be avoided in this case. Additionally, hysteretic behaviour at high Rayleigh numbers might be possible: for Rayleigh numbers greater than the value required for the dynamo onset, a "strong" dynamo solution (obtained by switching off the imposed field) might exist at the same parameters as a "weak" dynamo solution (obtained by small initial perturbation). This was indeed observed by Sarson et al. (1997Sarson et al. ( , 1999 and Sakuraba and Kono (2000) in this regime. We have also attempted to obtain hysteresis for B 0 = 4 and Ra ≥ 750, but this has been unfruitful so far and will require further study.
The large-scale field generation is greatly affected by changes in the imposed field strength. In our dynamo simulations at B 0 = 0, the axisymmetric poloidal field is generated by the "small-scale" convective flows that are relatively unaffected by the Lorentz force (corresponding to azimuthal orders m ∈ [8, 20] for Ek = 10 −5 , but we would expect these scales to become smaller for decreasing Ek based on linear stability). A clear change occurs for B 0 = 4, where only the large-scale flows (with m < 8) are significant in the production of the axisymmetric poloidal field. These scales carry most of the kinetic energy and are greatly affected by the Lorentz force. Notably, they lack attributes that are sometimes associated with the generation of axial dipoles (e.g. Olson et al. 1999, Sreenivasan andJones 2011): their relative kinetic helicity and columnarity are significantly reduced compared with cases at B 0 = 0. Another change occurs at larger B 0 (B 0 = 10), where a single non-axisymmetric mode is largely responsible for the production of the axisymmetric poloidal field together with contributions from the axisymmetric mode (at least up to Ra = 500, where this mode remains close to the value of the dominant wavenumber at onset). Our simulations are performed at moderate values of Ek, so the distinction between "small" and "large" scales is not strongly marked, but we would expect this distinction to become very pronounced in planetary core conditions at smaller Ek.
Magnetoconvection calculations help to elucidate the mechanisms for the production of magnetic fields in planets embedded in external fields of various strengths (e.g. the Jovian satellites and some exoplanets), as well as general aspects of field-generation in a strong-field regime. Our calculations complement the more prevalent studies of magnetoconvection with an imposed toroidal field. Some important differences are notedincluding the imposed field strength required to markedly affect the flow (Fearn 1979) and the effect of the imposed axial field to reduce the relative kinetic helicity (e.g. see Sreenivasan and Jones 2011), with possibly significant implications for magnetic field generation -which warrant further work to elaborate to what extent the effects noted here rely on the imposed axial field geometry. A disadvantage of using a uniform magnetic field is that the whole domain is subject to the presence of the ambient field, whereas recent dynamo simulations show that the magnetic field is heterogeneous within the bulk of the domain (Schaeffer et al. 2017). The presence of magnetic-free regions have not been considered here, but could be the subject of future magnetoconvection studies. Other work could further investigate different boundary conditions or parameters, or imposing an axial field that is not uniform (Sreenivasan and Gopinath 2017).
Work is currently in progress to extend this study to variable Pm, for values both smaller and larger than the Pm = 1 case considered here. The case Pm < 1 is relevant to planetary interiors (where Pm 1, based on molecular diffusivities); but as noted by Dormy (2016), the case Pm > 1 may allow calculations at relatively high Ek to attain a strong-field solution branch, more relevant to the expected magnetostrophic balance, than the weak-field and multipolar solutions found here for Pm = 1 (when considering B 0 = 0). The combined effects of varying Pm and B 0 are therefore of considerable interest for understanding planetary dynamos.
Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). DiRAC is part of the National e-Infrastructure. We thank the reviewers and editor whose comments helped clarify and improve the manuscript.

Disclosure statement
No potential conflict of interest was reported by the author(s). where we used L 2 Y m l = l(l + 1)Y m l , with L 2 the angular part of the Laplacian introduced in section 5.3. The solutions are of the form g m l = C α r α . Substituting in equation (A.6) gives α(α + 1) = l(l + 1), which has for solutions α = l (corresponding to external sources) and α = −l − 1 (corresponding to internal sources).
We consider the case where an external source at r > r o produces a uniform axial field, where we have introduced the relevant factors to match our choice of units and