*Article* **Onset of Inertial Magnetoconvection in Rotating Fluid Spheres**

**Radostin D. Simitev 1,\* and Friedrich H. Busse <sup>2</sup>**


**Abstract:** The onset of convection in the form of magneto-inertial waves in a rotating fluid sphere permeated by a constant axial electric current is studied in this paper. Thermo-inertial convection is a distinctive flow regime on the border between rotating thermal convection and wave propagation. It occurs in astrophysical and geophysical contexts where self-sustained or external magnetic fields are commonly present. To investigate the onset of motion, a perturbation method is used here with an inviscid balance in the leading order and a buoyancy force acting against weak viscous dissipation in the next order of approximation. Analytical evaluation of constituent integral quantities is enabled by applying a Green's function method for the exact solution of the heat equation following our earlier non-magnetic analysis. Results for the case of thermally infinitely conducting boundaries and for the case of nearly thermally insulating boundaries are obtained. In both cases, explicit expressions for the dependence of the Rayleigh number on the azimuthal wavenumber are derived in the limit of high thermal diffusivity. It is found that an imposed azimuthal magnetic field exerts a stabilizing influence on the onset of inertial convection and as a consequence magneto-inertial convection with azimuthal wave number of unity is generally preferred.

**Keywords:** rotating thermal magnetoconvection; linear onset; sphere

**Citation:** Simitev, R.D.; Busse, F.H. Onset of Inertial Magnetoconvection in Rotating Fluid Spheres. *Fluids* **2021**, *6*, 41. https://doi.org/doi:10.3390/ fluids6010041

Received: 15 December 2020 Accepted: 10 January 2021 Published: 13 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

## **1. Introduction**

Buoyancy-driven motions of rotating, electrically conducting fluids in the presence of magnetic fields represent a fundamental aspect of the dynamics of stellar and planetary interiors, see, e.g., in [1–4]. The problem of magnetic field generated and sustained by convection is rather difficult to attack both analytically and numerically because of its essential nonlinearity and scale separation [5,6]. Valuable insights can be gained by studying magnetoconvection, the simpler case of an imposed magnetic field, which has received much attention ever since the early work of Chandrasekhar [7], see in [8,9]. For instance, the propagation of rotating magnetoconvection modes excited in the deep convective region of the Earth's core has been proposed as a possible mechanism for explaining features of observed longitudinal geomagnetic drifts [10,11], see also the recent review of Finlay et al. [12]. A rather detailed classification of magnetoconvection waves in a rotating cylindrical annulus has been recently attempted by Hori et al. [13] and the authors of [14] who proceeded further to make useful comparisons with nonlinear spherical dynamo simulations and to provide estimates for the strength of the "hidden" azimuthal part of the magnetic field within the core. These authors used the rotating annulus model of Busse [15,16] and only considered values of the Prandtl number of the order unity. However, both spherical geometry as well as small values of the Prandtl number are essential features of a planetary or a stellar interior [17]. At sufficiently small values of the Prandtl number, a different style of convection exists that is sometimes called inertial or equatorially-attached convection or thermo-inertial waves [18–21]. In this limit, convection oscillates so fast that the viscous force does not enter the leading-order balance. The latter is then reduced to the Poincaré equation in a rotating spherical system [22–24]. On the longer time scale of the next order of approximation the buoyancy force maintains convection against the weak viscous dissipation. This regime of

convection thus represents a transition between thermal convection and wave propagation in rapidly rotating geometries. It is important to understand how this regime of inertial convection is affected by an imposed magnetic field.

With this in mind, we study in the present paper the onset of magneto-inertialconvection. In particular, we consider a rotating fluid sphere permeated by a constant axial electric current as proposed by Malkus [11] in the limit of low viscosity and high thermal diffusivity (small Prandtl number). A similarly configured problem was also investigated by Zhang and Busse [25] who derived an explicit dependence of the critical Rayleigh number on the imposed field strength but were not able to obtain an explicit dependence on the azimuthal wavenumber of the modes as this requires the evaluation of a volume integral of the temperature perturbation [25]. In an earlier paper, we proposed a Green's function method for the exact solution of the heat equation [26] which then allowed the analytical evaluation of the integral quantities needed to find a fully explicit expressions for the critical Rayleigh number and frequency for the onset of convection and to study mode competition. Here, we apply the same approach to the case of magneto-inertial convection and we consider both value and flux boundary conditions for the temperature.

In the following we start with the mathematical formulation of the problem in Section 2. The special limit of a high ratio of thermal to magnetic diffusivity will be treated in Section 3. The general case requires the symbolic evaluation of lengthy analytical expressions and will be presented in Section 4. A discussion of the results and an outlook on related problems will be given in the final Section 5 of the paper.

#### **2. Mathematical Formulation of the Problem**

We consider a homogeneously heated and self-gravitating sphere as illustrated in Figure 1. The sphere is filled with incompressible and electrically conducting fluid characterized by its magnetic diffusivity *η*, kinematic viscosity *ν*, thermal diffusivity *κ*, and density *\$*. The sphere is rotating with a constant angular velocity Ω**k** where **k** is the axial unit vector. The gravity field is given by **g** = −*gr*0**r**, where **r** is the position vector with respect to the center of the sphere, *r* is its length measured in fractions of the radius *r*<sup>0</sup> of the sphere, and *g* is the amplitude of the gravitational acceleration. Following Malkus [11], we assume that the fluid sphere is permeated by a toroidal magnetic field **B** ∼ **k** × **r**. As the Lorentz force like the centrifugal force can be balanced by the pressure gradient, a static state of no motion exists with the temperature distribution *T<sup>S</sup>* = *T*<sup>0</sup> − *βr* 2 0 *r* <sup>2</sup>/2. We employ the Boussinesq approximation and assume constant material properties *η*, *ν*, *κ*, and *\$* everywhere except in the buoyancy term where the density is assumed to have a linear dependence on temperature with a coefficient of thermal expansion *α* ≡ (*d\$*/*dT*)/*\$* = const.

**Figure 1.** Geometrical configuration of the problem. A part of the outer spherical surface is removed to expose the interior of sphere to which the conducting fluid is confined.

In order to study the onset of magnetoconvection in this system, we consider the linearized momentum, magnetic induction, heat, continuity, and solenoidality equations:

$$
\partial\_l \mathbf{\tilde{u}} + \tau \mathbf{k} \times \mathbf{\tilde{u}} + \nabla (\pi - \mathbf{\tilde{b}} \cdot \mathbf{j} \times \mathbf{r}) + (\mathbf{j} \times \mathbf{r}) \cdot \nabla \mathbf{\tilde{b}} - \mathbf{j} \times \mathbf{\tilde{b}} = \boldsymbol{\Theta} \mathbf{r} + P\_m \nabla^2 \mathbf{\tilde{u}},\tag{1a}
$$

$$
\partial\_t \tilde{\mathbf{b}} - (\mathbf{j} \times \mathbf{r}) \cdot \nabla \tilde{\mathbf{u}} + \mathbf{j} \times \tilde{\mathbf{u}} = \nabla^2 \tilde{\mathbf{b}},\tag{1b}
$$

$$
\hat{\mathbf{R}}\,\mathbf{r}\cdot\tilde{\mathbf{u}} + \nabla^2 \tilde{\Theta} - (P/P\_m)\partial\_l \tilde{\Theta} = 0,\tag{1c}
$$

$$
\nabla \cdot \tilde{\mathbf{u}} = 0, \quad \nabla \cdot \tilde{\mathbf{b}} = 0,\tag{1d}
$$

respectively, that govern the evolution of infinitesimal velocity perturbations **u**˜, temperature perturbations Θ˜ , and magnetic field perturbations **b**˜ away from the static state. The equations have been non-dimensionalized using the radius *r*<sup>0</sup> as a unit of length, *r* 2 0 /*η* as a unit of time, *η* <sup>2</sup>/*gαr* 4 0 as a unit of temperature, and √*µ\$η*/*r*<sup>0</sup> as a unit of magnetic flux density. The dimensionless magnetic field takes the form **<sup>j</sup>** <sup>×</sup> **<sup>r</sup>** <sup>+</sup> **<sup>b</sup>**˜, where **<sup>j</sup>** <sup>=</sup> *<sup>j</sup>***<sup>k</sup>** is the vector of the density of the imposed electric current. The problem is then characterized by five dimensionless parameters, namely, the Rayleigh number, the Coriolis parameter, the Prandtl number, the magnetic Prandtl number *Pm*, and the non-dimensional current density given by

$$\vec{R} = \frac{\arg \beta r\_0^6}{\eta \kappa}, \ \tau = \frac{2\Omega r\_0^2}{\eta}, \ P = \frac{\nu}{\kappa'} \ P\_m = \frac{\nu}{\eta'} \ j\_\nu \tag{2}$$

respectively. In fact, in the results obtained below the two Prandtl numbers enter only as their ratio *S* = *P*/*P<sup>m</sup>* = *η*/*κ*. To signify that in our definition of the Rayleigh number the magnetic diffusivity replaces the kinematic viscosity we have attached a hat to *R*ˆ.

#### **3. Perturbation Analysis Results**

Without loss of generality we assume that the velocity, the magnetic field, and the temperature perturbations have an exponential dependence on time *t* and on the azimuthal angle *φ*. Further, as both the velocity field and the magnetic field are solenoidal we use the poloidal-toroidal decomposition

$$\overrightarrow{\mathbf{u}} = \mathbf{u} \exp\left(i(\omega \tau t + m\phi)\right) = \left(\nabla \times (\nabla v \times \mathbf{r}) + \nabla w \times \mathbf{r}\right) \exp\left(i(\omega \tau t + m\phi)\right), \tag{3a}$$

$$\tilde{\mathbf{b}} = \mathbf{b} \exp\left(i(\omega \tau t + m\phi)\right) = \left(\nabla \times (\nabla h \times \mathbf{r}) + \nabla \mathbf{g} \times \mathbf{r}\right) \exp\left(i(\omega \tau t + m\phi)\right), \tag{3b}$$

$$\mathfrak{G} = \Theta \exp\left(i(\omega \,\tau t + m\phi)\right). \tag{3c}$$

Equation (1b) can now be written in the form

$$\mathbf{b} = \frac{m\gamma}{\omega}\mathbf{u} - \frac{i}{\omega\tau}\nabla^2\mathbf{b},\tag{4}$$

where the parameter *γ* is defined as *γ* = *j*/*τ*, and in the ∇-operator the *φ*-derivative is replaced by its eigenfactor *im*. This allows us to transform Equation (1a) in the form

$$\begin{split} i\omega \left(1 - \frac{m^2}{\omega^2} \gamma^2 \right) \mathbf{u} + \left(1 - \frac{m}{\omega} \gamma^2 \right) \mathbf{k} \times \mathbf{u} - \nabla \pi &= \frac{1}{\tau} \Theta \mathbf{r} + \frac{P\_m}{\tau} \nabla^2 \mathbf{u} + \frac{m^2 \gamma^2}{\omega^2 \tau} \nabla^2 \mathbf{u} \\ + \frac{2m\gamma^2}{i\omega^2 \tau} \mathbf{k} \times \nabla^2 \mathbf{u} + \frac{m\gamma}{\omega \tau} \nabla^2 \mathbf{b}\_b + \frac{2\gamma}{i\omega \tau} \mathbf{k} \times \nabla^2 \mathbf{b}\_b + \frac{P\_m}{\tau} \nabla^2 \mathbf{u}\_{b\prime} \end{split} \tag{5}$$

where *π*ˇ is the effective pressure. In Equation (5), the magnetic field **b** appears only in the form of the boundary layer correction **b***<sup>b</sup>* , which is required as the basic dissipationless solution does not satisfy all boundary conditions [25]. For the same reason, the Ekman layer correction **u***<sup>b</sup>* must be introduced [22].

Following the procedure of earlier papers [22,26], we use a perturbation approach and solve Equation (1a) in the limit of large *τ*, using the ansatz

$$\mathbf{u} = \mathbf{u}\_0 + \tau^{-1}\mathbf{u}\_1 + \dots, \quad \omega = \omega\_0 + \tau^{-1}\omega\_1 + \dots, \quad \mathbf{b} = \mathbf{b}\_0 + \tau^{-1}\mathbf{b}\_1 + \dots \tag{6}$$

The heat equation is solved unperturbed.

#### *3.1. Zeroth-Order Approximation*

In the following we shall assume the limit of large *τ* such that in zeroth order of approximation the right hand side of Equation (5) can be neglected. The left hand side together with the condition ∇ · **u** = 0 is of the same form as the equation for inertial modes [22,26]. In the nonmagnetic case, the inertial modes corresponding to the sectorial spherical harmonics yield the lowest critical Rayleigh numbers for the onset of convection [26]. We shall assume that this property continues to hold as long as the parameter *γ* is sufficiently small so that the nonmagnetic limit is approached in the left-hand side of Equation (5). The sectorial inertial modes are given by

$$w\_0 = P\_m^m(\cos \theta) f(r), \quad w\_0 = P\_{m+1}^m(\cos \theta) \psi(r), \tag{7a}$$

with

$$f(r) = r^m - r^{m+2}, \qquad \psi(r) = r^{m+1} \frac{2im(m+2)}{(2m+1)(\lambda\_0(m^2+3m+2)-m)}, \tag{7b}$$

where *λ*<sup>0</sup>

$$
\lambda\_0 = \frac{1}{m+2} \left( 1 \pm \sqrt{\frac{m^2 + 4m + 3}{2m + 3}} \right) \tag{7c}
$$

is the frequency of the inertial modes. The sectorial magneto-inertial modes are then described by the same velocity field (7a) and by a magnetic field **b**<sup>0</sup> = *mγ***u**0/*ω*0. In the above expressions, the subscript 0 refers to the dissipationless solution of Equation (1). The frequency *ω*<sup>0</sup> of the magneto-inertial waves is determined by

$$
\lambda\_0 = \frac{\omega\_0^2 - m^2 \gamma^2}{\omega\_0 - m\gamma^2}.\tag{8}
$$

which yields

$$
\omega\_0 = \frac{\lambda\_0}{2} \pm \sqrt{\frac{\lambda\_0^2}{4} + m\gamma^2 (m - \lambda\_0)}.\tag{9}
$$

With account of (7c), this dispersion relation allows for a total of four different frequencies *ω*0. For small values of *γ* 2 , these are given by

$$\begin{split} \omega\_{01,2} &= \frac{1}{m+2} \left( 1 \pm \sqrt{\frac{m^2 + 4m + 3}{2m + 3}} \right) \\ &+ m^2 \gamma^2 (m+2) \left( 1 \pm \sqrt{\frac{m^2 + 4m + 3}{2m + 3}} \right)^{-1} - m \gamma^2, \end{split} \tag{10a}$$

$$
\omega\_{03,4} = -m^2 \gamma^2 (m+2) \left( 1 \pm \sqrt{\frac{m^2 + 4m + 3}{2m + 3}} \right)^{-1} + m\gamma^2. \tag{10b}
$$

The upper sign in expression (10a) refers to retrogradely propagating modified inertial waves, while the lower sign corresponds to the progradely traveling variety. The effect of the magnetic field tends to increase the absolute value of the frequency in both cases. Expression (10b) describes the dispersion of the slow magnetic waves. The upper sign refers to the progradely traveling modified Alfven waves and the lower sign corresponds to retrogradely propagating modified Alfven waves.

#### *3.2. First-Order Approximation*

The magneto-inertial waves described by expressions (7a) satisfy the condition that the normal component of the velocity field vanishes at the boundary. This property implies that the normal component of the magnetic field vanishes there as well. Additional boundary conditions must be specified when the full dissipative problem described by (5) is considered. We shall assume a stress-free boundary with either a fixed temperature (case A) or a thermally insulating boundary (case B),

$$\mathbf{r} \cdot \mathbf{u} = \mathbf{r} \cdot \nabla (\mathbf{r} \times \mathbf{u}) / r^2 = 0 \quad \text{and} \quad \left\{ \begin{array}{ll} \Theta = 0 & \text{(case A)}\\ \partial\_r \Theta = 0 & \text{(case B)} \end{array} \right\} \quad \text{at} \quad r = 1. \tag{11}$$

Additionally, we shall assume an electrically insulating exterior of the sphere which requires

$$\mathbf{g} = \mathbf{0} \quad \text{at} \quad r = 1 \tag{12}$$

and the matching of the poloidal magnetic field to a potential field outside the sphere.

After the ansatz (6) has been inserted into Equation (5) such that terms with **u**<sup>1</sup> appear on the left hand side, while those with **u**<sup>0</sup> and *ω*<sup>0</sup> appear on the right hand side, we obtain the solvability condition for the equation for **u**<sup>1</sup> by multiplying it with **u** ∗ 0 and averaging it over the fluid sphere,

$$i\omega\_1 \left< |\mathbf{u}\_0|^2 \right> \left( 1 + \left( \frac{m^2}{\omega\_0^2} - \frac{m(\omega\_0^2 - m^2\gamma^2)}{\omega\_0^2(\omega\_0 - m\gamma^2)} \right) \gamma^2 \right) \tag{13}$$

$$= \left< \Theta \mathbf{r} \cdot \mathbf{u}\_0^\* \right> + \left( \langle \mathbf{u}\_0^\* \cdot \nabla^2 \mathbf{u}\_0 \rangle \frac{m\gamma}{\omega\_0} + \langle \mathbf{u}\_0^\* \cdot \nabla^2 \mathbf{b}\_{0b} \rangle \right) \left( \frac{m}{\omega} - \frac{\omega\_0^2 - m^2\gamma^2}{\omega\_0^2 - m\omega\_0\gamma^2} \right) \frac{\gamma}{\tau} \gamma$$

where the brackets h...i indicate the average over the fluid sphere and the ∗ indicates the complex conjugate. We have neglected all terms connected with viscous dissipation, i.e., we have assumed the vanishing of *Pm*, as we wish to focus on the effect of ohmic dissipation. The effects of viscous dissipation have been dealt with in the earlier paper [26]. Because h**u** ∗ 0 · ∇2**u**0<sup>i</sup> vanishes, as demonstrated in [27], we must consider only the influence of the boundary layer magnetic field **b**0*<sup>b</sup>* . It is determined by the equation

$$i\omega\_0 \tau \mathbf{b}\_{0\emptyset} = \nabla^2 \mathbf{b}\_{0\emptyset}.\tag{14}$$

As the solutions of this equation are characterized by gradients of the order √ *τ*, the boundary layer correction needed for the poloidal component is of the order √ *τ* smaller than the correction needed for the toroidal component. For large *τ* we need to take into account only the contribution *g*0*<sup>b</sup>* given by

$$g\_{0b} = -g\_0(r=1) \exp\left(-(1+is)(1-r)\sqrt{|\omega\_0|\tau/2}\right)$$

$$= -\frac{m\gamma}{\omega\_0} w\_0(r=1) \exp\left(-(1+is)(1-r)\sqrt{|\omega\_0|\tau/2}\right). \tag{15}$$

where *s* denotes the sign of *ω*0. The solvability condition thus becomes reduced to

$$\begin{split} \operatorname{Li}\_{\mathsf{T}}\left(|\mathbf{u}\_{0}|^{2}\right) \Bigg( \mathbbm{1} + \left(\frac{m\gamma^{2}(m-\omega\_{0})}{\omega\_{0}(\omega\_{0}-m\gamma^{2})}\right) \Bigg) \\ &= \frac{1}{\mathsf{T}} \langle \Theta \mathbf{r} \cdot \mathbf{u}\_{0}^{\*} \rangle - \frac{3}{2} \frac{m\gamma^{2}(m-\omega\_{0})(s+i)}{(\omega\_{0}-m\gamma^{2})\sqrt{2|\omega\_{0}|\pi}} \int\_{-1}^{1} |P\_{m}^{m+1}|^{2} \mathbf{d}(\cos\theta) \\ &\times (m+1)(m+2) \left| \frac{2m(m+2)}{(2m+1)\left(\frac{\omega\_{0}^{2}-m^{2}\gamma^{2}}{\omega\_{0}-m\gamma^{2}}(m+1)(m+2)-m\right)} \right|^{2} .\end{split}$$
