Differential Rotation in Convecting Spherical Shells with Non-Uniform Viscosity and Entropy Diffusivity

Abstract

Contemporary three-dimensional physics-based simulations of the solar convection zone disagree with observations. They feature differential rotation substantially different from the true rotation inferred by solar helioseismology and exhibit a conveyor belt of convective “Busse” columns not found in observations. To help unravel this so-called “convection conundrum”, we use a three-dimensional pseudospectral simulation code to investigate how radially non-uniform viscosity and entropy diffusivity affect differential rotation and convective flow patterns in density-stratified rotating spherical fluid shells. We find that radial non-uniformity in fluid properties enhances polar convection, which, in turn, induces non-negligible lateral entropy gradients that lead to large deviations from differential rotation geostrophy due to thermal wind balance. We report simulations wherein this mechanism maintains differential rotation patterns very similar to the true solar profile outside the tangent cylinder, although discrepancies remain at high latitudes. This is significant because differential rotation plays a key role in sustaining solar-like cyclic dipolar dynamos.

1. Introduction

Solar magnetic phenomena originate in the dynamic processes of thermal convection of electrically conducting plasma deep within the solar interior. High-resolution numerical simulations of thermal convection and magnetic field generation in rotating spherical shells based on fundamental physical laws have emerged as useful tools for understanding the complexities of solar dynamics [1,2]. Global simulations are increasingly compared to observational data and used to aid in observational interpretation. For instance, recent works have focused on leveraging time–distance helioseismology to constrain models of meridional circulation [3], while other studies have explored minimal models to explain solar cycle dynamics by emphasizing the significance of non-axisymmetric ($m=1$ or $m=2$) components of the magnetic field [4].

Despite the many advances made using numerical simulations, limitations and challenges persist in capturing several fundamental aspects of the solar dynamo. The primary issue we investigate in this article is these challenges, i.e., the so-called convective conundrum [5,6], which manifests as (a) a significant disparity between observed and simulated differential rotation and convective velocities and (b) the prominence of a conveyor belt of large-scale convective columns (Busse columns) in simulations with little corresponding observational evidence [2,7]. Differential rotation, a crucial part of the Sun’s global behaviour, is very well measured by helioseismic inversion [8,9]. Its primary features are the variation in rotation speed with latitude, with the equator rotating faster (around 25 days) and the poles rotating slower (around 36 days), displaying a conical profile with a nearly radial orientation of angular velocity isocontours at medium latitudes.

Simulations generally fail to produce isocontours of the angular velocity that are significantly inclined from the axis of rotation in the convection zone, and most simulations feature unrealistic isocontours parallel to the the axis of rotation instead [10,11,12]. Notable exceptions are the simulations reported in [13], where conically tilled contours were obtained by imposing a moderately large latitudinal entropy gradient. It is especially important to capture differential rotation correctly, as it plays a crucial role in magnetic field generation via the $\mathrm{\Omega }$ effect of poloidal-to-toroidal field conversion, in addition to inducing dynamo wave oscillations akin to the solar cycle, as demonstrated by flux-transport dynamo models [14] and some global simulations [4,15,16,17]. In the following sections, we demonstrate that a significantly large latitudinal entropy gradient can arise self-consistently without being additionally imposed, helping to improve agreement between our simulations and observations.

In contrast to planetary cores and atmospheres, where basic density and material properties of the fluid change only weakly, in the solar interior, including the convection zone, there are very significant radial variations of thermodynamic quantities and properties. Reference-state radial distributions of density, temperature and pressure in the convection zone can be estimated based on helioseismic arguments or solar evolution models (e.g., [18]). In theory, molecular viscosity and thermal diffusivity profiles can be inferred from molecular dynamics calculations (e.g., [19]), but because of the strong turbulence in the solar convection zone, it is unlikely that these molecular values adequately represent the profiles of the effective viscosity and thermal diffusivity in global convective dynamo models. Hence, the radial distributions of viscosity and entropy diffusivity remain largely modelling choices—although rather important ones, as they are expected to influence the style and the spatial location of convection and, by extension, the properties of the global dynamo process. For example, the radial profile of entropy diffusivity directly affects the entropy distribution, whose radial gradient consequently determines the local convective stability. Various choices of viscosity and entropy diffusivity profiles were investigated in the early study conducted by Glazmaier and Gilman [20], who found that under a moderate shell rotation rate, convection moves from the outer to inner regions as the diffusivities are increased in the outer regions and decreased in the inner regions. Several studies of solar convection have assumed that viscosity and diffusivity are functions involving a negative power of the mean density (e.g., [18,20,21,22]) as opposed to the more commonly used uniform profiles [12,16,17,23,24,25,26]. In a recent study conducted by Sasaki et al. [27], the effects of a radial distribution of entropy diffusivity on critical modes of anelastic thermal convection in rotating spherical shell were studied.They found that convection morphology and location are strongly affected, but they restricted their analysis to the linear onset and considered only uniform viscosity. There is a need to extend the work reported in [27] to non-linear regimes, especially with an eye to improving the agreement between computed and observed solar differential rotation. This is the goal of the present study.

To this end, we report a set of numerical simulations based on the model proposed in [17], which we extended to incorporate radially non-uniform profiles of viscosity and entropy diffusivity. To pinpoint the effect of these assumptions, we performed computations and compared the results with those of baseline reference simulations in which uniform profiles were used. We consider both high and low values of the Prandtl number to make sure our findings are robust. The structure of this paper is outlined as follows. In Section 2, we introduce the mathematical model used in this study. Section 3 describes our main results with respect to the effects of non-uniform material properties on the differential rotation structure. In Section 4, we propose a likely mechanism to help explain our findings, provide evidence in the form of additional results and mention limitations of the study. We finish with a brief conclusion in Section 5.

2. Mathematical Model

In this work, we extend the model of Jones at al. [28] by including radially non-uniform viscosity and entropy diffusivity. Our numerical implementation, methods of solution and analysis follow [17]. Details are recapitulated below for completeness.

We consider an electrically conducting perfect gas confined to a spherical shell. The shell rotates about the vertical axis with a fixed angular velocity (${\mathrm{\Omega }}_0\widehat{\mathit{k}}$), and an entropy contrast ($\mathrm{\Delta }S$) is imposed between its inner and outer surfaces. Assuming a gravitational field the inverse square of the radial distance from the centre ($r^{-2}$), a hydrostatic, polytropic reference state exists

$$\begin{array}{c}\overline{\rho }={\rho }_c{\zeta }^n,\overline{T}=T_c\zeta ,\overline{P}=P_c{\zeta }^{n+1},\end{array}$$

(1)

where $\zeta =c_0+c_{1}d/r$ is a radial profile with parameters $c_0=(2{\zeta }_o-\eta -1)/(1-\eta )$, $c_1=(1+\eta )(1-{\zeta }_o)/{(1-\eta )}^2$ and ${\zeta }_o=(\eta +1)/(\eta exp(N_{\rho }/n)+1)$. Constants ${\rho }_c$, $P_c$ and $T_c$ are reference values of density, pressure and temperature, respectively, in the middle of the shell. The gas polytropic index (n), the density scale height number ($N_{\rho }$) and the shell thickness ratio ($\eta$) are defined below.

Convection and magnetic field generation in this system are described by the evolution equations of continuity, momentum, energy and magnetic flux. Using the Lantz–Braginsky formulation of the anelastic approximation [28], these can be written in the following form:

$$\begin{array}{c}\nabla ·\overline{\rho }\mathbf{u}=0,\nabla ·\mathit{B}=0,\end{array}$$

(2a)

$$\begin{array}{c}{\partial }_t\mathbf{u}+(\nabla \times \mathbf{u})\times \mathbf{u}=-\nabla \mathrm{\Pi }-\tau (\widehat{\mathit{k}}\times \mathbf{u})+\frac{\mathrm{R}}{\mathrm{Pr}}\frac{S}{r^2}\widehat{\mathit{r}}+\frac{{\rho }_c}{\overline{\rho }}\nabla ·\widehat{\mathit{\sigma }}+\frac{1}{\overline{\rho }}(\nabla \times \mathit{B})\times \mathit{B},\end{array}$$

(2b)

$$\begin{array}{c}{\partial }_{t}S+\mathbf{u}·\nabla S=\frac{1}{\mathrm{Pr}\overline{\rho }\overline{T}}\nabla ·\overline{\kappa }\overline{\rho }\overline{T}\nabla S+\frac{c_1\mathrm{Pr}}{\mathrm{R}\overline{T}}\left(\widehat{\mathit{\sigma }}:e+\frac{1}{\mathrm{Pm}\overline{\rho }}{(\nabla \times \mathit{B})}^2\right),\end{array}$$

(2c)

$$\begin{array}{c}{\partial }_t\mathit{B}=\nabla \times (\mathbf{u}\times \mathit{B})+{\mathrm{Pm}}^{-1}{\nabla }^2\mathit{B},\end{array}$$

(2d)

where $\mathbf{u}$ is the velocity, B is the magnetic flux density, S is the entropy, $\nabla \mathrm{\Pi }$ includes all terms that can be written as gradients and $\widehat{\mathit{\sigma }}$ is the deviatoric stress tensor:

$$\begin{array}{c}{\sigma }_{ij}=2\overline{\nu }\overline{\rho }(e_{ij}-e_{kk}{\delta }_{ij}/3),e_{ij}=({\partial }_i{u}_j+{\partial }_j{u}_i)/2,\end{array}$$

where the double-dot symbol (:) denotes a component-wise inner product.

The governing equations are non-dimensionalized using the shell thickness ($d=r_o-r_i$) as a unit of length, $d^2/{\nu }_c$ as a unit of time, $\mathrm{\Delta }S$ as a unit of entropy, ${\nu }_c\sqrt{{\mu }_0{\rho }_c}/d$ as a unit of magnetic induction, ${\rho }_c$ as a unit of density and $T_c$ as a unit of temperature. Here, $r_i$ and $r_o$ are the inner and the outer radius, respectively; $\lambda$ and ${\mu }_0$ are the magnetic diffusivity and permeability, respectively; and $\nu$ and $\kappa$ are the viscosity and the entropy diffusivity, respectively. The latter pair is assumed to have radially non-uniform profiles:

$$\begin{array}{c}\nu \left(r\right)={\nu }_c{\left(\frac{\overline{\rho }}{{\rho }_c}\right)}^p,\kappa \left(r\right)={\kappa }_c{\left(\frac{\overline{\rho }}{{\rho }_c}\right)}^q,\end{array}$$

(3)

where ${\nu }_c$ and ${\kappa }_c$ are their values in the mid-shell, while p and q are real modelling constants. The formulation is then characterized by seven non-dimensional parameters: the radius ratio, the polytropic index, the density scale number, the Rayleigh number, the thermal Prandtl number, the magnetic Prandtl number and the Coriolis number, respectively defined as

$$\begin{array}{c}\eta :=\frac{r_i}{r_o},n,N_{\rho }:=ln\frac{\overline{\rho }\left(r_i\right)}{\overline{\rho }\left(r_o\right)},\mathrm{R}:=\frac{c_1{T}_c{d}^2\mathrm{\Delta }S}{{\nu }_c{\kappa }_c},\mathrm{Pr}:=\frac{{\nu }_c}{{\kappa }_c},\mathrm{Pm}:=\frac{{\nu }_c}{\lambda },\tau :=\frac{2{\mathrm{\Omega }}_0{d}^2}{{\nu }_c}.\end{array}$$

(4)

Equations (2a)–(2d) are supplemented by no-slip boundary conditions on the inner surface, with stress-free conditions on the outer surface:

$$\begin{array}{c}v=0,{\partial }_{r}v=0,w=0,\mathrm{at}r=r_i,\end{array}$$

(5a)

$$\begin{array}{c}v=0,{\partial }_r^{2}v-\frac{{\overline{\rho }}^{\prime }}{\overline{\rho }r}{\partial }_r\left(rv\right)=0,{\partial }_{r}w-\frac{{\overline{\rho }}^{\prime }}{\overline{\rho }}w=0\mathrm{at}r=r_o.\end{array}$$

(5b)

For the entropy, a fixed contrast is imposed between the inner and outer surfaces:

$$\begin{array}{c}S=1\mathrm{at}r=r_i,S=0\mathrm{at}r=r_o.\end{array}$$

(5c)

Lastly, “vacuum” boundary conditions for the magnetic field are derived from the assumption of an electrically insulating external region.

$$\begin{array}{c}g=0,h-h^{\left(e\right)}=0,{\partial }_r(h-h^{\left(e\right)})=0\mathrm{at}r=r_i,r_o.\end{array}$$

(5d)

In the above, v, w, h and g are poloidal and toroidal scalar fields in the following decompositions of of the momentum and magnetic fields, respectively:

$$\begin{array}{c}\overline{\rho }\mathit{u}=\nabla \times (\nabla \times \widehat{\mathit{r}}rv)+\nabla \times \widehat{\mathit{r}}{r}^{2}w,\end{array}$$

(6a)

$$\begin{array}{c}\mathit{B}=\nabla \times (\nabla \times \widehat{\mathit{r}}h)+\nabla \times \widehat{\mathit{r}}g,\end{array}$$

(6b)

The latter holds, since the mass flux ($\overline{\rho }\mathbf{u}$) and the magnetic flux density (B) are solenoidal vector fields. $\widehat{\mathit{r}}$ is the radial unit vector.

The boundary value problem expressed by Equations (2a)–(2d) and conditions (5a)–(5d) is three-dimensional, time-dependent, highly coupled and non-linear; for these reasons, it can only be solved numerically. Our method of numerical solution further exploits the poloidal–toroidal decomposition (6a) and (6b) as follows. Equation (2a) is satisfied automatically. Upon the application of $\widehat{\mathit{r}}·\nabla \times \nabla \times$ and $\widehat{\mathit{r}}·\nabla \times$ to Equation (2b), we obtain scalar equations for v and w and eliminate the pressure gradient. Similarly, applying $\widehat{\mathit{r}}·\nabla \times$ and $\widehat{\mathit{r}}·$ to Equation (2d), we find equations for h and g. The scalar unknowns (v, w, h, g and S) are then expanded in Chebychev polynomials in the radial variable (r) and in spherical harmonics in the angular variables $(\theta ,\phi )$. Their expansion coefficients are determined using a pseudospectral method adapted from [29]. The calculations reported below were performed with a truncations of 71, 193 and 193 for radial, zonal and latitudinal expansion coefficients, respectively.

3. Results

3.1. Design of the Study and Choice of Parameter Values

An evolution equation for the differential rotation generated by thermal convection in a rotating spherical shell can be obtained by averaging the zonal component of the momentum Equation (2b) over longitude. The result reveals that differential rotation is driven and maintained by a balance of Reynolds and Maxwell stresses modulated by stresses due to meridional circulation, magnetic tension and viscous dissipation [1,30]. Naturally, parameter values affect all of these stresses. Aiming to isolate the effects of radially non-uniform viscosity and entropy diffusivity profiles on the overall shape and structure of the differential rotation, we keep as many of the non-dimensional parameters as possible at fixed values. In particular, we set

$$\begin{array}{c}\eta =0.65,n=2,N_{\rho }=3,\mathrm{R}=3\times {10}^6,\tau =2000.\end{array}$$

(7)

Fixing parameter values is also helpful in restricting the large eight-dimensional parameter space of the problem to a manageable size. At $\eta =0.65$, the shell is slightly thicker than the solar convection zone ($\eta =0.7$), but this choice is made for ease of numerical simulation. The typical size of convective structures also depends on the thickness of the shell; thinner shells require spherical harmonic expansions of a higher order and degree to resolve the angular structure of flows. At $\tau =2000$, the Coriolis number is moderately but not excessively large, reflecting the model assumption that the flow in the deep convection zone is buoyancy- rather than rotation-dominated, although the effects of rotation are also essential. The value of the polytropic index ($n=2$) is larger than that for ideal gas ($n=3/2$), and the value of the density-scale height ($N_{\rho }=3$) is much smaller than that estimated for the solar convection zone. However, increasing $N_{\rho }$ much beyond 5 becomes computationally unfeasible. These choices are not expected to affect the dynamics significantly. They are also made for consistency and comparison with our earlier work [17].

The structure and intensity of the flow are most sensitive to the value of the Rayleigh number, which has been extensively studied (see, e.g., [17,31] and references therein). Here, we report a comparison of a few selected typical cases, which have an identical and moderately large value of the Rayleigh number. In Section 4.2, we discuss an extended sequence of simulations in which the latter parameter is increased systematically from the onset. While the Rayleigh number dependence is relatively well understood, the dependence on the Prandtl number is much less so, although we refer to [32,33]. Linear stability analyses [34,35], as well as finite-amplitude simulations of thermal convection in a rapidly rotating thick shell geometry under Boussinesq approximation [31,32], show that at significantly low Prandtl number values, convection develops in the so-called equatorially attached regime [35,36] for which the Coriolis force does not play a strong role in the force balance. At Prandtl number values close to unity, the balance starts to include Coriolis effects, and the regime of spiraling convection is entered, where convective structures become elongated, protruding from the inner to the outer surface with strong spiralling. Thus, in order to capture, to some extent, the dependence on the Prandtl number, we contrast simulations conducted with small ($\mathrm{Pr}=0.3$), moderate ($\mathrm{Pr}=1$) and somewhat larger values ($\mathrm{Pr}=5$).

To address the objectives of this study, we then compare a set of simulations performed using uniform radial profiles of viscosity and entropy diffusivity to a set of simulations in which these two profiles vary in the radial direction as $\overline{\nu }\propto {\rho }^{-0.9}$ and $\overline{\kappa }\propto {\rho }^{-0.5}$, respectively. In the absence of stringent observational or accurate theoretical constraints, these dependencies were selected so as to maximize the deviation of the two profiles from the uniform profiles at the fixed parameter values mentioned aboveand to ensure that well-resolved solutions are obtained. The two profiles, along with the radial profile of density in the model, are illustrated in Figure 1a.

The values of the Prandtl, Rayleigh and Coriolis numbers discussed above correspond to the fluid properties in the middle of the spherical shell (see definition (4)). Because of the non-uniform viscosity and the entropy diffusivity used in this study, the effective values of these quantities also vary in terms of radius, as illustrated in Figure 1b. This variation significantly affects the local morphology and intensity of the flows, as discussed further below.

In summary of this discussion, in

Table 1. Summary of model parameter values and computed energy densities for the six selected dynamo solutions discussed in the text.

	A	B	C	D	E	F
$\mathbf{Type}$	Uniform	Non-Uniform	Uniform	Non-Uniform	Uniform	Non-Uniform
$\eta$	0.65	0.65	0.65	0.65	0.65	0.65
R	$3\times {10}^6$	$3\times {10}^6$	$3\times {10}^6$	$3\times {10}^6$	$3\times {10}^6$	$3\times {10}^6$
$\tau$	2000	2000	2000	2000	2000	2000
$P_r$	0.3	0.3	1	1	5	5
$N_{\rho }$	3	3	3	3	3	3
n	2	2	2	2	2	2
$P_m$	1	1	4	4	10	10
${\overline{E}}_{tot}$	11,785.6	18,627.4	846.488	1006.96	29.1356	20.9883
${\overline{E}}_p$	23.7465	64.0577	0.97983	4.29216	0.0215626	0.0870345
${\overline{E}}_t$	5930.15	11,540.3	182.149	570.75	4.7595	1.32042
${\check{E}}_p$	1946.24	1746.21	271.194	175.511	11.0902	6.81947
${\check{E}}_t$	1574.74	2261.39	207.107	218.885	9.35735	11.5854
${\overline{M}}_{tot}$	1189.03	0.0662693	464.484	14.92	21.4386	0.638192
${\overline{M}}_p$	6.26582	0.00139628	4.48705	0.117454	0.463097	0.0718596
${\overline{M}}_t$	13.3781	0.0012872	6.38913	0.129885	0.43025	0.0194026
${\check{M}}_p$	256.771	0.0164442	111.957	4.44414	5.34042	0.181853
${\check{M}}_t$	290.496	0.0153525	105.607	3.82462	4.41429	0.0936899
$Rm$	153.5291	193.0150	164.5828	44.8767	76.3350	64.7890
$Ro$	0.1535	0.1930	0.0411	0.0448	0.0076	0.0064
$Lo$	0.0487	0.0003640	0.0304	0.0054	0.0065	0.001129

, we list the six main dynamo solutions on which we focus our attention for most of this article. The table also provides the values of some of their most most important global average characteristics, including energy density components that characterize the various components of the flow. The mean and fluctuating toroidal and poloidal components of the total kinetic energy ($E_{\mathrm{tot}}$) are defined as

$$\begin{array}{cc}\hfill {\overline{E}}_p=\langle {\left(\nabla \times (\nabla \overline{v}\times \mathit{r})\right)}^2/\left(2\overline{\rho }\right)\rangle ,& {\overline{E}}_t=\langle {\left(\nabla r\overline{w}\times \mathit{r}\right)}^2/\left(2\overline{\rho }\right)\rangle ,\end{array}$$

(8a)

$$\begin{array}{cc}\hfill {\check{E}}_p=\langle {\left(\nabla \times (\nabla \check{v}\times \mathit{r})\right)}^2/\left(2\overline{\rho }\right)\rangle ,& {\check{E}}_t=\langle {\left(\nabla r\check{w}\times \mathit{r}\right)}^2/\left(2\overline{\rho }\right)\rangle ,\end{array}$$

(8b)

where angular brackets denote averages over the spherical volume of the shell, bars denote axisymmetric parts and check marks denote non-axisymmetric parts of a scalar field. The total magnetic energy ($M_{\mathrm{tot}}$) can be split in a similar way, with components defined as in Equations (8a) and (8b) but with h and g replacing v and w, respectively, and without the ${\overline{\rho }}^{-1}$ factor within the angular brackets. The total energies are, of course, the sum of all components.

Figure 1. Radial profiles for different model parameters. (a) Non-uniform viscosity and entropy diffusivity vary relative to the density. (b) Non-dimensional parameters $\mathrm{R}$, $\mathrm{Pr}$ and $\tau$ vary when uniform and non-uniform viscosity and entropy diffusivity are considered. For display purposes, all parameters were normalized with respect to their values at the centre of the spherical shell. This allows for the visualization of how the parameters vary with respect to the uniform profile.

Figure 1. Radial profiles for different model parameters. (a) Non-uniform viscosity and entropy diffusivity vary relative to the density. (b) Non-dimensional parameters $\mathrm{R}$, $\mathrm{Pr}$ and $\tau$ vary when uniform and non-uniform viscosity and entropy diffusivity are considered. For display purposes, all parameters were normalized with respect to their values at the centre of the spherical shell. This allows for the visualization of how the parameters vary with respect to the uniform profile.

The parameter values reported in Table 1, along with the values for the solar mass ($M_{\odot }=1.98\times {10}^{30}\mathrm{kg}$), the Carrington sidereal rotation (${\mathrm{\Omega }}_0={\mathrm{\Omega }}_{\odot }=2.87\times {10}^{-6}\mathrm{rad}{\mathrm{s}}^{-1}$), the outer solar radius ($r_o=0.95r_{\odot }=6.59\times {10}^7\mathrm{m}$), the gravitational constant ($G=6.67\times {10}^{-11}{\mathrm{m}}^3{\mathrm{kg}}^{-1}{\mathrm{s}}^{-2}$), and a solar central density estimate (${\overline{\rho }}_c=1.622\times {10}^5{\mathrm{kg}/\mathrm{m}}^3$), can be used with polytropic profile (1) and definition (4) to estimate transport coefficients such as the viscosity, entropy and magnetic diffusivities in dimensional form, as well as derived thermodynamic quantities, such as polytropic constants and luminosities for each of the simulations. Unsurprisingly, comparison with known solar plasma values reveals differences, which are justified in the preceding discussion. The values we use must be interpreted as effective “eddy” coefficients that—at least in part—account for unresolved turbulent subgrid scales.

3.2. Direct Comparison between Simulated Differential Rotation and Observations

To assess the effects of non-uniform viscosity and entropy diffusivity profiles on the overall shape and structure of the differential rotation, it is desirable to compare the solutions of the model to the observed true profile of the solar differential rotation. The latter is well measured [8], as previously mentioned, and the angular velocity ($\mathrm{\Omega }$) is shown in Figure 2a courtesy of Howe [9]. The strongest rotation is concentrated near the equator, with contours of constant rotation closely following lines of oriented 25${}^{\circ }$ to the axis of rotation. It is these feature that give rise to what we refer to as “conical” differential rotation (as opposed to “cylindrical” rotation parallel to the axis of rotation). To aid in our numerical implementation, we chose to work with a simpler but, for our purposes, sufficiently accurate analytical approximation constructed by Kosovichev [37]. This is illustrated in Figure 2b, where we omitted the quiescent radiative interior, the thin solar near-surface shear layer and the tachocline layer, as these are not captured by our mathematical model. Since Equations (2a)–(2d) of our model are formulated with respect to a rotating coordinate system, we subtracted the rigid frame rotation of the Sun (${\mathrm{\Omega }}_{\odot }=870\pi$ nHz). Figure 2c shows the corresponding zonal (linear) velocity ($u_{\phi }=rsin\theta (\mathrm{\Omega }-{\mathrm{\Omega }}_{\odot })$, where the distance from the centre (r) is measured in the dimensional units of our model (shell thickness)).

Figure 3 shows a direct comparison of the azimuthally averaged zonal velocity (${\overline{u}}_{\phi }$) of our model solutions with the latter profile (Figure 2c). Our primary interest here is the shape of the differential rotation profile, so we rescaled the reference profile to obtain the same maximum at the equator as the simulation results. We scaled with respect to the equatorial maximum, since our model is primarily concerned with the profile in the tangent cylinder.

For the smallest value of the Prandtl number ($\mathrm{Pr}=0.3$) and with uniform profiles of $\overline{\nu }$ and $\overline{\kappa }$, Figure 3(Aa) shows that the differential rotation is geostrophic, with contour lines of ${\overline{u}}_{\phi }$ parallel to the rotation axis and a prominent prograde jet filling the region outside of the tangent cylinder (the coaxial cylinder with frame rotation axis and tangential to the inner core of the shell at the equatorial plane), while within the tangent cylinder, the zonal velocity vanishes. This pattern compares poorly with the solar profile shown in Figure 3(Ab), which is also reflected in Figure 3(Ac), where errors appear across all latitudes. The inclusion of radially non-uniform profiles of $\overline{\nu }$ and $\overline{\kappa }$ seems to have a negligible effect on the differential rotation, which retains the structure just described (Figure 3B). Somewhat surprisingly, little further change to the pattern is seen at the moderate value of the Prandtl number ($\mathrm{Pr}=1$) in both uniform and non-uniform profiles of $\overline{\nu }$ and $\overline{\kappa }$ (Figure 3C,D). At a large value of the Prandtl number ($\mathrm{Pr}=5$), the contours of zonal velocity start to deviate from a cylindrical shape (Figure 3(Ea)). The effect is strongly amplified in the case of non-uniform viscosity and entropy diffusivity (Figure 3(Fa)) and assumes a shape very similar to that of the observed solar rotation profile shown in Figure 3(Fb), resulting in vanishingly small discrepancies at mid latitudes (Figure 3(Fc)). It is notable that in the latter case, differential rotation develops within the tangent cylinder at high latitudes and in the polar regions. While the structure of the polar differential rotation is more complicated than the observed rotation, it is notable that it features a relatively large retrograde polar jet, as does the true solar zonal flow. We note that comparison in the polar regions is subject to greater uncertainty. Firstly, in simulations, numerical error increases when the distance from the axis of rotation tends toward zero, which occurs in the polar region. Secondly, the polar regions of the Sun are not in a direct line of sight, so observational measurements in these regions are less accurate.

In summary, the simulation that produces the most solar-like differential rotation is that described in column F of Table 1. Details of this simulation are displayed in Figure 3F, Figure 4F and Figure 5F and discussed further below.

4. Discussion including Further Results

4.1. Structure of the Flow and Thermal Wind Balance

In order to understand the increasing agreement between observed and measured zonal flows in the presence of non-uniform viscosity and entropy diffusivity at larger values of the Prandtl number, in Figure 4, we visualize the rest of the fluid flow components for the same simulations discussed in relation to Figure 3. Outside the tangent cylinder, convection occurs in the form of rolls aligned with the axis of rotation. These exhibit properties of thermal Rossby waves in that they drift in the prograde azimuthal direction. Since the variation in the direction of the axis of rotation is minimized by this form of convection, it is best visualized by the streamlines in the equatorial plane, as shown in the third and sixth columns (c) of Figure 4. These equatorial cross sections indicate the spiralling shape of the convection columns, which becomes less pronounced with the increase in the Prandtl number ($\mathrm{Pr}$). The convection columns extend to fill the entire tangent cylinder, as visible in the plots of the radial velocity projected on the spherical surface at mid-shell (second and fifth columns (b) of Figure 4). The most striking effect of the presence of non-uniform profiles of $\overline{\nu }$ and $\overline{\kappa }$ in comparison with the corresponding uniform-profile simulations is the development of dominant convection flow in the polar regions inside the tangent cylinder. As gravity and rotation vectors are nearly parallel in the polar regions, convection resembles that realized in a horizontal layer heated from below and rotated about a vertical axis but now modulated by interactions with the convection outside of the tangent cylinder. As noted in relation to Figure 3, differential rotation in the polar regions is retrograde, therefore tending to reduce the rotational constraint and forming a feedback loop, supporting the polar flow. At the largest value of the Prandtl number ($\mathrm{Pr}=5$), the polar flows become particularly strong and organized into elongated rolls in the shape of bicycle wheel spokes emanating from the axis of rotation and gently spiraling in the retrograde direction towards the periphery opposite to the drift of the equatorial structures.

The presence of relatively vigorous and well-organized polar convection suggests that entropy transport is enhanced in the polar regions compared to that in the equatorial regions. To verify this, in Figure 5, we plot the gradient of entropy with respect to latitude for all cases discussed in Figure 3 and Figure 4.

There is a clear distinction between the uniform cases (Figure 5A,C,E) and the non-uniform cases (Figure 5B,D,F), with a break in positive ${\langle S\rangle }_{\phi ,t}$ from mid to high/low latitudes. This result is significant because deviations of differential rotations form geostrophy (cylindrical profiles) are known to also be enhanced by the presence of non-vanishing latitudinal entropy gradients [1,30,38]. Indeed, when pressure gradients and Coriolis and buoyancy forces are in a dominant balance, as occured in our simulations and may indeed occur in the bulk of the solar convection zone, it can be shown that the zonal component of the curl of the momentum Equation (2b) reduces to the so-called thermal wind balance

$$\begin{array}{c}\widehat{k}·\nabla {\langle {\mathbf{u}}_{\phi }\rangle }_t\propto \frac{\partial {\langle S\rangle }_{\phi ,t}}{\partial \theta },\end{array}$$

(9)

which is a generalisation of the Taylor–Proudman theorem for rotating fluids in the presence of buoyancy. This relation shows that in or close to an adiabatic state, if $\partial {\langle S\rangle }_{\phi ,t}/\partial \theta \approx 0$, then the rotation profile must be close to cylindrical and when, on the other hand, significant latitudinal entropy gradients are present, as in the cases of non-uniform viscosity and diffusivity profiles shown in Figure 5, where non-cylindrical differential rotation is promoted.

To summarise, the radially non-uniform viscosity and entropy diffusivity profiles allow enhanced convection to develop at the poles, resulting in a non-vanishing entropy gradient with increased latitude. In turn, this helps to produce more solar-like differential rotation according to Equation (9).

4.2. Secondary Considerations

In this preliminary study, we did not fully explore the behaviour model solutions in the accessible parameter space, as we wished to focus our attention on the main new extensions of our model, i.e., the introduction of non-uniform $\overline{\nu }$ and $\overline{\kappa }$ profiles. We are, however, aware of other effects and processes that may alter the results described above, at least quantitatively. We briefly touch upon two of them now.

Solar–antisolar transition. Firstly, in Figure 3, Figure 4 and Figure 5 we restricted the discussion to a few cases with an identical but fixed value/profile of the Rayleigh number. The Rayleigh number measures the contribution of buoyancy to the overall balance of forces in the system; even in the presence of uniform viscosity and entropy diffusivity, it has a profound effect on convection. This dependence has been studied extensively using linear instability analysis, as well as finite-amplitude simulations under various conditions and approximations. Some references to our own works are [31,39,40], and recent review papers include [2]. Here, we only remark that with the increase in driving due to buoyancy, an abrupt transition occurs from prograde solar-like differential rotation to retrograde antisolar differential rotation in the equatorial region. Figure 6 illustrates the phenomenon for a set of cases with increasing Rayleigh number. This transition has been known since the early works reported in [41] and has recently attracted attention from a number of authors in the context of solar and stellar convection [17,33]. Since buoyancy is expected to predominate in the force balance of solar convection, the results reported above for a moderate value of the Rayleigh number may now appear transient. However, the critical value of the Rayleigh number ($\mathrm{R}$) for transition from solar to antisolar rotation is a function of the remaining parameter values. For instance, while we were able to observe a transition at $\mathrm{Pr}=0.3$, as shown in Figure 6, no such transition was observed at $\mathrm{Pr}=1$ and $\mathrm{Pr}=5$. The transition depends, for instance, on the value of the background density stratification parameter ($N_{\rho }$), which, in our simulations, was significantly far removed from estimates for the solar convection zone. Hence, we hypothesize that in the regimes of prograde equatorial rotation, the effects reported in our work still hold true. Another observation to note in Figure 6 is that the strong buoyancy driving convection towards turbulent states is not sufficient to produce conical differential rotation on its own; in fact, before and after, the transition rotation remains largely geostrophic in the illustrated uniform profile cases.

Effects of self-sustained magnetic fields. So far, we have omitted from the discussion the effects of self-sustained magnetic fields on the convective flow and the associated differential rotation. Figure 7 and Figure 8 show a comparison of a non-magnetic convection simulation with a dynamo solution obtained with identical parameter values and with a magnetic Prandtl number of $\mathrm{Pm}=10$. The value of the magnetic Prandtl number was selected so as to ensure that a non-decaying magnetic field was obtained. Figure 7, in particular, shows a direct comparison of the kinetic energy components of the two cases. The main effect of the self-generated and self-sustained magnetic field in the dynamo simulation is a reduction of the kinetic energy component associated with the zonal flow—by roughly 25% in this case. This is due to the fact that even in this significantly diffusive case, the poloidal components of the magnetic field behave, to some extent, as frozen in the fluid and thus act to oppose and slow down zonal jets. Zonal flow, by its nature, does not convect heat between the spherical boundaries but distorts and suppresses convective motions. Thus, because the magnetic field now acts to reduce zonal flow, the mean and fluctuating poloidal and toroidal kinetic energies are somewhat larger but not to a significant degree. Figure 7 further demonstrates the chaotic temporal behaviour of the solutions. Figure 8 shows a side-by-side comparison of the spatial structure of the non-magnetic and dynamo solutions in this case. Apart from a small difference in differential rotation isocontours, which are slightly more “pinched” in the non-magnetic case, the rest of the flow components appear nearly indistinguishable, including dominant azimuthal, radial and latitudinal wave numbers and typical length scales. Therefore, we must conclude that the self-sustained magnetic field does not substantially alter the convective and zonal flows. Additionally, Figure 9 illustrates the morphology of the generated magnetic field. The dynamo exhibits a dominant bipolar symmetry in the polar regions, with large patches of magnetic field of the opposite polarity situated at the north and south “caps” of the spherical shell, which remain largely unchanged over time (snapshots in time are not shown but available and analysed). These polar magnetic fluxes are also evident in the azimuthally averaged toroidal and poloidal field lines plotted in the meridional plane. In the vicinity of the equatorial plane up to mid latitudes, the magnetic field has a more complex quadrupolar symmetry with magnetic field patches of the same polarity across the equatorial plane. These quadrupolar patches exhibit a dominant azimuthal wave number ($m=3$) and drift in longitude and oscillate in latitude. Although this is, therefore, by no means a close model of the solar dynamo, we cannot help but point out promising similarities with the latter, for example, the overall bipolar structure, rudiments of cyclicity and active longitudes, not unlike those observed on the Sun [42]. The oscillation period is directly related to the amplitude of the differential rotation (see [43]).

5. Conclusions

In this work, we have presented an initial investigation into the effects of non-uniform viscosity and entropy diffusivity on differential rotation generated by convection in rotating spherical shells. Our motivation was to address the so-called convective conundrum, in which many models of the solar dynamo do not accurately reproduce the profile of solar differential rotation as measured using helioseismology. Our approach was to focus on a small but representative region of the parameter space in order to gain a foothold, establishing a relationship between non-uniform viscosity and entropy diffusivity and differential rotation. Upon keeping all parameter values constant, with the exception of the Prandtl number, our main result is that values of the Prandtl number somewhat larger than unity ($\mathrm{Pr}=5$ in our case) result in a differential rotation profile whose agreement with the observed profile is much improved compared to cases with lower values of the Prandtl number and cases using only uniform profiles of the viscosity and entropy diffusivity. The results are in agreement in the equatorial region up to mid latitudes, but discrepancies remain towards the poles, where both simulations and observations are less accurate.

For the non-uniform viscosity and entropy diffusivity profiles, the differential rotation profile is closely connected to the presence of stronger convection at the poles for higher values of $\mathrm{Pr}$. The increased entropy at the poles leads to a non-trivial entropy gradient in $\theta$, which, in turn, plays a key role in the thermal wind balance, which prevents a Taylor–Proudman state. The result is that conical rather than cylindrical profiles of differential rotation are preferred in the presence of sufficiently large radially non-uniform profiles of viscosity and entropy diffusivity. It is remarkable that spherically symmetric radial non-uniformity can lead to latitudinal non-uniformity, which, in contrast to the simulations reported in [13], develops self-consistently.

It is somewhat surprising that a large value of the Prandtl number leads to better agreement with observations. On one hand, the values of $\mathrm{Pr}$ in the Sun are expected to be very small based on molecular estimates [23], but on the other hand, turbulent mixing suggests values of the order of unity. The agreement may, in part, be due to the region of the parameter space under consideration, which was chosen so as to capture the presumed force balance in the solar convection zone, as well as for computational feasibility. Indeed, capturing appropriate force balances may be more meaningful than trying to adhere strictly to poorly estimated and difficult-to-achieve parameter values, as recognized in similar studies of the geodynamo [44]. Our results show that solar-like differential rotation is possible in simulations, providing a good starting point for more detailed parameter sweeps—both in terms of the fundamental non-dimensional parameters and the profiles of the non-uniform viscosity and entropy diffusivity. This will allow for a better understanding of how differential rotation affects the dynamo mechanism, which is a potential direction for further research.