1 Introduction↩︎

Convection in stellar and giant planet interiors is extremely turbulent, which poses a practical barrier to numerical simulation. To retain computational feasibility, simulations of this process are typically conducted in a parameter regime far from reality [1], [2], with the results subsequently extrapolated to real objects through scaling analysis and by appealing to theory [3]. In order to achieve this, it is desirable that the simulations exhibit the same scaling behaviour as would a real planet or star. It is typically assumed that the bulk properties of the convection should not depend on the molecular diffusivities, which are small in real astrophysical environments [4], [5]. This assumption underpins the derivation of mixing length theory (MLT; [6]), which has been applied with noteworthy success as a parametrisation of convective heat transport in stellar evolution models (e.g., [7]).

However, diffusivity independent (i.e., ‘diffusion-free’) dynamics are rarely obtained in numerical simulations. For example, the vast majority of shell convection simulations yield scaling behaviour for the radial heat transport that depends upon the diffusivity [8]–[13], although a notable exception is presented by [11], where a subset of their simulations do obtain a diffusion-free scaling for cases with very rapid rotation (paired with moderate supercriticality). The influence of diffusivity in simulations of shell convection also affects their bulk kinematic properties. For example, simulations of stellar convection with high diffusivity have been shown to generate prograde, ‘Solar-like’ equatorial zonal flows (or ‘differential rotation’), but they have a tendency to switch to a regime of retrograde, ‘anti-Solar’ flow as the diffusivities are reduced [14]–[18]. Moreover, [11] obtain a diffusivity dependent scaling for the Reynolds number even for simulations that are highly supercritical.

Driving convection using a conductive heat flux through the inner and outer radial boundaries is one cause of diffusivity dependence in numerical simulations. When convection is boundary-driven, thin thermal boundary layers form, and the heat transport through the domain is throttled by the conductive flux through the boundary layer ([11]; see also the review by [19]). Regarding the differential rotation, the transition from prograde to retrograde motion is understood to occur when the convective Rossby number, which for a fixed rotation rate generally increases as the diffusivities are reduced (in simulations with diffusion-dependent dynamics), exceeds unity [15], [20]. [21] provide a geometric interpretation for this transition, while other authors (e.g., [18]) argue it is directly related to the thinning of the thermal boundary layer. It is interesting to note that simulations of deep stellar convection driven entirely by internal heating and cooling also exhibit diffusivity dependent differential rotation [22]–[24]; however, in these models conductive heat transport is likely still important near the upper boundary, owing to the thinness of the cooling layer.

It can be argued ‘boundary-driven’ simulations are not particularly relevant to real astrophysical environments, as the boundary layers that dominate the dynamics over large swathes of parameter space likely bear little resemblance to those of the real objects [5], [18]. This has led some authors to pursue an alternative approach, where convection is driven by radiative heating and cooling (with a specified spatial structure) applied directly to the fluid. In the context of rapidly rotating \(f\)-plane convection in a Cartesian geometry, [5] demonstrate that this approach yields diffusion-free bulk scaling behaviour for both thermal (e.g., the convective heat transport) and kinematic (e.g., the fluctuating kinetic energy) properties. [25] generalise this result to the case where the directions of the acceleration due to gravity and the rotational axis are misaligned (as is the case for the non-polar latitudes of a sphere). However, [26] illustrate that it is possible to simultaneously obtain diffusion-free heat transport, but diffusion-dependent convective velocities (using simulations of internally heated and cooled, rotating convection in an idealised Cartesian geometry with \(2.5\) dimensions). Additionally, they show that the realisation of diffusion-free dynamics in their model is sensitive to the kinematic boundary condition (free-slip configurations more readily yield diffusion-free behaviour). [27] show that diffusion-free heat transport can be obtained in Cartesian simulations of non-rotating convection driven by internal heating and cooling (i.e., a realisation of the ‘ultimate regime’). Finally, laboratory experiments driven by internal heating and cooling have obtained diffusion-free heat transport for both the rotating case and the non-rotating case [28]–[30].

Building on these studies, our objective is to investigate the dynamics and scaling behaviour of convection driven by internal heating and cooling in a spherical geometry. We aim to determine whether internal heating and cooling gives rise to ‘diffusion-free’ dynamics, focusing on the convective heat transport, convective velocities, and the zonal flows (differential rotation). We also investigate the extent to which our simulation results can be described by standard scaling theories (whether diffusion-free or otherwise). In this work, we study the dynamics of a Boussinesq, hydrodynamic fluid, deferring inclusion of the effects of compressibility and magnetism to future work.

2 Methods↩︎

We study the dynamics of Boussinesq, rotating spherical shell convection driven by internal heating and cooling. Our simulations are conducted using the pseudo-spectral code Dedalus [31]. All variables are expended in terms of spherical harmonics in the horizontal directions and Chebyshev polynomials in the radial direction using a standard 3/2-dealisaing.

2.1 Model↩︎

The governing equations are \[\begin{align} \frac{\text{D}\boldsymbol{u}}{\text{D}t} + 2\mathbf{\Omega}\times\boldsymbol{u} &= -\mathbf{\nabla}p + \tilde{g}(r)T\mathbf{e}_{r} + \nu\nabla^{2}\boldsymbol{u} \\ \mathbf{\nabla}\cdot\boldsymbol{u}&=0\\ \frac{\text{D}T}{\text{D}t} &= q(r) + \kappa\nabla^{2}T, \end{align}\] where \(\boldsymbol{u}=(u_{r}, u_{\theta}, u_{\phi})\) is the fluid velocity in spherical polar coordinates, \(p\) is a pressure, and \(T\) is a scaled temperature variable, related to the real temperature by \(T = g\alpha \delta T_{\text{real}}\) (where \(g\) is the magnitude of gravitational acceleration and \(\alpha\) is the coefficient of thermal expansion, and \(\delta T_{\text{real}}\) is the Boussinesq temperature perturbation). The symbol \(\nu\) denotes the kinematic viscosity and \(\kappa\) denotes the thermal conductivity. The dimensionless gravity profile is given by \(\tilde{g}(r)=\left(r_{\text{o}}/r\right)^{2}\) where \(r_{\text{o}}\) is the outer radius of the convection zone (see below for definition). The unit vector in the radial direction is given by \(\mathbf{e}_{r}\), and the rotation vector \(\mathbf{\Omega}=\Omega_{0}(\cos\theta, -\sin\theta,\;0)\), where \(\theta\) is co-latitude and \(\Omega_{0}\) is the rotation rate.

Convection is driven by a prescribed heating and cooling function \(q(r)\), configured so that heat is deposited and removed from regions of depth \(\delta\) at the base and top of the simulated domain, respectively. We use the notation \(r_{\text{i}}\) and \(r_{\text{o}}\) to denote the radii of the inner and outer boundaries of the convection zone (CZ), within which there is no imposed heating or cooling. Mathematically, \(q\) is defined (following [5], [25]): \[q(r)=\frac{F}{\delta}\frac{r_{\text{i}}^{2}}{r^{2}}\begin{cases} 1+\cos\left[\frac{2\pi\left(r-r_{\text{i}}+\frac{\delta}{2}\right)}{\delta}\right]\kern-3pt; & r_{\text{i}}-\delta\leq r\leq r_{\text{i}}, \\ 0; & r_{\text{i}} < r < r_{\text{o}}, \\ -1-\cos\left[\frac{2\pi\left(r-r_{\text{o}}-\frac{\delta}{2}\right)}{\delta}\right]\kern-3pt; & r_{\text{o}}\leq r\leq r_{\text{o}}+\delta, \end{cases}\] where \(F\) is the flux injected at the base of the CZ (corresponding to a shell-integrated flux \(\mathscr{F}_{\text{tot}}=4\pi r_{\text{i}}^{2}F\); \(\mathscr{F}_{\text{tot}}\) is independent of \(r\) within the CZ). The geometric prefactor in the definition of \(q\) ensures that the cooling layer removes the same amount of heat as is injected by the heating layer.

At the inner and outer radii of the simulated domain (located at \(r_{\text{i}}-\delta\) and \(r_{\text{o}}+\delta\), respectively) we apply a free-slip, impenetrable, and insulating boundary condition (consistent with \(\iiint_{V} q(r)\,\text{d}V=0\)). Applying the heating and cooling function \(q(r)\) in combination with these boundary conditions has the effect of replacing the conductive boundary layers that exist in boundary driven simulations with heating and cooling regions of fixed size.

2.2 Non-dimensionalisation↩︎

To non-dimensionalise the governing equations, we use the depth of the convection zone \(d\equiv r_{\text{o}}-r_{\text{i}}\) as a reference length scale, a flux-based temperature scale, given by \(\Lambda\equiv Fd/\kappa\), and a timescale given by the free-fall time \(\tau=\sqrt{d/\Lambda}\).

Under this non-dimensionalisation, the dynamics of the system are determined by: a flux-based Rayleigh number, \(Ra_{\text{F}}\); the Taylor number, \(Ta\); the Prandtl number, \(Pr\); the CZ radius ratio, \(\eta\); and the non-dimensional heating/cooling depth \(\tilde{\delta}\). These numbers are defined by: \[\begin{align} &Ra_{\text{F}}\equiv\frac{d^{4}F}{\nu\kappa^{2}},\quad Ta\equiv\frac{4\Omega_{0}^{2}d^{4}}{\nu^{2}},\quad Pr\equiv\frac{\nu}{\kappa},\\ &\eta\equiv\frac{r_{\text{i}}}{r_{\text{o}}},\;\text{and}\quad \tilde{\delta}\equiv\frac{\delta}{d}. \end{align}\] In all simulations, we set \(Pr=1\), \(\eta=0.8\), and \(\tilde{\delta}=0.2\). From the definitions of \(Ra_{\text{F}}\), \(Ta\), and \(Pr\), it is possible to define a flux-based convective Rossby number (e.g., [8], [18], [32], [33]), given by: \[Ro_{\text{cv,F}}\equiv\left(\frac{Ra_{\text{F}}}{Ta^{\frac{3}{2}}Pr^{2}}\right)^{\frac{1}{3}}=\left(d F\right)^{\frac{1}{3}}\frac{1}{2\Omega_{0} d}. \label{eq:Rocv}\tag{1}\] which, notably, is independent of the diffusivities \(\nu\) and \(\kappa\). Therefore, for a fixed flux, variation of \(Ro_{\text{cv,F}}\) corresponds to varying \(\Omega_{0}\), while varying \(Ra_\text{F}\) and \(Ta\) so as to keep \(Ro_{\text{cv,F}}\) fixed is equivalent to varying \(\nu\) and \(\kappa\) (recall \(\kappa=\nu\) when \(Pr=1\)). We note that this definition of \(Ro_{\text{cv,F}}\) is closely related to estimates of the Rossby number commonly employed in stellar astronomy (e.g., [34]), \(Ro\,{\sim}\,1/(\Omega_{0}\tau_{\text{c}})\), assuming that the convective overturning time \(\tau_{\text{c}}\) is given by MLT.

2.3 Summary of experiments and numerical details↩︎

Results from 38 simulations are presented. The parameter survey has been constructed to produce sets of simulations with constant \(Ta\) (within which \(Ra_{\text{F}}\), and thus \(Ro_{\text{cv,F}}\), varies), as well as some sets with constant \(Ro_{\text{cv,F}}\) (within which \(Ra_{\text{F}}\) and \(Ta\) vary). The parameter survey spans a range of \(Ro_{\text{cv,F}}\) between \(0.01\) and \(3.1\), \(Ta\) between \(10^{5}\) and \(10^{9}\), and \(Ra_\text{F}\) between \(3\times10^{5}\) and \(9.5\times10^{9}\). A full list of our experiments is given in Appendix 5 (along with details of the resolution used), and a visual summary is shown in Figure 1.

Figure 1a shows the input parameters \(Ra_{\text{F}}\) and \(Ta\), Figure 1b shows \(Ra Ta^{-\frac{2}{3}}\) as a proxy for the supercriticality \(Ra/Ra_{\text{c}}\) (\(Ra_{\text{c}}\,{\sim}\,Ta^{\frac{2}{3}}\) for large \(Ta\); [35]), and Figure 1c shows the output Reynolds number to illustrate the degree of turbulence in each experiment. Above, \(Ra\) is the usual Rayleigh number, which is an output parameter for a fixed-flux set-up. It is defined as \[Ra \equiv \frac{d^{3}\Delta T}{\nu\kappa}, \label{eq:Ra}\tag{2}\] where \(\Delta T=T(r\,{=}\,r_{\text{o}}) - T(r\,{=}\,r_{\text{i}})\) is the output radial temperature contrast, evaluated between the boundaries of the CZ. The Reynolds number is defined \[Re \equiv \frac{UL}{\nu} \label{eq:Re}\tag{3}\] where we take \(U\) to be the root-mean-square (r.m.s.) of the output fluctuating (azimuthal average removed) velocity vector within the CZ, and \(L=d\). Also shown in Figure 1d are two different output Rossby numbers, defined \[Ro_{\text{bulk}} = \frac{U}{\Omega_{0} d}\quad\text{and}\quad Ro_\omega = \frac{\omega}{\Omega_{0}}, \label{eq:Ro95alt}\tag{4}\] where for \(Ro_{\omega}\), \(\omega\) is the r.m.s fluctuating vorticity within the CZ. These two output \(Ro\) are intended to indicate the degree of rotational constraint at the dominant scale of convection, and at small scales, respectively. Scaling results obtained from our experiments will show that it is useful to categorise them in terms of three regimes: (1) a ‘rotationally-constrained’ regime defined by \(Ro_{\text{cv,F}}<0.1\), for which both large and small scales are rotationally constrained (i.e., \(Ro_{\omega}\lesssim1\) and \(Ro_{\text{bulk}}\ll1\)); (2) a ‘rotationally-influenced’ regime defined by \(0.1\leq Ro_{\text{cv,F}}<1\), where only the large-scales are rotationally constrained (\(Ro_{\omega}\gtrsim1\) and \(Ro_{\text{bulk}}<1\)); and (3) a ‘rotationally-uninfluenced’ regime defined by \(Ro_{\text{cv,F}}\geq1\), for which the rotational constraint is lost entirely (\(Ro_{\omega}\gg1\) and \(Ro_{\text{bulk}}>1\)).

All simulations were integrated until the volume averaged kinetic energy and convective heat transport reached a steady state. In practise, this means that each simulation was run for multiple thermal diffusion times (corresponding to \(10^{3}\) – \(10^{4}\) free-fall times). To ensure numerical convergence, the simulations have been inspected to verify closure of the horizontally-averaged heat budget. The kinetic energy spectra have also been inspected to confirm there is no significant accumulation of energy at the smallest scales so that they are spatially resolved.

3 Results↩︎

3.1 Simulation results↩︎

Summary statistics for our experiments are presented in Figure 2. All quantities shown in this figure are pseudodimensional (specifically, they have been re-dimensionalised using \(F=1\) and \(d=1\) to set the units¹), and are plotted against \(Ro_{\text{cv,F}}\) (equivalent to \(1/(2\Omega_{0})\) for our choice of re-dimensionalisation; cf. Equation 1 ). In these panels, diffusivity-dependence manifests as scatter on the \(y\)-axis for a given value \(Ro_{\text{cv,F}}\) (note: not all \(Ro_{\text{cv,F}}\) have multiple simulations; those with only one entry are marked with a \(\dagger\) symbol). Conversely, if the dynamics are diffusion-free, then all the data for a given \(Ro_{\text{cv,F}}\) will collapse onto one point.

Figure 2a shows the shell-averaged temperature difference \(\Delta T\) between the inner and outer radii of the CZ. For \(Ro_{\text{cv,F}}<1\), \(\Delta T\) displays a strong dependence on \(Ro_{\text{cv,F}}\), increasing with \(\Omega_{0}\) as the rotational constraint acts to throttle convective heat transport. For each \(Ro_{\text{cv,F}}\), the scatter between experiments with different \(Ra_{\text{F}}\) is very small. In the rotationally-influenced regime, the dependence of \(\Delta T\) on \(Ro_{\text{cv,F}}\) (\(\Delta T\propto\Omega_{0}^{\frac{4}{5}}\,{\sim}\,Ro_{\text{cv,F}}^{-\frac{4}{5}}\)) is consistent with rotating mixing length theory (RMLT; see Section 3.2; [5], [36]), which is a diffusion-free theory. For the cases categorised as rotationally-constrained with \(Ro_{\text{cv,F}}<0.1\), \(\Delta T\) depends more strongly on \(Ro_{\text{cv,F}}\) than predicted by RMLT, while remaining independent of the diffusivity (for the values of \(Ro_{\text{cv,F}}\) for which there are multiple experiments). In the rotationally-uninfluenced regime where \(Ro_{\text{cv,F}}>1\), the rotational constraint is lost, causing \(\Delta T\) to become independent of \(Ro_{\text{cv,F}}\). Our simulations yield results that are consistent with this expectation, and thus by extension are consistent with diffusion-free behaviour.

The azimuthal velocity averaged within \(\pm5^{\circ}\) of the equator and throughout the depth of the CZ, \(U_{\text{eq}}\), is shown in Figure 2b. The dependence of \(U_{\text{eq}}\) on \(Ro_{\text{cv,F}}\) we obtain is consistent with previous work [15]. In the absence of a rotational constraint (\(Ro_{\text{cv,F}}>1\)), differential rotation does not develop (\(U_{\text{eq}}\approx0\)). In the rotationally-influenced regime (\(Ro_{\text{cv,F}}<1\); \(Ro_{\omega}\gtrsim1\)), anti-Solar (retrograde) differential rotation develops, with an amplitude that increases as the rotation rate is increased. In this regime, the amplitude of the equatorial zonal flow is independent of the diffusivity. When rotation becomes dominant (\(Ro_{\text{cv,F}}<0.1\); \(Ro_{\omega}\lesssim1\)), the differential rotation switches direction from retrograde to prograde (Solar-like). In this regime, the data in general do not collapse onto a single curve for \(U_{\text{eq}}\) vs. \(Ro_{\text{cv,F}}\) indicating that, at least for the supercriticalities considered in the present work, the strength of \(U_{\text{eq}}\) depends on the diffusivities. We note that there are only two values of \(Ro_{\text{cv,F}}<0.1\) for which we ran a significant number of experiments with varying \(Ra_{\text{F}}\) (5 experiments each for \(Ro_{\text{cv,F}}=3.11\times10^{-2}\) and \(6.70\times10^{-2}\)). For the series with \(Ro_{\text{cv,F}}=6.70\times10^{-2}\), \(U_{\text{eq}}\) does not vary significantly as the diffusivity is varied, but for the series with \(Ro_{\text{cv,F}}=3.11\times10^{-2}\) there is a diffusivity dependence.

Figures 2c and d show the kinetic energy (KE) averaged over the volume of the CZ. Figure 2c shows \(\text{KE}_{\text{fluc}}\), computed using the r.m.s. fluctuating velocity, and Figure 2d shows the ratio DRKE/KE, where KE is the total kinetic energy, and DRKE is the ‘differential rotation kinetic energy’, computed using the axially-symmetric component of the azimuthal velocity. The behaviour of \(KE_{\text{fluc}}\) in Figure 2c tells a similar story to \(U_{\text{eq}}\), namely that the kinematic properties of the flow display a weak dependence on the diffusivity for \(Ro_{\text{cv,F}}>0.1\) (rotationally-uninfluenced or rotationally-influenced), but a stronger dependence once the influence of rotation becomes dominant. There are two additional features of Figure 2c,d that are worth emphasising. First, we note that in the rotationally-influenced regime, the dependence of \(KE_{\text{fluc}}\) on \(Ro_{\text{cv,F}}\) is consistent with the diffusion-free ‘Coriolis-inertial-Archimedean’ (CIA) scaling (e.g., [11]; see Section 3.2). Second, we highlight that while \(\text{KE}_{\text{fluc}}\) (Figure 2c) and the DRKE depend on the diffusivity in the rotationally-constrained regime, they do so in such a way that the ratio DRKE/KE is roughly diffusion-free (Figure 2d). In this regime, the DRKE dominates the total KE so that their ratio saturates at \(\text{DRKE}/\text{KE}\approx1\).

The dependence of the summary statistics shown in Figure 2 on \(Ro_{\text{cv,F}}\) (and the diffusivities, where applicable) arise from the influence of these parameters on the convective flow patterns. To illustrate this, Figure 3 shows the perturbation temperature \(T^{\prime}\) (horizontal average removed) for six selected experiments. The upper row shows two experiments with \(Ro_{\text{cv,F}}=1.44\) that occupy the rotationally-uninfluenced regime, the middle row shows two experiments with \(Ro_{\text{cv,F}}=3.11\times10^{-1}\) that fall within the rotationally-influenced regime, and finally the bottom row shows two experiments with \(Ro_{\text{cv,F}}=6.70\times10^{-2}\) that occupy the rotationally-constrained regime.

In the rotationally-uninfluenced regime, convective plumes extend outwards in the radial direction at all latitudes and the effect of rotation is undetectable. In this regime, the temperature field is organised into large ‘warm’ and ‘cool’ patches that are dominated by spherical harmonic degree \(\ell=1\) or \(2\). The morphology in the rotationally-influenced cases with \(Ro_{\text{cv,F}}=3.11\times10^{-1}\) is similar, although the convection at higher latitudes now exhibits a weak alignment with the rotation axis and some low-latitude features are elongated in the azimuthal direction. However, the impact of rotation is hard to detect, consistent with \(Ro_{\text{bulk}}>1\). Figure 4 shows an additional case in the rotationally-influenced regime with \(Ro_{\text{cv,F}}=1.44\times10^{-1}\) (i.e., a slightly faster rotation rate), which more clearly shows that the convection is aligned cylindrically, and zonally-elongated at lower latitudes. We highlight that the pair of experiments with \(Ro_{\text{cv,F}}=3.11\times10^{-1}\) are morphologically very similar, as are the pair with \(Ro_{\text{cv,F}}=1.44\), which is consistent with the inference from Figure 2 that the diffusivity has a limited influence on the dynamics in the rotationally-uninfluenced and rotationally-influenced regimes.

The effect of rotation in the rotationally-constrained regime is clear. In mid-latitudes and the polar regions, convection is aligned cylindrically (i.e., with the rotation axis). Closer to the equator, large, modulated thermal Rossby wave-like features can be identified, that propagate in the prograde direction (with respect to the mean zonal flow). These waves are present in all experiments with \(Ro_{\text{cv,F}}<0.1\) but in none with \(Ro_{\text{cv,F}}\geq1\). The transition to a regime of equatorial dynamics dominated by wave-like features occurs when the small scales become rotationally constrained (\(Ro_{\omega}<1\)). For the relatively turbulent cases shown in Figure 3 there is a strong polar vortex. This is not present in simulations that are more strongly dissipative (an example is given by the lower panel of Figure 4).

Figure 5 shows cross-sections of the differential rotation, quantified by the azimuthally-averaged angular velocity normalised by the rotation rate \(\Omega_{0}\), \[\frac{\Omega}{\Omega_{0}} = 1 + \frac{\overline{u}_{\phi}}{\Omega_{0} r\sin\theta}, \label{eq:Om}\tag{5}\] for two series of experiments at \(Ro_{\text{cv,F}}=3.11\times10^{-1}\) (rotationally-influenced; shown in the upper row) and \(Ro_{\text{cv,F}}=6.70\times10^{-2}\) (rotationally-constrained; lower row). For a given \(Ro_{\text{cv,F}}\), larger \(Ra_{\text{F}}\) corresponds explicitly to reduced diffusivity (decreasing across panels moving from left to right). The experiments in the \(Ro_{\text{cv,F}}=3.11\times10^{-1}\) set exhibit ‘anti-Solar’ differential rotation with a slow equator and a fast pole (denoted SF). The differential rotation in this set of experiments has little dependence on the diffusivity, with the spatial structure and amplitude of \(\Omega/\Omega_{0}\) remaining essentially unchanged as \(Ra_{\text{F}}\) is increased. In contrast, the differential rotation in the \(Ro_{\text{cv,F}}=6.70\times10^{-2}\) set does depend on the diffusivity. All of these experiments have a fast equator, but the structure of the zonal flow at higher latitudes varies. The most diffusive experiment (left-most panel in Figure 5, bottom row) features a Solar-like fast equator and slow pole (denoted FS), while the remaining four experiments have prograde flow at the equator, flanked by retrograde flow in mid-latitudes, and then finally a polar vortex (cf. Figure 3d) at high-latitudes (i.e., fast–slow–fast; denoted FSF). Once the transition from FS to FSF differential rotation occurs, the \(Ro_{\text{cv,F}}=6.70\times10^{-2}\) experiments appear to converge towards an asymptotic ‘diffusion-free’ structure and amplitude as the diffusivity is reduced further (see also the small spread in \(U_{\text{eq}}\) when \(Ro_{\text{cv,F}}=6.70\times10^{-2}\) in Figure 2).

Each of the scatter charts in this article use the marker shape to illustrate the occurrence of the SF, FSF, and FS configurations of differential rotation in parameter space (circles for SF, pentagons for FSF, and squares for FS). From Figure 1a–c, we identify that solar-like differential rotation (FS) is only obtained for the most diffusive experiments (\(Ra_{\text{F}}\lesssim10^{6}\)), while more turbulent simulations with a strong rotational constraint exhibit FSF differential rotation. Figure 1d shows that the transition from a fast equator (either FS or FSF) to a slow equator occurs at the boundary between the rotationally-constrained and rotationally-influenced regimes, when \(Ro_{\text{cv,F}}=0.1\), which corresponds to \(Ro_{\omega}\approx1\), consistent with previous work [15]. The fact that all of the FS or FSF cases have \(Ro_{\text{cv,F}}<0.1\), and all of the SF cases have \(Ro_{\text{cv,F}}\geq0.1\), implies that the transition from a fast equator to a slow equator is diffusion-free (even if the zonal flow structure in the fast-equator cases themselves displays a dependence on the diffusivity). We believe that this transition is associated with the existence of thermal Rossby wave-like features at low latitudes, which are present in all simulations in the rotationally-constrained regime but absent otherwise. However, a complete analysis of the zonal momentum budget is left for future work.

Finally, we note that the experiments we categorise as rotationally-constrained were run with a lower supercriticality than those categorised as rotationally-influenced or rotationally-uninfluenced (Figure 1b). While we do not observe it here, we would expect the anti-Solar differential rotation in these regimes to display a diffusivity-dependence for sufficiently low supercriticality (i.e., higher diffusivity).

3.2 Comparison with scaling predictions↩︎

3.2.1 Nusselt number scaling↩︎

Radial heat transport is often characterized using the Nusselt Number, defined as the ratio of the shell-integrated total heat flux (convective plus conductive) to the conductive flux. In this work, we evaluate this as an integral over the CZ, which yields: \[Nu\equiv\frac{\mathscr{F}_{\text{tot}}}{\mathscr{F}_{\text{cond}}}=\frac{r_{\text{i}}^{2}}{\int_{r_i}^{r_o}\partial_{r}\langle T\rangle\,r^{2}\text{d}r}. \label{eq:Nu}\tag{6}\] All quantities above are non-dimensional, \(\langle T\rangle\) is the shell-averaged temperature, and the second relation follows from the fact that \(\mathscr{F}_{\text{tot}} = 4\pi r_{\text{i}}^{2}\) within the CZ. Mixing length theory, which essentially equates the kinetic energy of fluid parcels to the buoyancy work done over a characteristic length \(l\) (the mixing length), states that the Nusselt number should scale as: \[Nu\propto Ra^{\frac{1}{2}}, \label{eq:NuMLT}\tag{7}\] for \(Pr=1\). A central assumption in MLT is that the total dimensional heat transport is independent of the diffusivities, and thus the scaling above is diffusion-free. We note that a simple dimensional analysis for a non-rotating system yields the scaling above under this assumption². For a rotating system, the analogous theory to MLT is rotating mixing length theory (RMLT; [36]), which assumes that convection is dominated by modes with a growth rate well described by linear theory and amplitudes determined by a balance between their growth rate and the non-linear cascade rate. RMLT predicts that \(Nu\) scales as \[Nu\propto\frac{Ra^{\frac{3}{2}}}{TaPr^{\frac{1}{2}}}. \label{eq:NuRMLT}\tag{8}\] As with Equation 7 , this is a diffusion-free scaling. It scaling can also be obtained by combining the ‘diffusion-free heat transport assumption’ with a further assumption that the heat transport depends only on the supercriticality of the system \(Ra/Ra_{\text{c}}\) (recall that \(Ra_{\text{c}}\,{\sim}\,Ta^{\frac{3}{2}}\) for large \(Ta\)) [4], [5]. We note that \(Ra\approx Ra_{\text{F}}/Nu\), which yields the scaling \(Nu\propto\left(Ra_{\text{F}}Ta^{-2/3}\right)^{3/5} Pr^{-1/5}\). Using the definition of \(Ro_{\text{cv,F}}\), this can be re-written \(Nu\propto Ra_{\text{F}}^{1/3}Pr^{1/3}Ro_{\text{cv,F}}^{4/5}\). Therefore, with \(Pr\) constant, RMLT implies \[Nu\propto Ra^{\frac{1}{2}}\;\;\text{for constant}\;Ro_{\text{cv,F}} \label{eq:NuRMLTRoCV}\tag{9}\] (i.e., as in MLT). Previous studies of boundary-driven Boussinesq convection in a spherical shell have not yielded MLT-like scaling for the Nusselt number in the regime where rotation is unimportant. Instead, shallower scaling exponents were obtained, due to the throttling of convection by the conductive heat transport in the boundary layers (e.g., [10]). Some studies have obtained RMLT-like scaling in the rapidly rotating regime, but only for extreme parameter values (e.g., [11], wherein experiments with \(Ta\gtrsim10^{12}\) and \(RaTa^{-2/3}\simeq10\) exhibit a Nusselt number scaling in agreement with RMLT).

Figure 6a and 6b show \(Nu\) plotted against \(RaTa^{-2/3}\) in panel a) and \(Ra\) in panel b). Dashed lines corresponding to the MLT and RMLT scaling predictions are included in each panel (note that the MLT scaling in panel a, and the RMLT scaling in panel b, apply to series of experiments with constant \(Ta\)). In panel a), the marker colour is used to show \(Ro_{\text{cv,F}}\).

Inspection of Figure 6a reveals that our experiments fall into three groups distinguished by different \(Nu\) vs. \(RaTa^{-2/3}\) scaling behaviours. When \(Ro_{\text{cv,F}}\geq1\) (indicated by a green marker colour), the Nusselt number follows the MLT scaling prediction (Equation 7 ). These experiments are those that we have categorised as rotationally-uninfluenced. For, \(0.1\leq Ro_{\text{cv,F}}<1\) (indicated by blue marker colours), i.e., experiments in the rotationally-influenced regime, a steeper scaling is obtained that is consistent with RMLT (Equation 8 ). The transition between the MLT and RMLT scalings occurs when \(Ro_{\text{cv,F}}=1\), which roughly corresponds to the point where \(Ro_{\text{bulk}}=1\). The agreement between the scaling obtained by the rotationally-influenced experiments and RMLT is striking, given that the effect of rotation on the convective morphology in this regime can be relatively weak (e.g., Figure 3). These results (as well as those presented in the next subsection) inform our definition of the rotationally-uninfluenced and rotationally-influenced regimes. Finally, in the rotationally-constrained regime for which \(Ro_{\text{cv,F}}<0.1\) (orange marker colours), the Nusselt number scaling departs from the RMLT branch occupied by the rotationally-influenced cases. This departure from RMLT is manifest as an offset from the RMLT branch (reduced \(Nu\) for a given \(RaTa^{-2/3}\)) that increases as \(Ro_{\text{cv,F}}\) is reduced. For constant \(Ro_{\text{cv,F}}\), experiments in this regime roughly follow the \(Nu\propto Ra^{1/2}\) scaling noted above (Equation 9 ; for simulations with sufficient supercriticality, \(RaTa^{-2/3}\gtrsim5\)). This indicates that the convective heat transport is still independent of the diffusivity (consistent with the pseudodimensional \(\Delta T\) shown in Figure 2a), and that the offset from the rotationally-influenced RMLT branch captures a dependence of the convective heat transport on \(Ro_{\text{cv,F}}\) (a diffusion-free parameter) that is not accounted for by RMLT. Hereafter, we refer to this scaling as ‘offset RMLT’. As noted previously, the transition from the rotationally-influenced regime to this new regime occurs when \(Ro_{\omega}\lesssim1\), which indicates that all spatial scales are now affected by rotation, motivating our choice to refer to this regime as rotationally-constrained.

When the same data are plotted against \(Ra\) (Figure 6b), the MLT scaling obtained by the rotationally-uninfluenced regime can be identified for experiments with \(Ta=10^{5}\), which unambiguously follow the \(Ra^{1/2}\) prediction. The \(Ta=10^{6}\) series follows the RMLT scaling for moderate supercriticality, before transitioning towards the MLT scaling at higher supercriticality, and the experiments at \(Ta=10^{7}\) are approaching the RMLT scaling. Finally, the offset RMLT scaling identified for the rotationally-constrained regime is manifest as a steepening of the \(Nu(Ra)\) relation for fixed \(Ta\) (but, again, note that series of experiments with constant \(Ro_{\text{cv,F}}\), joined by black solid lines, follow the RMLT prediction given by Equation 9 ).

A preliminary analysis (not shown) of the latitudinal dependence of the heat budget and Nusselt number (cf. [37]) suggests that the transition to the offset RMLT scaling that characterises the rotationally-constrained regime may be related to the change in the dynamics that from a flow dominated by turbulent convection to a wave-dominated flow (we find wave-like features are always present when \(Ro_{\text{cv,F}}<0.1\), but are absent when \(Ro_{\text{cv,F}}\geq0.1\)). However, the offset scaling could alternatively arise due to interaction between the convection and the zonal flow, although this would have to be a ‘signed’ effect, since strong (anti-Solar) zonal flows are also present in the rotationally-influenced regime for which the standard RMLT scaling is recovered. We will explore this further in future work. We note that the offset scaling persists even for the rotationally-constrained experiments with the highest supercriticality. As a result, we do not believe that it arises because the rotationally-constrained simulations generally have lower supercriticality than those categorised as rotationally-influenced.

3.2.2 Reynolds number scaling↩︎

A common metric of the degree of turbulence in fluid flow is the Reynolds number \(Re\) (given by Equation 3 ), which quantifies the ratio of the inertial and viscous terms in the momentum equation. In the non-rotating limit, dimensional analysis predicts the ‘ultimate’ scaling \[Re\propto Ra^{\frac{1}{2}}\label{eq:ReMLT}\tag{10}\] under the assumption that the velocity \(U\) (used to compute \(Re\)) and the temperature difference \(\Delta T\) (used to compute \(Ra\)) are independent of the diffusivities. When instead the system is rotating, then a diffusion-free balance between the Coriolis, inertial, and Archimedean (buoyancy driving) terms in the momentum equation (CIA balance) suggests the scaling [5], [11], [38]: \[Re\propto \left(RaNuPr^{-2}Ta^{-\frac{1}{4}}\right)^{\frac{2}{5}},\label{eq:ReCIA}\tag{11}\] sometimes referred to as the inertial scaling of rotating convection. If alternatively, a balance is sought between the Coriolis, viscous, and Archimedean terms (i.e., VAC balance; [39]), then the following scaling: \[Re\propto\left(RaNuPr^{-2}Ta^{-\frac{1}{3}}\right)^{\frac{1}{2}}\label{eq:ReVAC}\tag{12}\] is obtained. Equation 12 implies that \(U\) depends on the diffusivity. We are not aware of any study that recovers either of the diffusion-free scalings (Equations 10 and 11 ) from simulations of spherical shell convection that include rotation. For example, even the simulations presented in [11], that achieve an RMLT-like scaling for \(Nu\) in the rapidly rotating regime, exhibit a Reynolds number scaling that is not diffusion-free. There the rapidly rotating cases follow an alternate scaling in which CIA balance is obtained in the bulk, but kinetic energy dissipation is dominated by the boundary layers, while their weakly-nonlinear simulations follow the VAC scaling described above. However, we note that [12] do find diffusion-free behaviour for the kinetic energy in simulations that are non-rotating.

Figure 6c and 6d present a comparison between the Reynolds number scaling obtained in our simulations, and the theoretical relations given by Equations 10 –12 . As with the Nusselt number, the three regimes we have identified correspond to three distinct scaling relations. Experiments with \(Ro_{\text{cv,F}}>1\) in the rotationally-uninfluenced regime, for which small scales are rotationally-unconstrained, follow the ultimate scaling given by Equation 10 (green points in Figure 6d). Meanwhile, the experiments that are rotationally-influenced (\(Ro_{\text{bulk}}<1\) but \(Ro_{\omega}\gtrsim1\)) follow the inertial scaling (blue points in Figure 6d). By contrast, the rotationally-constrained \(Ro_{\text{cv,F}}<0.1\) cases display a Reynolds number scaling that closely follows the diffusivity-dependent VAC balance scaling (orange points in Figure 6c), which is consistent with the dependence of the KE on the diffusivities in Figure 2.

4 Summary and Discussion↩︎

We have presented 3D simulations of rotating convection in a spherical shell, driven by an internal heating and cooling function. This function was designed so that the same flux is extracted by cooling (over a fixed region near the top of the computational domain) as is input by heating (over an identically-sized region near the base). A key finding of our work is that in some regimes, the heat transport and other aspects of the flow dynamics are well described by diffusion-free theory, with no evident dependence on the viscosity or thermal diffusivity. Below, we briefly summarise the main results of this study, relate them to prior work, and comment on their astrophysical implications.

The scaling results presented in the previous section suggest that our simulations can be grouped into three distinct regimes, which we term rapidly rotating (\(Ro_{\text{cv,F}}<0.1\)), rotationally-influenced (\(0.1\leq Ro_{\text{cv,F}}<1\)), and rotationally-uninfluenced (\(Ro_{\text{cv,F}}\geq1\)). In the rotationally-constrained regime, both the small scales and the bulk are constrained by rotation. The rotationally-influenced regime is an intermediate regime where large scales are constrained by rotation, but small scales are not (Figure 1d). Finally, there is no discernible impact of rotation on the dynamics of simulations in the rotationally-uninfluenced regime. The key features of each regime are summarised in Table [tab:1] and described below.

The main result of this paper is that all of our experiments (aside from those that are very weakly supercritical) exhibit Nusselt number scalings and a radial temperature contrast that are independent of the diffusivities (Figure 2a, Figure 6a,b). Simulations in the non-rotating and rotationally-influenced regimes follow scalings that are consistent with MLT and RMLT, respectively; the rotationally-constrained cases also follow a scaling that is RMLT-like, but offset by a factor of the (diffusion-free) convective Rossby number (\(Ro_{\text{cv,F}}\)). In addition, diffusion-free scaling is also obtained for the Reynolds number in the non-rotating and rotationally-influenced regimes (Figure 6d), but in the rotationally-constrained regime the experiments follow a Reynolds number scaling consistent with VAC balance (which depends on the diffusivity; Figure 6c). Finally, in the rotationally-constrained regime, the morphology of the differential rotation depends on the supercriticality; weakly supercritical simulations (\(RaTa^{-2/3}\lesssim4\)) feature Solar-like differential rotation (denoted FS), whereas moderately supercritical simulations additionally feature a strong, prograde polar vortex (denoted FSF) (as in Figure 5, bottom row). In the rotationally-influenced and non-rotating regimes, anti-Solar differential rotation (denoted SF) is obtained (as in Figure 5, top row), with a structure and amplitude that is independent of the diffusivity.

An important result from our simulations is that it is possible to obtain zonally-averaged zonal flow statistics that are independent of the diffusivity. This result is demonstrated by rotationally-influenced cases with finite amplitude anti-Solar differential rotation (Figure 2b and Figure 5). In the rotationally-constrained regime, the differential rotation does depend on the diffusivity, although we suggest that too few of our experiments have sufficient supercriticality (only six experiments with \(RaTa^{-2/3}>10\) and \(Ro_{\text{cv,F}}<0.1\)) to explore the possibility for diffusion-free prograde differential rotation. Notwithstanding the diffusivity-dependence of the differential rotation in the rotationally-constrained regime, we have found that the sign of the equatorial zonal flow (i.e., prograde – FS or FSF, or retrograde – SF), does not depend on the diffusivity.

Previous work has shown that ‘local Nusselt numbers’ defined for the equatorial and polar regions can yield different scaling results that are set by the local dynamics [37]. Therefore, the fact that RMLT does not fully capture the convective heat transport scaling in the rapidly-rotating experiments is not surprising, given that the morphology of the convection depends strongly on latitude (e.g., wave-like dynamics at low latitudes but convective turbulence at high latitudes) but this effect is not accounted for by RMLT. A full analysis of the latitudinal heat budget, as well as the dependence of the size of the equatorial region on \(Ro_{\text{cv,F}}\), is required to form a complete understanding of the \(Nu\) scaling in the rapidly-rotating regime. A second, puzzling feature of the rapidly-rotating regime is that the convective velocity scaling (i.e., the \(Re\) scaling) is diffusivity-dependent, but the heat transport scaling is not. Our analysis in this article has focused on the CZ in isolation; however, by analogy with the influence that boundary layers can have on the bulk in boundary-driven convection [11], [40], [41], it is plausible that the dynamics of the heated and cooled regions can exert an influence on the CZ. We note that forming a volume-integrated relation between \(Re\) and the viscous dissipation rate is one route to obtaining the inertial scaling for \(Re\) (Equation 10 ; [11]). As we are interested in the dynamics integrated over the CZ, this integral relation will not be exact, but instead contain ‘boundary terms’ that describe the exchange of kinetic energy between the CZ and the heated and cooled regions. Understanding the details of this exchange are likely important to understanding the scaling behaviour of \(Re\) in the CZ.

Several aspects of the flows here are unlike those realised in a star like the Sun. We believe most of these discrepancies result from choices and simplifications made in our modeling (rather than from a direct dependence on diffusivity), but testing this assertion will require future work. For example, none of the zonal flows in our calculations are particularly ‘Solar-like’; those that exhibit a fast equator typically also have a ‘polar vortex’, consisting of retrograde flow at high latitudes. This may reflect the narrow aspect ratio of the convective shell considered here (see, e.g., [42]). More generally, we have restricted ourselves to the simplest case of Boussinesq, hydrodynamic convection. This set-up has the advantage of being easily relatable to theory. However, both magnetism and compressibility have important effects on the dynamics of rotating, spherical shell convection (e.g., [9], [43]–[47]). Therefore, incorporating these effects is an obvious and important follow-up to this work. In addition, we have not analysed the sensitivity of our results to the functional form of the prescribed internal heating and cooling profile. We intend to study these effects in future work.

In spite of these limitations, we argue that our findings have very promising implications for the numerical modeling of convection in stellar (and giant planet) interiors. All simulations of interior convection operate with diffusivities (whether explicit or numerical) that are orders of magnitude larger than those of real objects, so it is vital to understand how these diffusivities alter the results of the simulations. Our demonstration here that aspects of the heat transfer and flow dynamics are diffusion-free (i.e., independent of the diffusivities) will, we believe, enable more confident extrapolation from our results to real environments. To be explicit, we think that if a simulation in the rotationally-uninfluenced or rotationally-influenced regimes could be conducted with vastly lower diffusivities (but at the same value of \(Ro_{cv, F}\), and in a otherwise identical setup), its temperature gradient, integrated fluctuating kinetic energy, and differential rotation would not differ from those achieved in our simulations at high supercriticality. Moreover, we believe that some aspects of the rotationally-constrained cases are approaching a diffusion-free state for which a similar conclusion could be drawn, but simulations with a higher supercriticality are required to demonstrate this conclusively.

5 Simulation Input Parameters↩︎

References↩︎