2 Mathematical model

Our model of a stably stratified stellar radiative zone is based on the anelastic approximation (Gough Reference Gough1969; Braginsky & Roberts Reference Braginsky and Roberts1995; Lantz & Fan Reference Lantz and Fan1999) which is widely adopted for numerical simulations of convection in Solar and stellar interiors (Jones et al. Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011). Accordingly, we consider a perfect gas confined to a spherical shell rotating with a fixed angular velocity $\unicode[STIX]{x1D6FA}\hat{\boldsymbol{k}}$ , where $\hat{\boldsymbol{k}}$ is the unit vector in the direction of the axis of rotation. A positive entropy contrast $\unicode[STIX]{x0394}S$ is imposed between its outer and inner surfaces at radii $r_{o}$ and $r_{i}$ , respectively. We assume a gravity field proportional to $g/r^{2}$ . To justify this choice consider the Sun, the star with the best-known physical properties. The Solar density drops from $150$ at the centre to $20~\text{g}~\text{cm}^{-3}$ at the core-radiative zone boundary (at $0.25R_{\odot }$ ) to only 0.2 at the tachocline (at $0.7R_{\odot }$ ). A crude piecewise linear interpolation shows that most of the mass is concentrated within the core. In this setting a hydrostatic polytropic reference state exists with profiles of density, temperature and pressure given by the expressions

(2.1a-c ) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70C}_{c}\unicode[STIX]{x1D701}^{n},\quad \bar{T}=T_{c}\unicode[STIX]{x1D701},\quad \bar{P}=P_{c}\unicode[STIX]{x1D701}^{n+1},\end{eqnarray}$$

respectively, where

(2.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D701}=c_{0}+{\displaystyle \frac{c_{1}}{r}},\quad \unicode[STIX]{x1D701}_{o}={\displaystyle \frac{\unicode[STIX]{x1D702}+1}{\unicode[STIX]{x1D702}\exp (N_{\unicode[STIX]{x1D70C}}/n)+1}},\\ c_{0}={\displaystyle \frac{2\unicode[STIX]{x1D701}_{o}-\unicode[STIX]{x1D702}-1}{1-\unicode[STIX]{x1D702}}},\quad c_{1}={\displaystyle \frac{(1+\unicode[STIX]{x1D702})(1-\unicode[STIX]{x1D701}_{o})}{1-\unicode[STIX]{x1D702}}}^{2},\end{array}\right\}\end{eqnarray}$$

see (Jones et al. Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011). The parameters $\unicode[STIX]{x1D70C}_{c}$ , $P_{c}$ and $T_{c}$ are reference values of density, pressure and temperature at mid-shell. The gas polytropic index $n$ , the density scale height $N_{\unicode[STIX]{x1D70C}}$ and the shell radius ratio $\unicode[STIX]{x1D702}$ are defined further below. Since a rigidly rotating stellar interior cannot exist, the deviation from such a state is described by the baroclinic parameter $\mathit{Z}$ , introduced below. This parameter is associated with the centrifugal force which balances the baroclinic torques (Busse Reference Busse1982). Following Rieutord (Reference Rieutord2006) we neglect the distortion of the isopycnals caused by the centrifugal force. The governing anelastic equations of continuity, momentum, energy (entropy) and magnetic induction assume the form

(2.3a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D70C}}\boldsymbol{u} & = & \displaystyle 0,\end{eqnarray}$$

(2.3b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B} & = & \displaystyle 0,\end{eqnarray}$$

(2.3c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+(\unicode[STIX]{x1D735}\times \boldsymbol{u})\times \boldsymbol{u} & = & \displaystyle -\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}-{\displaystyle \frac{2}{\mathit{Ek}}}(\hat{\boldsymbol{k}}\times \boldsymbol{u})+{\displaystyle \frac{R}{\mathit{Pr}}}{\displaystyle \frac{S}{r^{2}}}\hat{\boldsymbol{r}}+\boldsymbol{F}_{\unicode[STIX]{x1D708}}\nonumber\\ \displaystyle & & \displaystyle +\,{\displaystyle \frac{1}{\bar{\unicode[STIX]{x1D70C}}}}(\unicode[STIX]{x1D735}\times \boldsymbol{B})\times \boldsymbol{B}-\mathit{Z}(\overline{S}+S)\;\hat{\boldsymbol{k}}\times (\boldsymbol{r}\times \hat{\boldsymbol{k}}),\end{eqnarray}$$

(2.3d ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}S+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\overline{S}+S) & = & \displaystyle {\displaystyle \frac{1}{\mathit{Pr}\bar{\unicode[STIX]{x1D70C}}\bar{T}}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\bar{\unicode[STIX]{x1D70C}}\bar{T}\unicode[STIX]{x1D735}S+{\displaystyle \frac{c_{1}\mathit{Pr}}{R\bar{T}}}\left(Q_{\unicode[STIX]{x1D708}}+{\displaystyle \frac{1}{\mathit{Pm}\bar{\unicode[STIX]{x1D70C}}}}Q_{j}\right),\end{eqnarray}$$

(2.3e ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{B} & = & \displaystyle \unicode[STIX]{x1D735}\times (\boldsymbol{u}\times \boldsymbol{B})+\mathit{Pm}^{-1}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{B},\end{eqnarray}$$

$S$

(2.4) $$\begin{eqnarray}\overline{S}={\displaystyle \frac{\unicode[STIX]{x1D701}(r)^{-n}-\unicode[STIX]{x1D701}_{o}^{-n}}{\unicode[STIX]{x1D701}_{o}^{-n}-\unicode[STIX]{x1D701}_{i}^{-n}}},\end{eqnarray}$$

$\boldsymbol{u}$ is the velocity vector, $\boldsymbol{B}$ is the magnetic flux density, $\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}$ includes all terms that can be written as gradients and $\boldsymbol{r}=r\hat{\boldsymbol{r}}$ is the position vector with respect to the centre of the sphere (Jones et al. Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011; Simitev, Kosovichev & Busse Reference Simitev, Kosovichev and Busse2015). The viscous force, the viscous heating and the ohmic heating, given by

(2.5a-c ) $$\begin{eqnarray}\boldsymbol{F}_{\unicode[STIX]{x1D708}}=(\unicode[STIX]{x1D70C}_{c}/\bar{\unicode[STIX]{x1D70C}})\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64E},\quad Q_{\unicode[STIX]{x1D708}}=\unicode[STIX]{x1D64E}\boldsymbol{ : }\boldsymbol{e},\quad Q_{j}=(\unicode[STIX]{x1D735}\times \boldsymbol{B})^{2},\end{eqnarray}$$

respectively, are defined in terms of the deviatoric stress tensor

(2.6a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D61A}_{ij}=2\bar{\unicode[STIX]{x1D70C}}(e_{ij}-e_{kk}\unicode[STIX]{x1D6FF}_{ij}/3),\quad e_{ij}=(\unicode[STIX]{x2202}_{i}u_{j}+\unicode[STIX]{x2202}_{j}u_{i})/2,\end{eqnarray}$$

where double-dots (:) denote a Frobenius inner product, and $\unicode[STIX]{x1D708}$ is a constant viscosity. The governing equations (2.3) have been non-dimensionalised with the shell thickness $d=r_{o}-r_{i}$ as unit of length, $d^{2}/\unicode[STIX]{x1D708}$ as unit of time, $\unicode[STIX]{x1D708}\sqrt{\unicode[STIX]{x1D707}_{0}\unicode[STIX]{x1D70C}_{c}}/d$ as a unit of magnetic induction and $\unicode[STIX]{x0394}S$ , $\unicode[STIX]{x1D70C}_{c}$ and $T_{c}$ as units of entropy, density and temperature, respectively. The system is then characterised by eight dimensionless parameters: the radius ratio, the polytropic index of the gas, the density scale number, the Prandtl number, the magnetic Prandtl number, the Rayleigh number, the Ekman number and the baroclinicity parameter,

(2.7) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D702}=r_{i}/r_{o},\quad n,\quad N_{\unicode[STIX]{x1D70C}}=\ln (\bar{\unicode[STIX]{x1D70C}}(r_{i})/\bar{\unicode[STIX]{x1D70C}}(r_{o})),\quad \mathit{Pr}=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705},\quad \mathit{Pm}=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D706},\\ R=-gd^{3}\unicode[STIX]{x0394}S/(\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}c_{p}),\quad \mathit{Ek}=\unicode[STIX]{x1D708}/(\unicode[STIX]{x1D6FA}d^{2}),\quad \mathit{Z}=\unicode[STIX]{x1D6FA}^{2}d^{4}\unicode[STIX]{x0394}S/(\unicode[STIX]{x1D708}^{2}c_{p}),\end{array}\right\}\end{eqnarray}$$

respectively, where $\unicode[STIX]{x1D705}$ is a constant entropy diffusivity, $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D707}_{0}$ are the magnetic diffusivity and permeability, and $c_{p}$ is the specific heat at constant pressure. Note, that the Rayleigh number assumes negative values in the present problem since the basic entropy gradient is reversed with respect to the case of buoyancy-driven convection.

The poloidal–toroidal decomposition

(2.8) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\bar{\unicode[STIX]{x1D70C}}\boldsymbol{u}=\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D735}\times \hat{\boldsymbol{r}}rv)+\unicode[STIX]{x1D735}\times \hat{\boldsymbol{r}}r^{2}w,\\ \boldsymbol{B}=\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D735}\times \hat{\boldsymbol{r}}h)+\unicode[STIX]{x1D735}\times \hat{\boldsymbol{r}}g,\end{array}\right\}\end{eqnarray}$$

is used to enforce the solenoidality of the mass flux $\bar{\unicode[STIX]{x1D70C}}\boldsymbol{u}$ and of the magnetic flux density. This has the further advantages that the pressure gradient can be eliminated and scalar equations for the poloidal and the toroidal scalar fields, $v$ and $w$ , are obtained by taking $\hat{\boldsymbol{r}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D735}\times$ and $\hat{\boldsymbol{r}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times$ of equation (2.3c ). Similarly equations for the poloidal and the toroidal scalar fields, $h$ and $g$ , are obtained by taking $\hat{\boldsymbol{r}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times$ and $\hat{\boldsymbol{r}}\boldsymbol{\cdot }$ of equation (2.3e ). Except for the term with the baroclinicity parameter $\mathit{Z}$ the resulting equations are identical to those described by Simitev et al. (Reference Simitev, Kosovichev and Busse2015). Assuming that entropy fluctuations are damped by convection in the region above $r=r_{o}$ we choose the boundary condition

(2.9a ) $$\begin{eqnarray}S=0\quad \text{at }r=\left\{\begin{array}{@{}l@{}}r_{i}\equiv \unicode[STIX]{x1D702}/(1-\unicode[STIX]{x1D702}),\\ r_{o}\equiv 1/(1-\unicode[STIX]{x1D702}),\end{array}\right.\end{eqnarray}$$

while the inner and the outer boundaries of the shell are assumed stress free and impenetrable for the flow

(2.9b-d ) $$\begin{eqnarray}v=0,\quad \unicode[STIX]{x2202}_{r}^{2}v={\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}^{\prime }}{\bar{\unicode[STIX]{x1D70C}}r}}\unicode[STIX]{x2202}_{r}(rv),\quad \unicode[STIX]{x2202}_{r}(rw)={\displaystyle \frac{\bar{\unicode[STIX]{x1D70C}}^{\prime }}{\bar{\unicode[STIX]{x1D70C}}}}w\quad \text{at }r=\left\{\begin{array}{@{}l@{}}r_{i},\\ r_{o}.\end{array}\right.\end{eqnarray}$$

The boundary conditions for the magnetic field are derived from the assumption of an electrically insulating external region. The poloidal function $h$ is then matched to a function $h^{(e)}$ , which describes an external potential field,

(2.9e-g ) $$\begin{eqnarray}g=0,\quad h-h^{(e)}=0,\quad \unicode[STIX]{x2202}_{r}(h-h^{(e)})=0\quad \text{at }r=\left\{\begin{array}{@{}l@{}}r_{i},\\ r_{o}.\end{array}\right.\end{eqnarray}$$

3 Methods for numerical solution

A pseudo-spectral method described by Tilgner (Reference Tilgner1999) was employed for the direct numerical simulation of problem (2.3–2.8). We adapted a code developed and used by us for a number of years (Busse, Grote & Simitev Reference Busse, Grote, Simitev, Jones, Soward and Zhang2003; Busse & Simitev Reference Busse and Simitev2011; Simitev et al. Reference Simitev, Kosovichev and Busse2015) that has been extensively benchmarked for accuracy (Marti et al. Reference Marti, Schaeffer, Hollerbach, Cébron, Nore, Luddens, Guermond, Aubert, Takehiro and Sasaki2014; Simitev et al. Reference Simitev, Kosovichev and Busse2015; Matsui et al. Reference Matsui, Heien, Aubert, Aurnou, Avery, Brown, Buffett, Busse, Christensen and Davies2016). Adequate numerical resolution for the simulations was chosen as described in Simitev et al. (Reference Simitev, Kosovichev and Busse2015). To analyse the properties of the solutions we decompose the kinetic energy density into poloidal and toroidal components (respectively denoted by subscripts as in $X_{p}$ and $X_{t}$ in equations (3.1) below and where $X$ denotes an appropriate quantity) and further into mean (axisymmetric) and fluctuating (non-axisymmetric) components (respectively denoted by bars and tildes as in $\overline{X}$ and $\widetilde{X}$ below) and into equatorially symmetric and equatorially antisymmetric components (respectively denoted by superscripts as in $X^{s}$ and $X^{a}$ below),

(3.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\overline{E}_{p}=\overline{E}_{p}^{s}+\overline{E}_{p}^{a}=\langle (\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D735}(\overline{v}^{s}+\overline{v}^{a})\times \boldsymbol{r}))^{2}/(2\bar{\unicode[STIX]{x1D70C}})\rangle ,\\ \overline{E}_{t}=\overline{E}_{t}^{s}+\overline{E}_{t}^{a}=\langle (\unicode[STIX]{x1D735}r(\overline{w}^{s}+\overline{w}^{a})\times \boldsymbol{r})^{2}/(2\bar{\unicode[STIX]{x1D70C}})\rangle ,\\ \widetilde{E}_{p}=\widetilde{E}_{p}^{s}+\widetilde{E}_{p}^{a}=\langle (\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D735}(\widetilde{v}^{s}+\widetilde{v}^{a})\times \boldsymbol{r}))^{2}/(2\bar{\unicode[STIX]{x1D70C}})\rangle ,\\ \widetilde{E}_{t}=\widetilde{E}_{t}^{s}+\widetilde{E}_{t}^{a}=\langle (\unicode[STIX]{x1D735}r(\widetilde{w}^{s}+\widetilde{w}^{a})\times \boldsymbol{r})^{2}/(2\bar{\unicode[STIX]{x1D70C}})\rangle ,\end{array}\right\}\end{eqnarray}$$

where angular brackets $\langle \,\,\rangle$ denote averages over the volume of the spherical shell. Since in our code the spectral representation of all fields $X$ is given by the set of coefficients $\{X_{l}^{m}\}$ of their expansions in spherical harmonics $Y_{l}^{m}$ , it is easy to extract the relevant components, i.e. coefficients with $m=0$ and with $m\neq 0$ represent axisymmetric and non-axisymmetric components, respectively, while coefficients with even $(l+m)$ and with odd $(l+m)$ represent equatorially symmetric and equatorially antisymmetric components, respectively. The magnetic energy density is similarly decomposed into components.

4 Parameter values and initial conditions

In the simulations presented here fixed values are used for some governing parameters, namely

(4.1a-d ) $$\begin{eqnarray}\unicode[STIX]{x1D702}=0.3,\quad n=2,\quad N_{\unicode[STIX]{x1D70C}}=2,\quad R=-5\times 10^{4}.\end{eqnarray}$$

The value for the shell thickness represents a stellar radiative core geometrically similar to that of the Sun. The values for $n$ and $N_{\unicode[STIX]{x1D70C}}$ are not significantly different from current estimates for the Solar radiative zone, $n=1.5$ and $N_{\unicode[STIX]{x1D70C}}=4.6$ . The Rayleigh number is set to a negative value in order to model a convectively stable configuration.

A large variation is known to exist in the rates of stellar rotation, e.g. van Saders & Pinsonneault (Reference van Saders and Pinsonneault2013). Here, we use three different values of the Ekman number, namely $\mathit{Ek}=2/300$ , $1/300$ and $1/500$ , selected so that the effect of rotation is strong enough to govern the dynamics of the system, but not too strong to cause a significant increase in computational expense; similar values are used e.g. by Simitev et al. (Reference Simitev, Kosovichev and Busse2015) and in the cases F1–4 of Käpylä et al. (Reference Käpylä, Käpylä, Olspert, Warnecke and Brandenburg2017). This variation in the values of the Ekman number allows us to observe the effect of rotation.

Unresolved subgrid scales in convective envelopes are typically modelled by the assumption of approximately equal turbulent eddy diffusivities. As a consequence the choice of $\mathit{Pr}=1$ is often made in the modelling and simulation literature (Miesch et al. Reference Miesch, Matthaeus, Brandenburg, Petrosyan, Pouquet, Cambon, Jenko, Uzdensky, Stone and Tobias2015). However, estimates of the Prandtl number values based on molecular diffusivities are minute and it is thus unlikely that eddy diffusivities could increase the effective Prandtl number to a value of the order of unity. Furthermore, in a stably stratified system turbulence is expected to be anisotropic (Zahn Reference Zahn1992). With this in mind we have decreased the value of the Prandtl number from $0.1$ used in Simitev & Busse (Reference Simitev and Busse2017) to smaller values, namely $\mathit{Pr}=0.08$ , $0.05$ and $0.03$ . This variation in the value of $\mathit{Pr}$ suffices to produce pronounced differences in the simulation results. Further reductions of $\mathit{Pr}$ prove to be difficult numerically.

The molecular values of the magnetic Prandtl number are estimated to increase from approximately $10^{-5}$ at the surface of the Solar convection zone to approximately $10^{-1}$ at its bottom (Brandenburg & Subramanian Reference Brandenburg and Subramanian2005; Rieutord & Rincon Reference Rieutord and Rincon2010). Invoking the eddy diffusivity argument mentioned above, it may be expected that the effective values of the magnetic Prandtl number are somewhat further increased by turbulent mixing. Considering this, the value of $\mathit{Pm}=1.5$ for which a dynamo is demonstrated in the present paper is certainly large but, perhaps not excessively large. We remark, that it is possible to decrease $\mathit{Pm}$ further as discussed in relation to figure 10 later in the paper. With the aim of establishing the possibility of self-sustained dynamo action in all of the baroclinically driven flow states found (see below) we also report simulations with values of $\mathit{Pm}$ larger than 1.5. Finally, values similar to the ones used by us are also used by other numerical simulations reported in the literature e.g. many of the cases in Käpylä et al. (Reference Käpylä, Käpylä, Olspert, Warnecke and Brandenburg2017) have values significantly greater than unity even though theirs is a model of the Solar convection zone where $\mathit{Pm}$ is estimated to be smaller than in the radiative zone.

With these choices of parameter values, we find it convenient to organise and present our results in terms of several sequences of cases in which the parameter $\mathit{Z}$ is increased in small steps to study radiative zones of various degrees of baroclinicity while all other parameter values in a sequence are kept fixed, see for instance figure 1. We remark that the strength of the baroclinic forcing, measured by $\mathit{Z}$ , is limited from above so that

(4.2) $$\begin{eqnarray}\mathit{Z}<(1-\unicode[STIX]{x1D702})^{3}\;|R|/\mathit{Pr}.\end{eqnarray}$$

This restriction guarantees that the apparent gravity does not point outward such that the model explicitly excludes the well-studied case of buoyancy-driven thermal convection.

Initial conditions of no fluid motion are used at vanishingly small values of $\mathit{Z}$ , while at finite values of $\mathit{Z}$ the closest equilibrated neighbouring case is used as initial condition to help convergence and reduce transients. To ensure that transient effects are eliminated from the sequences presented below, all solutions have been continued for at least 15 time units. Similarly, several dynamo cases have been started with small random seeds for the magnetic field, while most other cases were subsequently started from neighbouring simulations with equilibrated dynamos.

Figure 1. Time-averaged kinetic energy densities as functions of baroclinicity $\mathit{Z}$ for parameter values (4.1) and $\mathit{Pr}=0.05$ , $\mathit{Ek}=2/300$ (a,b), $\mathit{Pr}=0.05$ , $\mathit{Ek}=1/300$ (c,d) and $\mathit{Pr}=0.03$ , $\mathit{Ek}=1/500$ (e,f). Full and empty symbols indicate equatorially symmetric and asymmetric energy components, respectively. Black circles, red squares, green triangles-up and blue triangles-down indicate the energy components $\overline{E}_{p}^{s,a}$ , $\overline{E}_{t}^{s,a}$ , $\widetilde{E}_{p}^{s,a}$ , $\widetilde{E}_{t}^{s,a}$ , respectively. Axially symmetric and axially asymmetric components are plotted in the left and the right panels, respectively. Vertical dash-dotted lines indicate transition points. The ranges over which distinct states are observed are indicated by arrows near the bottom abscissa, with some states co-existing as indicated. Energy components not shown are at least 10 orders of magnitude smaller than the ones shown.

5 Instabilities of baroclinically driven flows

The baroclinically driven problem is invariant under a group of symmetry operations including rotations about the polar axis i.e. invariance with respect to the coordinate transformation $\unicode[STIX]{x1D711}\rightarrow \unicode[STIX]{x1D711}+\unicode[STIX]{x1D6FC}$ , reflections in the equatorial plane $\unicode[STIX]{x1D703}\rightarrow \unicode[STIX]{x03C0}-\unicode[STIX]{x1D703}$ and translations in time $t\rightarrow t+a$ , see Gubbins & Zhang (Reference Gubbins and Zhang1993). As baroclinicity $\mathit{Z}$ is increased the available symmetries of the solution are broken resulting in a sequence of states ranging from simpler and more symmetric flow patterns to more complex flow patterns of lesser symmetry. In this respect the system resembles Rayleigh–Bénard convection (Busse Reference Busse2003).

Three sequences of non-magnetic simulations with gradually increasing values of the baroclinicity $\mathit{Z}$ and fixed values of the other parameters were obtained and are shown in figure 1. The first sequence is the most detailed one and allows us to capture the transitions that occur as the flow is more strongly driven while the other two sequences allow us to describe the effects of the variation in $\mathit{Ek}$ and $\mathit{Pr}$ .

The sequence with $\mathit{Pr}=0.05$ and $\mathit{Ek}=2/300$ is illustrated in the top two panels of figure 1 and in figure 2. The sequence starts with the basic axisymmetric, equatorially symmetric and time-independent state with a dominant wavenumber $k=1$ in the radial direction labelled A in figures 1 and 2. An instability occurs at about $\mathit{Z}=9\times 10^{4}$ leading to a state B characterised by a dominant azimuthal wavenumber $m=2$ in the expansion of the solution in spherical harmonics $Y_{l}^{m}$ , i.e. while the axisymmetry is broken, a symmetry holds with respect to the transformation $\unicode[STIX]{x1D711}\rightarrow \unicode[STIX]{x1D711}+\unicode[STIX]{x03C0}$ . Simultaneously, the symmetry about the equatorial plane is also broken in state B. At $\mathit{Z}=13\times 10^{4}$ a further transition to a pattern labelled C occurs. This state looks much like state A, i.e. it is axisymmetric, equatorially symmetric and time independent but differs from state A in that it keeps a dominant radial wavenumber $k=2$ . State C ceases to exists above $19\times 10^{4}$ . At $\mathit{Z}=22\times 10^{4}$ a further transition to a pattern labelled E occurs. State E seems to be related to state B in a way similar as state C is related to state B. Similarly to B, state E has a $m=2$ azimuthal symmetry and is asymmetric with respect to the equatorial plane. E differs from B, however, in that a dominant radial wavenumber $k=2$ develops. Remarkably, state C coexists with a pattern D that can be found for a range of values of $\mathit{Z}$ from $15\times 10^{4}$ to $22\times 10^{4}$ . Like C, state D is equatorially symmetric, has a dominant radial wavenumber $k=2$ and is time independent. Unlike C it is not axisymmetric, but exhibits a $m=2$ structure. Which one of the coexisting patterns C and D will be found in a given numerical simulation depends on the initial conditions.

Figure 2. Flow structures with increasing baroclinicity $\mathit{Z}\times 10^{-4}=1$ , $12$ , $18$ , $18$ and $30$ from top to bottom for $\mathit{Pr}=0.05$ , $\mathit{Ek}=2/300$ and fixed values (4.1). Bistability occurs at $\mathit{Z}\times 10^{-4}=18$ . The first plot in each row shows isocontours of $\overline{u}_{\unicode[STIX]{x1D711}}$ (left half) and streamlines $r\sin \unicode[STIX]{x1D703}(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\overline{v})=$ const. (right half) in the meridional plane. The second plot shows isocontours of $u_{r}$ at $r=r_{i}+0.7$ mapped to the spherical surface using isotropic Aitoff projection. The third plot shows isocontours of $u_{r}$ in the equatorial plane. The isocontours are equidistant with positive isocontours shown by solid lines, negative isocontours shown by broken lines and the zeroth isocontour shown by a dotted line in each plot. All contour plots are snapshots at a fixed representative moment in time. Letters at each row denote corresponding flow states as indicated in figure 1.

Figure 3. Flow structures with increasing baroclinicity $\mathit{Z}\times 10^{-4}=30$ , $32$ , $34$ for $\mathit{Pr}=0.05$ , $\mathit{Ek}=1/300$ and fixed values (4.1). The same quantities are plotted in each row as in figure 2 and letters on the left-hand side denote the corresponding flow state as in figure 1.

Figure 4. Flow in the case $\mathit{Z}\times 10^{-4}=120$ , for $\mathit{Pr}=0.03$ , $\mathit{Ek}=1/500$ and fixed values (4.1). The same quantities are plotted in each row as in figure 2 and letters denote the corresponding flow state as in figure 1.

Figure 5. Time series of kinetic energy densities in the case shown in figure 4. Equatorially symmetric components are shown in (a) and equatorially antisymmetric components are shown in (b). The second panel in each group shows an enlarged segment of the time axis. The components $\overline{E}_{p}$ , $\overline{E}_{t}$ , $\widetilde{E}_{p}$ , $\widetilde{E}_{t}$ , are labelled in the plot and are also indicated by black lines, red lines, green lines and blue lines, respectively.

The Ekman number is decreased to $\mathit{Ek}=1/300$ in our second sequence while we keep $\mathit{Pr}=0.05$ . The sequence is illustrated in the middle two panels of figure 1 and in figure 3. States C and E are not observed. State B becomes rather more pronounced with patches detaching from each other and with arms shooting towards the poles. State D loses its twofold azimuthal symmetry $m=2$ . In addition, state D becomes time periodic at the largest values of $\mathit{Z}$ , e.g. at $\mathit{Z}=34\times 10^{4}$ shown in the last row of figure 3.

The Prandtl number is decreased to $\mathit{Pr}=0.03$ and the Ekman number is further decreased to $\mathit{Ek}=1/500$ in our third sequence. The kinetic energy components are plotted in the bottom two panels of figure 1. In comparison with the preceding sequence, state D is not observed. For values of the baroclinicity larger than $\mathit{Z}=110\times 10^{4}$ the spatial structure of state B becomes irregular as shown in figure 4 and exhibits a chaotic time dependence as illustrated by the time series of its kinetic energy components in figure 5.

A summary of the symmetry properties of all distinct flow states identified is listed in table 1 as is their capability of self-sustained dynamo action, see next section.

We wish to close this section by a brief comparison of baroclinically driven flows with the extensively studied problem of convective flows in rotating systems. Typical convection-driven finite-amplitude instabilities are described for instance by Sun, Schubert & Glatzmaier (Reference Sun, Schubert and Glatzmaier1993), Simitev & Busse (Reference Simitev and Busse2003), Busse & Simitev (Reference Busse and Simitev2005b ) and it is apparent that the baroclinically driven flows described in the article are essentially different. One feature that appears similar at first sight is the finding that baroclinically driven flows similarly to buoyancy-dominated convective turbulent flows (or equivalently slowly rotating convection) are characterised by retrograde (anti-Solar) differential rotation. However, in rotating turbulent convection anti-Solar rotation occurs after a sharp transition from prograde (Solar) rotation as the buoyancy is gradually increased (Simitev et al. Reference Simitev, Kosovichev and Busse2015), see also Gastine et al. (Reference Gastine, Yadav, Morin, Reiners and Wicht2013), Käpylä, Käpylä & Brandenburg (Reference Käpylä, Käpylä and Brandenburg2014). Moreover, convection in the Solar rotation regime is characterised by elongated convective columns, known as ‘Busse columns’, while convection in the anti-Solar regime is strongly disorganised, see Simitev et al. (Reference Simitev, Kosovichev and Busse2015). Furthermore, anti-Solar rotation is reversed to Solar rotation in the presence of magnetic field (Simitev et al. Reference Simitev, Kosovichev and Busse2015) and other effects of the magnetic field on this transition are reported in Fan & Fang (Reference Fan and Fang2014), Karak et al. (Reference Karak, Käpylä, Käpylä, Brandenburg, Olspert and Pelt2015). In contrast, in baroclinically driven flows anti-Solar rotation occurs at all values of the baroclinicity starting from the basic flow state A, the baroclinic flow is regular and well organised, and the magnetic fields generated by baroclinic dynamos do not reverse anti-Solar rotation as seen in the figures of § 6 below.

Table 1. Summary of symmetry properties and dynamo capability of states in the case $\mathit{Pr}=0.05$ , $\mathit{Ek}=2/300$ and fixed values (4.1).

6 Dynamos generated by baroclinically driven flows

The possibility of dynamos generated in stably stratified stellar radiative regions has received only limited support in the literature, e.g. Braithwaite (Reference Braithwaite2006), because it is well known that dynamo action does not exist in a spherical system in the absence of a radial component of motion (Bullard & Gellman Reference Bullard and Gellman1954; Kaiser & Busse Reference Kaiser and Busse2017). The latter is, indeed, rather weak in the basic state A of low poloidal kinetic energy as discussed above. However, with increasing baroclinicity $\mathit{Z}$ , the growing radial component of the velocity field strengthens the likelihood of dynamo action. Indeed, the existence of dipolar dynamos sustained by the equatorially symmetric flow state D was demonstrated in Simitev & Busse (Reference Simitev and Busse2017) for $\mathit{Pr}=0.1$ , $\mathit{Ek}=2/300$ and $\mathit{Pm}=16$ . Below we investigate further the existence of dynamo action in all other baroclinically driven flow states described in the preceding section.

A quadrupolar dynamo solution sustained by a baroclinically driven flow in state B at $\mathit{Pr}=0.03$ , $\mathit{Ek}=1/500$ and $\mathit{Z}=100\times 10^{4}$ and for a magnetic Prandtl number value $\mathit{Pm}=4$ is presented in figure 6. This is a steady dynamo. The ratio of poloidal to toroidal magnetic energy is $E_{p}^{\text{magn}}/E_{t}^{\text{magn}}=0.033$ , the ratio of $E_{\text{dipole}}^{\text{magn}}/E_{\text{quadrupole}}^{\text{magn}}=0.596$ and the ratio of magnetic to kinetic total energy is $E^{\text{magn}}/E^{\text{kin}}=0.009$ . The magnetic field has a weaker dipolar component and a large-scale quadrupolar topology that does not change in time. The surface structure of the magnetic field is characterised by relatively strong inward magnetic flux at both poles. Four localised patches of inward magnetic flux appear in the equatorial region. The azimuthally averaged toroidal field shows two large-scale strong flux tubes near the pole. The azimuthally averaged poloidal field is rather remarkable in that it remains almost completely confined to the spherical shell giving rise to an ‘invisible dynamo’.

This invisibility, however, does not persist as illustrated by a more strongly driven quadrupolar dynamo at $\mathit{Z}=110\times 10^{4}$ and $\mathit{Pm}=2$ shown in figure 7. The dynamo exhibits regular oscillations in which the quadrupole component reverses polarity.

Figure 6. A dynamo solution for a flow in state B in the case $\mathit{Z}\times 10^{-4}=100$ , with $\mathit{Pr}=0.03$ , $\mathit{Ek}=1/500$ , $\mathit{Pm}=4$ and values (4.1). (a) The first plot shows isocontours of radial magnetic field $B_{r}$ at $r=r_{o}+0.1$ in isotropic Aitoff projection. The second plot shows contours of $B_{r}$ in the equatorial plane. The third plot shows isocontours of $\overline{B}_{\unicode[STIX]{x1D711}}$ (left half) and meridional field lines $r\sin \unicode[STIX]{x1D703}\,\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\overline{h}=$ const. (right half). (b) The first plot shows isocontours of radial velocity $u_{r}$ at $r=r_{i}+0.7$ in isotropic Aitoff projection. The second plot shows contours of $u_{r}$ in the equatorial plane. The third plot shows isocontours of $\overline{u}_{\unicode[STIX]{x1D711}}$ (left half) and streamlines $r\sin \unicode[STIX]{x1D703}(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\overline{v})=$ const. (right half) in the meridional plane. All contour plots are snapshots at a fixed representative moment in time. The letter B denotes the corresponding flow state as indicated in figure 1.

Figure 7. A dynamo solution for a flow in state B in the case $\mathit{Z}\times 10^{-4}=110$ , $\mathit{Pr}=0.03$ , $\mathit{Ek}=1/500$ , $\mathit{Pm}=2$ and values (4.1). The same quantities are plotted as in figure 6 and the letter B denotes the corresponding flow state as indicated in figure 1.

Figure 8. A dynamo solution for a flow in state D in the case $\mathit{Z}\times 10^{-4}=32$ , $\mathit{Pr}=0.05$ , $\mathit{Ek}=1/300$ , and $\mathit{Pm}=5$ and values (4.1). The same quantities are plotted as in figure 6 and the letter D denotes the corresponding flow state as indicated in figure 1.

Figure 9. A dynamo solution for a flow in state E in the case $\mathit{Z}\times 10^{-4}=34$ , $\mathit{Pr}=0.05$ , $\mathit{Ek}=2/300$ , and $\mathit{Pm}=1.5$ and values (4.1). The same quantities are plotted as in figure 6 and the letter E denotes the corresponding flow state as indicated in figure 1.

A multipolar dynamo solution sustained by a baroclinically driven flow in state D at $\mathit{Pr}=0.05$ , $\mathit{Ek}=1/300$ and $\mathit{Z}=32\times 10^{4}$ and $\mathit{Pm}=5$ is presented in figure 8. This dynamo exhibits a steady time dependence. The ratio of poloidal to toroidal magnetic energy is $E_{p}^{\text{magn}}/E_{t}^{\text{magn}}=0.384$ , the ratio of $E_{\text{dipole}}^{\text{magn}}/E_{\text{quadrupole}}^{\text{magn}}=0.0004$ and the ratio of magnetic to kinetic total energy is $E^{\text{magn}}/E^{\text{kin}}=0.003$ . The spatial structure of the magnetic field of this dynamo is highly unusual in that it is not dominated by a large-scale dipolar component ( $Y_{1}^{0}$ ) or quadrupolar component ( $Y_{2}^{0}$ ) as is the case with most other large-scale dynamos reported to date but rather resembles approximately the structure of the $Y_{6}^{4}$ spherical harmonic function as seen in the leftmost plot in the upper row of figure 8. Because of the weak $m=0$ contribution to the magnetic field the azimuthally averaged components plotted in the rightmost plot in the upper row of figure 8 are insignificant. It is remarkable that in this case the magnetic field significantly alters the fluid flow pattern as evident by comparison with the non-magnetic simulation at identical parameter values shown in the second row of figure 3, that was used as an initial condition for the dynamo run. The new magnetic flow pattern retains profiles of the differential rotation and the meridional circulation similar to those of the non-magnetic case. However, the dominant azimuthal wavenumber increases from 2 to 4 with a stronger contribution of $m=0$ . None of the other dynamos we report alter their generating flow states so significantly.

A dipolar dynamo generated by a baroclinically driven flow in state E for $\mathit{Pr}=0.05$ , $\mathit{Ek}=2/300$ , $\mathit{Z}=34\times 10^{4}$ and $\mathit{Pm}=1.5$ is shown in figure 9. The case has been started from an equilibrated neighbouring case to help convergence and reduce transients. The dynamo solution is steady. The ratio of poloidal to toroidal magnetic energy is $E_{p}^{\text{magn}}/E_{t}^{\text{magn}}=0.078$ , the ratio of $E_{\text{dipole}}^{\text{magn}}/E_{\text{quadrupole}}^{\text{magn}}=5.97$ and the ratio of magnetic to kinetic total energy is $E^{\text{magn}}/E^{\text{kin}}=0.344$ . Furthermore, the energy density of the magnetic field $E^{\text{magn}}$ is comparable to the kinetic energy density of the poloidal component of the velocity field, $E_{p}^{\text{kin}}$ , with a ratio $E^{\text{magn}}/E_{p}^{\text{kin}}=1.42$ . The magnetic field has a weaker quadrupolar component and a large-scale dipolar topology with prominent polar magnetic fluxes that does not change in time, as shown in figure 9, resembling in this respect the surface fields of Ap–Bp stars (Donati & Landstreet Reference Donati and Landstreet2009). The azimuthally averaged toroidal field shows a pair of hook-shaped toroidal flux tubes largely filling each hemisphere of the shell. A large-scale dipolar component emerges outside of the spherical shell. The structure of the fluid flow remains largely unaffected by the presence of the magnetic field.

Figure 10. (a) Time-averaged total magnetic energy density as a function of the magnetic Prandtl number for $\mathit{Ek}=2/300$ in the cases $\mathit{Pr}=0.1$ , $\mathit{Z}=17\times 10^{4}$ (blue triangles), $Pr=0.08$ , $\mathit{Z}=20\times 10^{4}$ (black circles), $P=0.05$ , $\mathit{Z}=32\times 10^{4}$ (red squares). Small black rectangles indicate the boundary of a region of no dynamo excitation for the respective sequence. (b) Components of the magnetic and kinetic energy densities (as labelled) as a function of time in the case $Pr=0.08$ , $\mathit{Z}=20\times 10^{4}$ , $\mathit{Ek}=2/300$ , $\mathit{Pm}=12$ .

An important issue for application of these results to astrophysical objects is to determine whether dynamo action persists at sufficiently small values of the magnetic Prandtl number. Figure 10(a) shows the dependence of the time-averaged total magnetic energy density on the magnetic Prandtl number for three sequences of self-sustained dynamo cases with $\mathit{Ek}=2/300$ . In each sequence no active dynamo is found for smaller values of $\mathit{Pm}$ ; thus the smallest value of $\mathit{Pm}$ is an estimate of the critical value of $\mathit{Pm}$ for dynamo onset at fixed values of the other parameters. In particular, the figure shows that $\mathit{Pm}_{\text{crit}}$ decreases as the value of the ordinary Prandtl number $\mathit{Pr}$ is decreased and the value of the baroclinicity parameter $\mathit{Z}$ is increased reducing $\mathit{Pm}_{\text{crit}}$ to $1.5$ at $P=0.05$ and $\mathit{Z}=32\times 10^{4}$ . While it is numerically challenging to proceed in the direction of stronger baroclinic driving and smaller ordinary Prandtl number this trend, which we believe is likely to continue, indicates that baroclinically driven dynamos at more realistic values of $\mathit{Pm}$ eventually exist. In this respect baroclinically driven dynamos seem to behave similarly to both Taylor–Couette dynamos (Willis & Barenghi Reference Willis and Barenghi2002) and to randomly forced small-scale dynamos (Schekochihin et al. Reference Schekochihin, Cowley, Maron and McWilliams2004, Reference Schekochihin, Haugen, Brandenburg, Cowley, Maron and McWilliams2005) where evidence is found that for fixed other parameter values a critical value of $Pm$ exists below which dynamo action is not possible. Surprisingly, figure 10(a) also shows that there is an optimal value of $\mathit{Pm}$ for magnetic field generation, e.g.  $\mathit{Pm}_{\text{opt}}=5$ at $\mathit{Pr}=0.08$ , $\mathit{Z}=20\times 10^{4}$ and that dynamo action is lost when sufficiently large values of $\mathit{Pm}$ are used. This appears to be due to a decline of the non-axisymmetric poloidal and toroidal kinetic energy components, $\widetilde{E}_{p}$ and $\widetilde{E}_{t}$ , which are most strongly depleted by the magnetic field. This is related to the specific mechanism of dynamo excitation as discussed below.

A common feature of all dynamo-capable baroclinic flow states is that they exhibit significant non-axisymmetric flow components, $\widetilde{E}_{p}$ and $\widetilde{E}_{t}$ as seen in the left panels of figure 1. These non-axisymmetric flow components appear to be essential for dynamo action as strikingly illustrated in figure 10(b). Here a solution for $\mathit{Pr}=0.08$ , $\mathit{Z}=20\times 10^{4}$ , $\mathit{Ek}=2/300$ and for $\mathit{Pm}=12$ is shown. In the initial part of the simulation an oscillatory dipolar dynamo field is sustained by a flow in state D which is characterised by significant non-axisymmetric flow components. Flow state D, however, co-exists with flow state C which has negligibly small non-axisymmetric flow components and does not seem to support dynamo action. As $\widetilde{E}_{p}$ and $\widetilde{E}_{t}$ are weakened by the magnetic field a flow transition from state D to state C occurs after which the dynamo rapidly decays. On their own, the mean components of the flow which in the basic state A consist of radially decreasing differential rotation (‘subrotation’) coupled with counterclockwise meridional circulation in the northern hemisphere do not lead to dynamo excitation as also shown in kinematic dynamo models with prescribed laminar flows by Dudley & James (Reference Dudley and James1989), see also Latter & Ivers (Reference Latter and Ivers2010). However, since the differential rotation remains the dominant flow component in all states the baroclinically driven dynamos reported here are of $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FA}$ type.