On differential rotation and overshooting in solar-like stars

TL;DR

This study uses numerical simulations to explore how changes in rotation rate affect the dynamics, differential rotation states, and overshooting in the convective envelopes of solar-like stars across spectral types, revealing multiple rotation regimes.

Contribution

It introduces a comprehensive analysis of differential rotation states and overshooting in stars with varying rotation rates, using 3D simulations and scaling laws to map stellar behavior.

Findings

Identification of three main rotation states: anti-solar, solar-like, and cylindrical.

Discovery of a transition to Jupiter-like rotation profiles under strict rotational constraints.

Proposal of a Rossby number threshold for the existence of anti-solar differential rotation.

Abstract

We seek to characterize how the change of global rotation rate influences the overall dynamics and large scale flows arising in the convective envelopes of stars covering stellar spectral types from early G to late K. We do so through numerical simulations with the ASH code, where we consider stellar convective envelopes coupled to a radiative interior with various global properties. As solar-like stars spin down over the course of their main sequence evolution, such change must have a direct impact on their dynamics and rotation state. We indeed find that three main states of rotation may exist for a given star: anti-solar-like (fast poles, slow equator), solar-like (fast equator, slow poles), or a cylindrical rotation profile. Under increasingly strict rotational constraints, the latter profile can further evolve into a Jupiter-like profile, with alternating prograde and retrograde zonal jets. We have further assessed how far the convection and meridional flows overshoot into the radiative zone and investigated the morphology of the established tachocline. Using simple mixing length arguments, we are able to construct a scaling of the fluid Rossby number $R_{of} = \tilde{\omega}/2\Omega_* \sim \tilde{v}/2\Omega_* R_*$, which we calibrate based on our 3-D ASH simulations. We can use this scaling to map the behavior of differential rotation versus the global parameters of stellar mass and rotation rate. Finally, we isolate a region on this map ($R_{of} \gtrsim 1.5-2$) where we posit that stars with an anti-solar differential rotation may exist in order to encourage observers to hunt for such targets.

Tables

Table 1. Table 1 Global properties on the main sequence of the 4 stars used in our ASH models

Mass	Radius	$L_{\ast }$	$T_{eff}$	Sp. T.	$M_{bcz}$	$R_{bcz}$	${\overline{T}}_{bcz}$	${\overline{\rho }}_{bcz}$	${\Delta }_{cz}\overline{\rho }$	${\Delta }_f\overline{\rho }$
$\left(M_{\odot }\right)$	$\left(R_{\odot }\right)$	$\left(L_{\odot }\right)$	$\left(K\right)$		$(M_{\odot },M_{\ast })$	$(R_{\odot },R_{\ast })$	$\left(K\right)$	$\left(gc{m}^{-3}\right)$	-	-
0.5	0.44	0.046	4030	K7	0.18, 0.36	0.25,0.56	$4.3\times {10}^6$	14.0	42	193
0.7	0.64	0.15	4500	K4/K5	0.079, 0.11	0.42,0.66	$3.0\times {10}^6$	2.1	50	605
0.9	0.85	0.55	5390	G8	0.042, 0.046	0.59,0.69	$2.6\times {10}^6$	0.51	67	1013
1.1	1.23	1.79	6030	G0	0.011, 0.010	0.92,0.75	$1.6\times {10}^6$	0.048	81	830

Table 2. Table 2 Parameters for the Stellar ASH Models

Mass	$R_{bot}$	$R_{top}$	$R_{cz}$	Name	Rotation ${\Omega }_{\ast }$
$\left(M_{\odot }\right)$	$\left(R_{\ast }\right)$	$\left(R_{\ast }\right)$	$\left(R_{\ast }\right)$		$\left({\Omega }_{\odot }\right)$
0.5	0.13	0.95	0.56	M05 S	0.125
				M05 R1	1
				M05 R3	3
				M05 R5	5
0.7	0.32	0.97	0.66	M07 S	0.3
				M07 R1	1
				M07 R3	3
				M07 R5	5
0.9	0.38	0.97	0.69	M09 S	0.5
				M09 R1	1
				M09 R3	3
				M09 R5	5
1.1	0.50	0.97	0.75	M11 R1	1
				M11 R3	3
				M11 R5	5

Table 3. Table 3 Representative Velocities

	Mid Convective Zone
Case	${\tilde{v}}_r$	${\tilde{v}}_{\theta }$	${\tilde{v}}_{\mathrm{pol}}$	${\tilde{v}}_{\mathrm{mc}}$	${\tilde{v}}_{\varphi }$	${\tilde{v}}_{\varphi }^{\prime }$	$\tilde{v}$	${\tilde{v}}^{\prime }$
M05 S	15	14	20	3.3	25	18	32	27
M05 R1	8	8	11	1.1	40	8	41	14
M05 R3	6	7	9	0.5	53	7	54	11
M05R5	5	5	7	0.3	46	5	46	9
M07 S	27	26	37	3.2	43	36	57	50
M07 R1	22	21	30	1.3	77	24	83	39
M07 R3	17	18	25	0.9	116	23	119	34
M07 R5	15	16	22	0.3	131	20	133	30
M09 S	61	47	77	4.7	70	63	105	99
M09 R1	50	46	68	2.5	100	52	121	85
M09 R3	44	44	62	1.4	240	53	248	82
M09 R5	27	38	47	1.9	259	47	263	66
M11 R1	149	133	200	22.2	225	133	300	238
M11 R3	103	96	141	3.6	404	107	428	176
M11 R5	83	82	117	3.0	605	93	616	149

Table 4. Table 4 Non-dimensional numbers and Stellar, Convective and Fluid Rossby numbers

Name	$R_e$	$R_a$	$T_a$	$P_e$	$R_{of}$	$R_{oc}$	$R_{os}$
		$\left({10}^6\right)$	$\left({10}^6\right)$
M05 S	107	0.80	0.09	14.8	1.77	1.32	2.44
M05 R1	131	9.98	33.5	18.7	0.35	0.82	0.16
M05 R3	179	54.1	907.7	24.4	0.16	0.64	0.04
M05 R5	189	97.7	4202.5	26.2	0.09	0.51	0.02
M07 S	58	0.12	0.041	7.9	1.23	0.73	1.94
M07 R1	72	0.96	1.66	10.2	0.42	0.62	0.39
M07 R3	109	5.81	44.9	13.6	0.17	0.51	0.10
M07 R5	124	27.0	208.1	15.5	0.11	0.66	0.05
M09 S	60	0.13	0.07	9.3	1.29	0.82	1.79
M09 R1	64	0.43	0.41	9.4	0.67	0.74	0.73
M09 R3	106	2.54	10.9	14.2	0.28	0.60	0.21
M09 R5	110	5.48	50.9	11.2	0.14	0.53	0.08
M11 R1	63	0.18	0.08	9.9	1.40	2.45	1.80
M11 R3	81	1.13	2.1	11.9	0.54	2.05	0.41
M11 R5	89	2.86	9.6	12.4	0.34	1.96	0.20

Table 5. Table 5 Kinetic energies normalized by the shell volume, Latitudinal angular velocity, Temperature and Entropy latitudinal contrasts

Name	KE	DRKE	MCKE	CKE	$\Delta \Omega$	$\Delta T$	$\Delta S$
	${10}^{6}erg$	${10}^{6}erg$	${10}^{3}erg$	${10}^{6}erg$	$nHz$	$K$	${10}^{3}erg/g/K$
M05 S	16.0	3.07 (19.2%)	241.4 (1.5%)	12.7 (79.3%)	-24	-1.9	-0.005
M05 R1	31.2	27.9 (89.4%)	24.4 (0.1%)	3.28 (10.5%)	129	2.9	0.2
M05 R3	78.8	74.8 (94.9%)	8.92 (0.01%)	4.01 (5.09%)	85	2.8	0.2
M05 R5	58.4	55.9 (95.7%)	5.80 (0.01%)	2.48 (4.29%)	146	12.9	0.9
M07 S	3.59	0.775 (21.5%)	56.4 (1.6%)	2.76 (76.9%)	-32	-2.8	-0.01
M07 R1	8.89	7.17 (80.6%)	13.8 (0.2%)	1.71 (19.2%)	120	2.5	0.2
M07 R3	21.9	20.2 (91.7%)	6.85 (0.1%)	1.81 (8.2%)	187	13.8	1.1
M07 R5	34.9	30.9 (88.7%)	5.05 (0.01%)	3.94 (11.29%)	223	27.8	2.3
M09 S	3.25	0.519 (16.0%)	32.6 (1.0%)	2.69 (83.0%)	-25	-2.5	-0.1
M09 R1	4.59	2.58 (56.2%)	16.9 (0.4%)	1.99 (43.4%)	108	8.4	0.9
M09 R3	19.9	18.4 (92.1%)	9.61 (0.05%)	1.57 (7.85%)	288	57.5	5.9
M09 R5	5.49	28.2 (51.5%)	7.31 (0.01%)	26.6 (48.49%)	338	94.9	9.75
M11 R1	3.08	1.22 (39.6%)	25.5 (0.8%)	1.84 (59.6%)	-131	-10.2	-1.52
M11 R3	9.96	8.79 (88.3%)	6.89 (0.07%)	1.16 (11.63%)	291	133.5	20.7
M11 R5	22.8	21.7 (94.9%)	4.61 (0.02%)	1.16 (5.08%)	435	328.9	52.1

Table 6. Table 6 Kinetic energy exchanges

	$Q_V$		$Q_C^r$	$Q_C^l$	$Q_R^r$	$Q_R^l$
	$\left[L_{\star }\right]$	$\left[L_{\odot }\right]$	$\left[L_{\star }\right]$	$\left[L_{\star }\right]$	$\left[L_{\star }\right]$	$\left[L_{\star }\right]$
M05 S	-0.053	-0.002	0.065	-0.012	0.096	-0.104
M05 R1	-0.163	-0.007	0.056	-0.021	0.149	0.023
M05 R3	-0.195	-0.009	0.108	-0.094	0.184	0.055
M05 R5	-0.229	-0.011	0.095	-0.074	0.157	0.046
M07 S	-0.042	-0.006	0.019	-0.016	0.067	-0.017
M07 R1	-0.129	-0.020	0.002	0.019	0.094	0.024
M07 R3	-0.170	-0.026	0.020	0.010	0.118	0.037
M07 R5	-0.184	-0.028	0.004	-0.023	0.138	0.026
M09 S	-0.029	-0.016	0.006	-0.004	0.043	-0.003
M09 R1	-0.063	-0.035	0.001	0.011	0.052	0.016
M09 R3	-0.170	-0.093	0.008	0.026	0.113	0.025
M09 R5	-0.155	-0.085	0.042	0.025	0.056	0.020
M11 R1	-0.035	-0.063	0.033	-0.014	0.040	-0.018
M11 R3	-0.093	-0.166	0.011	0.006	0.077	0.015
M11 R5	-0.119	-0.213	0.039	-0.010	0.098	0.022

Table 7. Table 7 Diffusivity profile and nuclear heating source parameters

Mass	Name	${\nu }_{top}$	$r_t$	${\sigma }_t$	${ϵ}_0$	$n_c$	a	b
$\left(M_{\odot }\right)$		$\left(c{m}^2{s}^{-1}\right)$	$\left(cm\right)$	$\left(cm\right)$
0.5	M05 S	$10.5\times {10}^{11}$	$1.65\times {10}^{10}$	$4.0\times {10}^8$	$2.20\times {10}^{-7}$	7.3	9.78e-04	-3.57e-09
	M05 R1	$4.47\times {10}^{11}$
	M05 R3	$2.58\times {10}^{11}$
	M05 R5	$2.00\times {10}^{11}$
0.7	M07 S	$5.31\times {10}^{12}$	$2.67\times {10}^{10}$	$1.0\times {10}^9$	$5.56\times {10}^{-9}$	8.8	9.59e-03	-9.56e-09
	M07 R1	$2.91\times {10}^{12}$
	M07 R3	$1.68\times {10}^{12}$
	M07 R5	$1.30\times {10}^{12}$
0.9	M09 S	$1.43\times {10}^{13}$	$3.94\times {10}^{10}$	$7.0\times {10}^8$			1.37e-02	-3.78e-08
	M09 R1	$1.01\times {10}^{13}$
	M09 R3	$5.81\times {10}^{12}$
	M09 R5	$4.50\times {10}^{12}$
1.1	M11 R1	$3.80\times {10}^{13}$	$6.20\times {10}^{10}$	$8.0\times {10}^8$			1.09e-02	-3.01e-07
	M11 R3	$2.20\times {10}^{13}$
	M11 R5	$1.70\times {10}^{13}$

Equations

$$v=c_1{\left(\frac{L_{\ast }}{{\rho }_{bcz}{R}_{\ast }^2}\right)}^{1/3}$$

$$\displaystyle\mbox{\boldmath$\nabla$}\cdot(\bar{\rho}{\bf v})$$

$$\displaystyle\frac{\partial{\bf v}}{\partial t}+({\bf v}\cdot\mbox{\boldmath$\nabla$}){\bf v}$$

$$\displaystyle-\frac{1}{\bar{\rho}}\mbox{\boldmath$\nabla$}\cdot\mbox{\boldmath$\cal D$}$$

$$\bar{\rho }\bar{T}\frac{\partial S}{\partial t}$$

$$\displaystyle+\mbox{\boldmath$\nabla$}\cdot[\kappa_{r}\bar{\rho}c_{p}\mbox{\boldmath$\nabla$}(\bar{T}+T)$$

$$\displaystyle{\cal D}_{ij}=-2\bar{\rho}\nu[e_{ij}-1/3(\mbox{\boldmath$\nabla$}\cdot{\bf v})\delta_{ij}]\,,$$

$$\frac{\rho }{\bar{\rho }}=\frac{P}{\bar{P}}-\frac{T}{\bar{T}}=\frac{P}{\gamma \bar{P}}-\frac{S}{c_p},$$

$$\bar{P}=\mathcal{R}\bar{\rho }\bar{T},$$

$$\nu (r)={\nu }_{bot}+{\nu }_{top}{f}_{step}(r),$$

$$f_{step}(r)$$

$$f(r)$$

$$\beta$$

$$\displaystyle{\bar{\rho}\bf v}=\mbox{\boldmath$\nabla$}\times\mbox{\boldmath$\nabla$}\times(W\hat{\bf e}_{r})+\mbox{\boldmath$\nabla$}\times(Z\hat{\bf e}_{r})\,.$$

$$\Omega ={\Omega }_{\ast }+A{e}^{-\frac{1}{2}{\left[\frac{(x-x_c)}{w}\right]}^2},$$

$$\mathrm{DRKE}=\frac{1}{2}\int \bar{\rho }{⟨v_{\varphi }⟩}^2{r}^2\sin \theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi .$$

$$|\Delta \Omega |$$

$$\propto$$

$$|\Delta \Omega |$$

$$r\sin \theta ⟨\bar{\rho }{v}_r⟩=-\frac{1}{r}\frac{\partial \Psi }{\partial \theta }\text{mbox}andr\sin \theta ⟨\bar{\rho }{v}_{\theta }⟩=\frac{\partial \Psi }{\partial r}.$$

$$\left<\bar{\rho}{\bf v}_{m}\right>\mbox{\boldmath$\cdot$}\mbox{\boldmath$\nabla$}{\cal L}={\cal T}~{}~{}~{},$$

$$\mathcal{L}={\varpi }^2\Omega =\varpi \left(⟨v_{\varphi }⟩+\varpi {\Omega }_{\ast }\right),$$

$$\mathrm{MCKE}=\frac{1}{2}\int \bar{\rho }\left({⟨v_r⟩}^2+{⟨v_{\theta }⟩}^2\right)r^2\sin \theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi ,$$

$$\mathrm{MCKE}$$

$$\bar{\rho}\frac{\partial L}{\partial t}=\mbox{\bf{$\tau$}}$$

$$\mbox{\boldmath$\nabla\mbox{\boldmath$\cdot$}$}\mbox{\boldmath${\cal F}$}=\mbox{\boldmath$\nabla\mbox{\boldmath$\cdot$}$}\left({\cal F}_{r}\mbox{\boldmath$\hat{\bf e}_{r}$}+{\cal F}_{\theta}\mbox{\boldmath$\hat{\bf e}_{\theta}$}\right)=0$$

$${\mathcal{F}}_r$$

$${\mathcal{F}}_{\theta }$$

$$I_r(r)$$

$$I_{\theta }(\theta )$$

$$2{\Omega }_{\ast }\frac{\partial ⟨v_{\varphi }⟩}{\partial z}$$

$$\text{centering}\frac{\partial ⟨v_{\varphi }⟩}{\partial z}=\frac{g}{2{\Omega }_{\ast }r{c}_p}\frac{\partial ⟨S⟩}{\partial \theta }.\text{@add@centering}$$

$$\Delta S$$

$$\Delta T$$

$$\Delta S$$

Full text

1. Introduction

Most solar-like stars possess a deep turbulent convective envelope. Understanding their dynamical and nonlinear properties is fundamental in order to characterize the transport of heat, energy and angular momentum in these stars. Of particular importance are the large-scale flows, like the differential rotation and meridional circulation, that are achieved in these turbulent convective envelopes. Indeed, such flows are considered as essential ingredients to explain stellar magnetism (Weiss, 1994; Saar & Brandenburg, 1999; Jouve et al., 2010). For instance, the differential rotation can convert poloidal magnetic fields into toroidal fields via the so-called $\Omega$-effect (Parker, 1955a; Moffatt, 1978). Such toroidal magnetic fields are possibly at the origin of star spots (Parker, 1955b). Furthermore, the meridional circulation can transport both the surface poloidal field from low latitudes to the pole as well as the magnetic field at the base of the convective envelope in the so-called flux transport solar models (Wang et al., 1989; Schrijver, 2001; Charbonneau, 2005; Dikpati, 2011; Brun et al., 2015a), some of which consider self-consistent dynamo action.

Despite substantial efforts, there are still many unknowns related to how the large-scale flows and stellar convection change with global stellar parameters such as rotation $\Omega_{*}$ and mass $M_{*}$. Over the years, several observational campaigns have been pursued using either ground based or space-born instruments, to derive useful constraints on rotation of solar like stars (see for instance the recent studies by Bouvier, 2013; Reinhold & Reiners, 2013; do Nascimento et al., 2014; Garcia et al., 2014; Reinhold & Arlt, 2015, and references therein). What is clear from these observational studies is that young solar-like stars are fast rotators with rotational periods on the order of days and old stars are slow rotators, with periods closer to months. This trend of increasing rotation period, or decreasing rotation rate $\Omega_{*}$, is clearly seen on the main sequence. Indeed, as was first suggested by Skumanich (1972), it appears that $\Omega_{*}\propto t^{-1/2}$, with $t$ the stellar age.

We have a relatively good understanding about how the rotational braking of solar-like stars occurs. It is due to the continuous action of a stellar wind driven by thermal pressure taking away mass and angular momentum from the aging star (Parker, 1958; Schatzman, 1962; Weber & Davis, 1967; Kawaler, 1988; Matt et al., 2012). The relation between stellar age and rotation has been termed gyrochronology (Barnes, 2003, 2010; Meibom et al., 2015) and is a useful marker of stellar evolution. Note however that in their recent asteroseismic study with data coming from the Kepler satellite, van Saders et al. (2016) have questioned Skumanich’s law and the accuracy of gyrochronology for solar-like stars older than the Sun (e.g., $t>4.5Gyr$), but these techniques have their own accuracy issues (Aigrain et al., 2015). Associated with this decreasing influence of rotation is a lower level of stellar activity (Wilson, 1978; Noyes et al., 1984; Pizzolato et al., 2003; Reiners, 2012), also called magnetochronology (Vidotto et al., 2014; Folsom et al., 2016). There is thus a positive feedback loop between rotation, magnetism and wind in solar-like stars that is also important to characterize (Brown, 2014; Meibom et al., 2015; Matt et al., 2015; Réville et al., 2015).

It is to be expected that the global rotation rate of a star will influence its large-scale mean flows. How those flows vary with rotation rate is still poorly known. Observational studies of large-scale flows in solar-like stars are rare and even less accurate than studies of the stellar rotation rate itself. Most studies have focused on the surface differential rotation, for the meridional circulation is too weak to be detected. However, several interesting trends for the differential rotation in solar-like stars have been uncovered (Donahue et al., 1996; Barnes et al., 2005; Collier Cameron, 2007; Saar, 2009). First, the observed differential rotation $\Delta\Omega$ appears to increase with stellar mass, or more precisely with a higher effective temperature. Second, some variations of the surface differential rotation have been found with rotation rate: $\Delta\Omega\propto\Omega_{*}^{n}$, with $n$ a positive exponent. However, observers currently disagree on the exact value of the exponent $n$. Donahue et al. (1996), Messina & Guinan (2003), Saar (2009) advocate for a value close to 0.6-0.7, whereas Barnes et al. (2005) and Collier Cameron (2007) argue for a much smaller value around 0.15, implying only a weak dependence on $\Omega_{*}$. Studies based on asteroseismic data have more recently started providing trends that somewhat lie in between whose values ($n\sim 0.3$ in Reinhold & Gizon (2015); see also (Reinhold & Reiners, 2013; Garcia et al., 2014; Rauer et al., 2014; Reinhold & Gizon, 2015; Balona & Abedigamba, 2016)). Actually, Balona & Abedigamba (2016) even propose that the exponent $n$ is different for each spectral types, confirming a strong dependency of the differential rotation amplitude with effective temperature, massive stars having larger shear. Theoretical interpretation are thus in order to explain those seemingly different trends.

Another physical property of key interest in solar-like stars is the extent and amount of overshooting at the base of their convective envelope, where surrounding layers of extra mixing are triggered by convective plumes at the boundaries of convection zones (Roxburgh, 1978; Zahn, 1991). It is quite unknown how the overshoot of convection varies with global stellar parameters such as mass or rotation rate, yet it has a direct impact on mixing and transport of chemicals and magnetic fields. The amount and extent of overshooting is usually chosen to be proportional to the local pressure scale height and to be a small fraction of the convective heat flux (Maeder & Meynet, 1989; Zahn, 1991; Browning et al., 2004). Those model parameters are often calibrated using stellar isochrones to match observations of surface chemical abundances (see for instance Meynet et al. (1993)). However, it has been shown that the overshooting depth can be overestimated if the stellar models used to compute the isochrones are neglecting other physical processes such as rotational mixing (Ekström et al., 2012). Hence, we must find complementary means to constrain stellar overshooting.

It is extremely difficult to observationally infer information about the meridional circulation and overshooting layers in stars, and the observational data for surface differential rotation are often inadequately robust. Hence, theoretical developments and numerical simulations are useful to help characterize the large-scale mean flows and to assess the degree of coupling between convection and radiative zones. For instance, multi-cellular meridional flows, which are now confirmed to exist within the Sun using local helioseismic analysis (Haber et al., 2002; Mitra-Kraev & Thompson, 2007; Zhao et al., 2013), have been the natural outcome of solar convection 3-D simulations for many years (Miesch et al., 2000; Brun & Toomre, 2002; Aurnou et al., 2007; Käpylä et al., 2011; Gastine et al., 2013; Guerrero et al., 2013; Featherstone & Miesch, 2015), long before they became accepted as a plausible solar flow pattern by the community.

Thus, along side observational campaigns, several groups have attempted to make progress in understanding rotating convection and the establishment of large-scale flows in stars through numerical simulations. Many of these groups have benefited from the pioneering work of Gilman & Glatzmaier where they modeled the largest scales in the Sun’s convection zone (see the series of papers listed in Glatzmaier & Gilman, 1982). For such simulations, a full spherical geometry is required since we wish here to characterize a star’s global-scale flows. Such global-scale simulations of stellar convection have been performed using various numerical codes based on either finite volumes or on spectral methods (see Brun et al., 2015b, and references there in for a recent review). These various anelastic codes, such as ASH (Clune et al., 1999), have been benchmarked internationally (Jones et al., 2011). They generally agree on the main properties of global rotating convection when identical setups and nondimensional numbers are chosen. These codes have been used to simulate very different type of stars: from massive ones (Browning et al., 2004; Featherstone et al., 2009; Augustson et al., 2016) to dwarf stars (Dobler et al., 2006; Browning, 2008) and of course solar-like ones (Ballot et al., 2007; Brown et al., 2008; Bessolaz & Brun, 2011; Matt et al., 2011; Käpylä et al., 2011; Augustson et al., 2012; Guerrero et al., 2013; Fan & Fang, 2014; Karak et al., 2015). It is found that large-scale flows are sensitive to the intensity of the convective driving and to the rotation rate of the simulation. Yet almost none of these studies have systematically taken into account the coupling to a stably stratified interior (see Miesch et al., 2000; Guerrero et al., 2016, for counter-examples) and/or looked at the influence of the aspect ratio of the convection zone on the resulting convection and its mean flows, as we have simultaneously done in this work.

Another very important motivation for running global convection simulations is to understand stellar magnetism. This implies the need to compute non-ideal magnetohydrodynamic (MHD) simulations of magnetized convection. In particular, the focus of many of these simulation has been to find and understand how convective dynamo solutions can become cyclic and to further analyze them in terms of mean-field dynamo theory. Recently, significant progress has been made in that direction with many global MHD solutions now possessing a cyclic behavior. We refer to the following recent papers for a discussion of stellar magnetism and how it may arise (see for instance Brun et al., 2004; Browning et al., 2006; Ghizaru et al., 2010; Brown et al., 2010, 2011; Racine et al., 2011; Käpylä et al., 2012; Augustson et al., 2013; Käpylä et al., 2013; Nelson et al., 2013; Fan & Fang, 2014; Karak et al., 2015; Augustson et al., 2015; Lawson et al., 2015; Simitev et al., 2015; Guerrero et al., 2016). InVarela et al. (2016), we have started computing and studying the MHD equivalent to the 15 simulations here, which cover four mass bins and several rotation rates.

In this paper, we report on novel 3–D numerical experiments in spherical geometry designed to investigate how the complex, nonlinear dynamics occurring in the convective envelope of solar-like stars changes with stellar parameters such as mass and rotation rate. Our approach differs from other recent studies listed above by taking into account the dynamical influence of a deep and stable radiative interior. It is complementary to the study of Guerrero et al. (2013), since we further consider various aspect ratios and stellar spectral types. We also propose a model that utilizes a simple mixing length theory scaling to identify the possible states of differential rotation for various stellar spectral types. This model is then calibrated using the set of 15 models discussed below.

The paper is organized as follows. In §2, we derive a scaling relationship for the Rossby number based on mixing length arguments. In §3 we describe our equations, numerical models, and the various ingredients employed to model four stars of differing spectral type at a selection of rotation rates. In sections 4, 5, and 6, we discuss the properties of convection, penetration, and the large-scale flows of our models as well as derive scaling relationships for all key quantities. In §7, we perform a detailed analysis of angular, energy, and heat transport in our simulations. Finally, we conclude in §8.

2. Hints from mixing length theory for states of stellar differential rotation

It is worth noting that one can already guess by using simple mixing length scaling and 1-D stellar structure models the outcome of the simulations in terms of the overall surface differential rotation by distinguishing two states: fast vs slow equator rotation. Our 3-D numerical simulations are key as they will help us characterizing the differential rotation profile as a function of depth (and latitude) and how the coupling to a radiative interior may tilt the iso-contour of omega by so-called thermal wind effect (Miesch et al., 2006).

One can make an educated guess of the differential rotation state realized in a star, by evaluating the convective velocity from mixing length arguments (Kippenhahn & Weigert, 1994; Augustson et al., 2012; Brun et al., 2015a):

[TABLE]

with the typical values taken for the stellar luminosity, radius, and density at the base of the convection zone listed in Table 1 for the stellar spectral range considered in this study, and with $c_{1}$ a proportionality factor. Classical stellar evolution indicates that $L_{*}\sim M_{*}^{4}$ and $R_{*}\sim M_{*}^{0.9}$ (Kippenhahn & Weigert, 1994). Assuming that $\rho_{bcz}\sim M_{*}^{n}$, with $n<0$ but undetermined for now, one directly sees that $v\sim M_{*}^{(2.2-n)/3}$. Regression fits to the values listed in Table 2 allow us to obtain a scaling for $\rho_{bcz}$, and so we can refine the stellar mass dependence of $L_{*}$ and $R_{*}$. In particular, we find that $\rho_{bcz}\sim M_{*}^{-6.9}$, $L_{*}\sim M_{*}^{4.6}$, and $R_{*}\sim M_{*}^{1.3}$. Replacing these scalings in equation 1, yields $v\sim M_{*}^{3}$. As expected, more massive stars have faster convective flows. Knowing how $v$ is expected to scale with stellar mass, we can now compute an approximated fluid Rossby number $R_{of}=v/2\Omega_{*}R_{*}=c_{1}M_{*}^{1.7}/\Omega_{*}$ as a function of the stellar rotation rate $\Omega_{*}$. We choose to use the fluid Rossby number instead of the stellar or convective ones. We defer the reader to Appendix  B for a further discussion of the various definitions of Rossby numbers.

In Figure 1, we display the fluid Rossby number as a function of stellar rotation rate and mass. One can easily identify 2 regimes separated by the Rossby equals 1 white line: small Rossby number ($\lesssim 0.66$) and large Rossby number ($\gtrsim 1.5$). From the various studies published over the last decade (Ballot et al., 2007; Brown et al., 2008; Bessolaz & Brun, 2011; Augustson et al., 2012; Takehiro et al., 2013; Gastine et al., 2014b; Käpylä et al., 2014), it is clear that simulations of rotating convective shells with small Rossby numbers often possess fast prograde flows at the equator and slowly rotating poles, possibly ressembling a solar-like differential rotation. In contrast, those simulations with a large Rossby number typically have retrograde flows at the equator and fast rotatin polar regions e.g. something like an anti-solar-like differential rotation. With $0.66\lesssim R_{of}\lesssim 1.5$, it is more difficult to anticipate the result given the crudeness of our derivation, but one may expect that for $R_{of}\gtrsim 1$ the simulation will likely be slowly rotating at the equator. In the remainder of the paper we will use the following terminology for characterizing differential rotation profile: solar-like will mean fast equator, slow poles, anti-solar like will mean slow equator and fast poles, cylindrical, will mean that the iso-contours of $\Omega(r,\theta)$ are constant along cyclinders aligned with the rotation axis and Jupiter-like differential rotation will mean that the profile is cylindrical but non monotonic, with alternance of prograde and retrogade jets. Analyzing this figure further, one might expect that a 1.1 solar mass star could be anti-solar-like at a rotation rate around the solar rate. Such a state may also occur for respectively 0.9, 0.7 and 0.5 solar mass stars at rotation rates around 0.6, 0.4, and 0.25 times the solar rate. One also sees that for a solar-mass star the solar rotation rate corresponds to a Rossby number less than unity and hence it likely has prograde equatorially, as observed. In the following, we will characterize further the profile of differential rotation and compare the outcome of nonlinear 3-D numerical simulations of rotating convection with this simple analysis based on the mixing length. We will in particular show that for low values of the Rossby number, there is a shift from a conical to cylindrical profile, where this third state is akin to Jupiter’s alternating zonal jets.

3. Modelling stars in 3-D with ASH

We present the simulation setup used to model the various spectral type stars considered in our study with the ASH code.

3.1. Model Equations

We use the ASH code (see Clune et al., 1999; Miesch et al., 2000; Brun et al., 2004) to model solar-like stars with mass ranging from 0.5 to 1.1 $M_{\odot}$. The domain of each of these simulations is large enough to encompass both a portion of the deep radiative zone and most of the overlying convective envelope that is representative of the targeted star. Hence, we are self-consistently capturing the nonlinear interactions that seamlessly couple those two zones. ASH solves the full set of 3–D anelastic equations of motion in a rotating, convective and radiative spherical shell utilizing massively-parallel computing architectures. These equations are fully nonlinear in the velocity variable. However, under the anelastic approximation, the thermodynamic variables are linearized with respect to a spherically symmetric and evolving mean state having a density $\bar{\rho}$, pressure $\bar{P}$, temperature $\bar{T}$ and specific entropy $\bar{S}$. Fluctuations about this reference state are denoted by $\rho$, $P$, $T$, and $S$. The resulting equations are (Glatzmaier, 1984; Clune et al., 1999):

[TABLE]

where ${\bf v}=(v_{r},v_{\theta},v_{\phi})$ is the local velocity in spherical coordinates in the frame rotating at constant angular velocity ${\bf\Omega_{*}}$, ${\bf g}$ is the gravitational acceleration, $c_{p}$ is the specific heat per unit mass at constant pressure, $\varpi=P/\bar{\rho}$ is the reduced or kinematic pressure, $\kappa_{r}$ is the radiative diffusivity, and ${\bf\cal D}$ is the viscous stress tensor. The components of ${\bf\cal D}$ are given by

[TABLE]

where $e_{ij}$ is the strain rate tensor, and $\nu$, $\kappa$ and $\kappa_{0}$ are effective eddy diffusivities. A volumetric heating term $\bar{\rho}\epsilon$ is also taken into account to mimic generation of energy by nuclear reactions. The nuclear reactions are modelled very simply by assuming that $\epsilon=\epsilon_{0}\bar{T}^{n_{c}}$. By enforcing that the integrated luminosity of the star match its known surface value, we can determine $\epsilon_{0}$ and $n_{c}$ as listed in Table 7. Note that only M05 and M07 series of models require that heating source term since their computational domain include a portion of the nuclear energy generation core. Here we solve the energy conserving anelastic equations, which have a momentum equation (3) that takes a slightly different form than the compressible Navier-Stokes equations. These anelastic equations have been shown to properly conserve energy in both convection zones and stably-stratified regions like the tachoclines that we study here (Brown et al., 2012; Vasil et al., 2013).

To complete the set of equations, we use the linearized equation of state

[TABLE]

where $\gamma$ is the adiabatic exponent, and assume the ideal gas law

[TABLE]

where ${\cal R}$ is the gas constant. The reference state is derived from a 1–D solar structure model (Brun et al., 2002, cf. §3.2) and is continuously updated with the spherically-symmetric components of the thermodynamic fluctuations as the simulation proceeds. It begins in hydrostatic balance so the bracketed term on the right-hand-side of equation (3) initially vanishes. However, as the simulation evolves, turbulent pressure drives the reference state slightly away from hydrostatic balance.

Due to limitations in computing resources, no simulation achievable now or in the near future can hope to directly capture all scales of stellar convection from global to molecular dissipation scales. The simulations reported here resolve nonlinear interactions among a large range of scales both in the convective and radiative zones. The nonlinear coupling of the two zones plus the use of a realistic stratification in the radiative interior is what sets these 3-D global simulations of solar-like stars apart from previous work. Motions and waves must exist in the solar-like stars on scales smaller than our grid resolution. In this sense, our models should be regarded as large-eddy simulations (LES) with parameterizations to account for subgrid-scale (SGS) motions.

Thus the effective eddy diffusivities $\nu$, $\kappa$ and $\kappa_{0}$ represent momentum and heat transport by motions which are not resolved by the simulation. They are allowed to vary with radius but are independent of latitude and longitude, and vary only slightly with time for a given simulation as the reference density evolves. Their amplitudes and radial profiles are varied depending on the resolution and objectives of each simulation. In the simulations reported on here, the radial profiles of $\nu$ and $\kappa$ are given by

[TABLE]

where

[TABLE]

and with $\nu_{top}$ in cm2 s*-1* and $r_{t}$ and $\sigma_{t}$ in cm are given in Table 7 (see Appendix), $\alpha$ is -0.5 for all cases. All models assumed a Prandtl number of 0.25, such that $\kappa$ can be directly obtained from the amplitude and profile of $\nu$. These tapered profiles are chosen in order to take into account the much smaller sub-grid scale transport expected in the convectively stable radiative interior. Their profile is shown in Figure 2.

To maintain a high degree of supercriticality of the convective instability in our simulations, we have lowered the diffusivities as we increase the rotation rate (see Table 7). We have chosen to scale $\nu_{top}\propto 1/\Omega_{*}^{0.5}$. This is a compromise between keeping the diffusivity constant but making the convective patterns too laminar and the exact scaling as $\propto 1/\Omega_{*}^{2}$ which would otherwise implies too wide a parameter range to cover given our computer resources.

The diffusivity $\kappa_{0}$ is set such as to have the unresolved eddy flux carrying the stellar flux (which depends on the spectral type considered see Table 1) outward at the top of the domain (see Figure 6). It drops off exponentially with depth in order to avoid a large inward heat flux in the stable zone (see Miesch et al., 2000; Brun et al., 2011; Alvan et al., 2014). Of course there is some arbitrariness in choosing the exact shape and amplitude of our diffusivity profiles and we optimize their profiles such as to limit their influence on the results reported here.

The velocity and thermodynamic variables are expanded in spherical harmonics $Y_{\ell m}(\theta,\phi)$ for their horizontal structure and in Chebyshev polynomials $T_{n}(r)$ for their radial structure (see Glatzmaier, 1984; Clune et al., 1999, for more details on the numerical method and anelastic approximation).

Given that the convective and radiative zones are nonlinearly and dynamically coupled, internal waves can easily be excited by the pummeling of convective plumes on top of the radiative interior (Brun et al., 2011; Alvan et al., 2014). The Brunt-Väisälä frequency $N$ of the models are very close to that deduced from 1–D stellar models computed with the CESAM code (Morel, 1997; Brun et al., 2002, see Figure 4). The transition at the base of the convective envelope has just been slightly soften, as can be seen in the Figure 4 when comparing solid and dash lines. We are thus expecting the propagation of the internal waves to be realistic, aside from the shallower cavity due to our choice of $r_{bot}\neq 0$ and the enhanced thermal and viscous diffusion present in the model that translates into an enhanced damping (Zahn et al., 1997). While present in the simulations, we will not study in details the internal waves and their spectra in this paper. We choose instead to focus our study on convection and the generation and maintenance of the large-scale mean flows and how they vary as a function of spectral type. As shown in Alvan et al. (2014), fully spherical models are more adequate to realistically model internal waves and gravity modes, but they are much more expensive to run, which makes a comprehensive parameter study impractical.

In order to ensure that the mass flux remains divergenceless, we use a toroidal–poloidal decomposition as:

[TABLE]

This system of hydrodynamic equations requires 8 boundary conditions in order to be well-posed. Since assessing the angular momentum redistribution in our simulations is one of the main goals of this work, we have opted for torque-free velocity conditions:

impenetrable top and bottom: $v_{r}=0|_{r=r_{bot},r_{top}}$ 2. 2.

stress free top and bottom:

$\frac{\partial}{\partial r}\left(\frac{v_{\theta}}{r}\right)=\frac{\partial}{\partial r}\left(\frac{v_{\phi}}{r}\right)=0|_{r=r_{bot},r_{top}}$ 3. 3.

constant entropy gradient at top and bottom: $\frac{\partial\bar{S}}{\partial r}=a|_{r=r_{bot}}\mbox{ and }b|_{r=r_{top}}$

the values of $a$ and $b$ depend on the modelled star (see Table 7).

3.2. Numerical Experiments

Our numerical model is a simplified portrayal of convection and radiative zones in solar like-stars: typical values deduced from 1-D stellar evolution models are taken for the heat flux, rotation rate, mass and radius, and a perfect gas is assumed (Morel, 1997). The anelastic reference state is based on a 1–D standard stellar structure model discussed in Brun et al. (2002). We list the main stellar parameters in Table 1. We also show in Figures 3 and 4 the density, pressure, temperature and entropy gradient radial profiles used for the four stellar spectral types studied. We initialize the reference state of the 3–D model by specifying the entropy gradient $d\bar{S}/dr$ and gravitational acceleration $g$ based on the 1–D model. The steep negative entropy gradient near the stellar surface is artificially suppressed to avoid the driving of small-scale convective motions that cannot be resolved in this model. We then solve the equation of hydrostatic balance for the reference density $\bar{\rho}$ with a Newton-Raphson method, assuming an ideal gas equation of state and using the density profile of the 1–D structure model as an initial guess to initiate the iterative solve. The resulting reference state is close to the stellar structure model, with slight departures due to the modified entropy profile in the convection zone and the ideal gas equation of state. Similarly, the radial profile of the radiative diffusivity $\kappa_{r}$ is based on the 1–D stellar structure model, slightly adjusted to accommodate the small changes between the reference state and the structure model and later on to compensate for the inward enthalpy flux at the base of the convection zone. The radial profiles of $\kappa_{r}$ are shown in Figure 2 for all four stellar models.

Thus the departure of the reference entropy gradient from the stellar structure model near the top of the convection zone evident in Figure 4b is largely imposed. However, the departure near the base of the convection zone is established by the convection itself, as downflow plumes deposit low entropy material just before entering the stable radiative zone. The initial reference state in this region follows more closely the 1–D structure model.

The computational domain for each of the 4 stellar masses considered is listed in Table 7 and shown in Figure 5. All models use a numerical resolution of $N_{r}\times N_{\theta}\times N_{\phi}=770\times 256\times 512$. The depth of the convection zone $r_{bcz}$ is defined by the change of sign of the initial mean entropy gradient $d\bar{S}/dr$, and it is listed in Table 2. This depth is slightly modified by the convective motions as the simulation evolves and matures. We also list the density contrasts both in the convective envelope and over the whole radial extent of the domain. We note that several scale heights are present in the convective envelope and that overall the models have large density contrasts. The resolution at the base of the convection zone is $\delta r=1.5\times 10^{-3}\,R_{*}$.

4. Convection in spherical shells of various aspect ratios and rotation rates

Once the 1-D structure of the model is established, a 3-D random perturbation of the entropy field is introduced such as to trigger the convective instability. Our simulations begin with a Rayleigh number of the order of $\sim 10^{6}$ to $10^{7}$, which is substantially larger than the critical Rayleigh number that is typically about $10^{4}$ for the value of the Taylor numbers used here (see Jones et al., 2009). Following a linear phase of exponential growth, the convective instability non linearly saturates by reducing the entropy gradient in the bulk of the domain except for intense thermal boundary layers created near the surface and at the base of the convective envelope. After several convective overturning times the radial transport of energy reaches an equilibrium as shown in Figure 6. As can be seen, the enthalpy flux, due to the correlation between radial velocity and temperature fluctuations, is dominant throughout the convective zone. It peaks near the surface where the negative (radially-inward) kinetic energy flux forces it to become locally greater than the stellar luminosity. At the bottom of the convective envelope, we note the presence of a negative enthalpy flux associated with overshooting. This requires some adjustment of the radiative diffusivity near the base of the convective envelope is necessary to ensure the full transport of the stellar luminosity, as is shown in the right panel of Figure 6. Such adjustment is done for all models (as in Brun et al., 2011). We will discuss in more details the properties of the overshooting layer in section § 4.2.

In the model shown in Figure 6, the luminosity $L(r)$ increases from the lower boundary until it reaches $L_{*}$. This is due to the presence of a volumetric heating source that mimics the production of internal energy by nuclear reactions. Only the stellar models with masses of 0.5 and 0.7 $M_{\odot}$ have such a heating source term as the models are deep enough that they include part of the nuclear core in their radiative zone.

4.1. Convective patterns

The development of the convective instability in the spherical shells leads to the classical network of downflow lanes surrounding broad upflows as can be seen in Figure 7. As the rotation rate increases (from top to bottom) we note that a higher azimutal wavenumber $m$ characterizes the convection. We clearly see that there are more downflow lanes near the equator, which are narrower, and smaller convective patches in the polar caps. At the fastest rates convection can be longitudinally modulated at low latitude, even exhibiting the so-called active nests of convection (Brown et al., 2008). Such active nests dominate the local transport processes and imprint their motions onto the top of the radiative zone. Except for the localized convective patterns, their overall behavior is consistent with models having a $R_{of}<1$, e.g. Reynolds stresses accelerate the equator (see Sections 5.1 and 7.1).

Notice that the low-latitude downflow lanes are bent into a banana-like shape, with its retrograde tails occurring at mid-latitudes. This is consistent with the statistical correlation of the radial and longitudinal motions that produce an efficient angular momentum transport through Reynolds stresses, as will be discussed in § 5 and § 7.1. The slowly rotating cases show a modest deflection in the opposite direction due to the prograde, high-latitude shear of the differential rotation that modifies the structure of the convective motions. At the low rotation rates of cases M05 S, M07 S, M09 S (not shown), and M11 R1, the convective patterns look more alike at all latitudes, much like a soccer-ball, and the mean flows are more randomly oriented. As demonstrated in Chandrasekhar (1961); Brun & Palacios (2009), at low rotation rate a dominant $\ell=1$ mode develops, as is clearly seen in the fluctuating temperature field. It wobbles around, making longitudinal averages hardly meaningful.

As anticipated in our introduction of a Rossby number scaling based on mixing length theory arguments, our numerical simulations show that for higher mass the convective rms velocity increases. For instance, going from about 10 $\mathrm{ms^{-1}}$ in a $0.5\,M_{\odot}$ star to hundreds of meters per second in a $1.1\,M_{\odot}$ star (see Table 3). The convective velocity amplitude realized in the 3-D models are higher than those evaluated using mixing length theory. This difference (of order 10-20%) in the model presented in this study is due to the inward kinetic energy flux found in the models that results in a slightly overluminous convective enthalpy flux and hence larger convective velocity (which we recall scales like $\sim\sqrt[3]{L_{*}}$, see section 2). This is a well documented physical mechanism in 3-D compressible turbulent convection whose amplitude depends on both turbulence degree and stratification, see for instance discussion in (Cattaneo et al., 1991; Brun & Palacios, 2009). We also note that as the rotation rate is increased the fluctuating rms velocity decreases due to a slightly decreased super-criticality of the system, whereas the full rms velocity increases due to the presence of a differential rotation. When considering only the fluctuating rms velocity, the flow field is found to be close to isotropy ($\tilde{v}^{\prime}_{r},\,\tilde{v}^{\prime}_{\theta}$ and $\tilde{v}^{\prime}_{\phi}$ of the same order). Note that for $\tilde{v}_{r}$ and $\tilde{v}_{\theta}$, the weak meridional circulation does not make much difference between the full and fluctuating components of the rms velocity. This is of course not true for $\tilde{v}_{\phi}$, which increases due to the differential rotation.

We also see that as the mass increases (from left to right in Figure 7), the convective flows are less constrained by rotation. This is clearly seen by comparing their respective Rossby number. In Table 4, we list three different flavors of this important number (see Appendix), and they each show that, for an identical rotation rate, the more massive stars possess the largest Rossby number, and thus the lowest level of rotational influence. Still, for a fixed Rossby number, we found that models of the more massive stars exhibit smaller convective patterns, which is likely due to their shallower convective zones.

The time dependence of the convective flow is very rich with continuous emergence, merging, cleaving of the convective cells and strong vortical downflows at the interstices of the downflow lanes for the low Rossby number cases. A clear advection to the right at low latitudes and to the left at high latitudes (e.g., respectively faster and slower than the global rotation rate) is evident for the cases rotating at the intermediate regime that possess a monotonic differential rotation. For the slowly rotating ones, convection is more isotropic and less influenced by mean flows. In the case with the highest rotation rate (Rossby number $\leq$ 0.1), the advection by the mean flow is less systematic due to the alternating prograde and retrograde jets (see § 5).

4.2. Overshooting in solar-like stars

At the base of the convective envelope of all the models presented in this work, a shallow region of mixing develops. These overshooting regions exist because the downward plumes do not immediatdiately halt their descending motions as they go through the swift transition in stratification, which is from unstable to convection (negative entropy gradient) to the stable radiative interior (positive entropy gradient). Indeed, these vortical downward-flowing structures require a finite distance before they overturn, since they have some incoming inertia as they cross the transition point. This distance depends upon the stratification, degree of turbulence, the rotation rate, and the thermal diffusivity (Zahn, 1991; Brummell et al., 2002; Browning et al., 2004; Rempel, 2004; Rogers et al., 2006; Chan et al., 2011; Brun et al., 2011; Masada et al., 2013). As they decelerate due to the action of buoyancy breaking, they induce extra mixing and turbulent flows across this finite layer. The superadiabatic to subadiabatic change in the stratification implies that the correlation between convective velocity and temperature fluctuations should reverse. This is illustrated in Figure 6, where the enthalpy flux is clearly negative at the base of the convective envelope. Since the amount of overshooting depends on the stiffness of the stratification of the radiative interior, this has motivated our choice of using a realistic stratification that is directly deduced from 1-D stellar evolution models. Theoretical studies have revealed that the Péclet number ($Pe=vl/\kappa$, with $l$ here a typical size of the plumes) of the individual plumes is the key quantity to assess the properties of the overshooting (see for instance Zahn, 1991; Brummell et al., 2002). Flows with a small $Pe$ number tend to overshoot, leading to a stratification that remains convectively stable. In contrast, flows with a high Péclet number require that the stratification become nearly adiabatic, being locally modified by efficient turbulent mixing, and the convection extends deeper into the radiative zone. Hence the use of convective penetration to describe this regime. It is clear that the degree of turbulence of our convective zone is mild, and hence our Péclet number small. In real solar-like stars, the Péclet number of turbulent plumes is much higher. We can thus expect that our simulations give an upper limit to the amount of penetration occurring in solar-like stars.

In the upper panels of Figure 8, we display the temporally and longitudinally averaged radial enthalpy flux profiles for models M07 with different rotation rates. In these meridional cuts we can see that the enthalpy flux is mostly concentrated in the convective envelope, and it is predominantly positive at all latitudes. Some inhomogeneities are apparent. They are likely due to the moderate degree of turbulence of the simulations, as discussed in Miesch et al. (2000). A negative enthalpy flux is observed at the base of the convection zone (delineated by the dashed black line). As the rotation rate increases (left to right), the radial enthalpy fluxes maxima drifts from the equator to the poles. In the lower panels, we show radial cuts of the radial enthalpy flux in the northern hemisphere. On these panels, we identify the location where the overshoot begins between the convective and radiative regions as the radial position ($r_{c}$), which is where the radial enthalpy flux crosses zero. The radial location where the overshooting ends is the radial position ($r_{0}$), which is where the radial enthalpy flux is only a $10\%$ of its local minima (Brun et al., 2011). The overshooting motions show a clear dependence upon latitude. This is quantified in Figure 9, where $r_{0}$ (solid lines) and $r_{c}$ (dashed line) are displayed between the equator to the latitude 60*∘* for the fast rotators in our sample (left panel) and the slow rotators (right panel). The colors label the masses of the models, as in previous figures. The overshoot region is wider near the poles in the slowly rotating models than it is in the solar-like cases. The fast rotating cases exhibit an interesting localized increase of the overshooting depth at smaller and smaller latitude when the mass of the star increases. The shape of the overshooting region is, as a result, sensitive to the rotation rate of the star, with slow rotators favoring a wider overshooting region near the poles and fast rotators at mid-to-low latitude.

We define the average overshooting depth $\bar{{\rm d}}_{\rm ov}$ by averaging the difference between $r_{c}$ and $r_{0}$ between the equator and latitude 55∘. We limit ourselves to these relatively modest latitudes to avoid any spurious averaging effects associated with the large temporal variations of the enthalpy flux at high latitudes. We display $\bar{{\rm d}}_{\rm ov}$ as a function of stellar mass and Froude number in Figure 10, where the rotation of the modelled star is indicated by the different symbols, and its mass by their color. The normalized overshooting depth generally decreases with mass for a given rotation rate, as is clearly seen on the left panel. The Froude number is defined as $F_{r}=\left(2\Omega_{\star}/N\right)^{2}$, where $N$ is the Brunt-Väisälä frequency. When the rotation rate is held constant, increasing the Froude number corresponds to a decrease of the Brünt-Vaisala frequency and thus eases the overshooting of convective plumes impacting the stable region. We indeed observe a clear trend (see left panel in Figure 10) with an overshooting depth increasing when the Froude number increases in all of our models. We finally note that some modulation of the overshooting depth may occur when magnetic fields are taken into account and will be fully characterized in a subsequent study (see Varela et al., 2016).

5. Various states of Differential Rotation

We now turn to a discussion regarding how differential rotation in rotating convection zone is established in various solar-like stars in order to interpret recent observational surveys (Saar, 2009; Reinhold & Reiners, 2013; do Nascimento et al., 2014; Garcia et al., 2014; Reinhold & Arlt, 2015).

5.1. Differential rotation profiles and amplitudes

As is evident in Figure 11, various states of differential rotation have been achieved in our parametric study. Differential rotation profiles and amplitudes are found to be influenced by both rotation rate and spectral type. First we see that for each given row (for respectively ’S’ models, one, three, and five times the solar rotation rate going from the top to the bottom of the figure), the angular velocity patterns change. Moreover, we observe the existence of three main states of differential rotation for which the location of prograde and retrograde longitudinal flows differ from one another. For instance, these three states are found near the solar rotation rate (second row), where there is an anti-solar-like profile (slow equator/fast poles) for the 1.1 $M_{\odot}$ model, a solar-like profile (fast equator/slow poles) for the 0.9 and 0.7 $M_{\odot}$ models, and a Jupiter-like profile (alternating zonal jets and a cylindrical angular velocity profile) for the 0.5 $M_{\odot}$ model at 3 times solar (third row). For the 0.5, 0.7 and 0.9 $M_{\odot}$ models, we also find anti-solar-like profiles when further reducing the rotation rate below the solar rate (models ’S’ in the upper row). In particular, the 0.9 $M_{\odot}$ model shows an anti-solar rotation profile if the rotation rate is half the solar rotation rate, 0.7 $M_{\odot}$ for $1/4$ solar rotation rate, and 0.5 $M_{\odot}$ for $1/8$ solar rotation rate, hence for rotation rates smaller than those deduced in § 2 using mixing length. Looking more closely at these simulations we notice the retrograde flow near the surface extends up to the tangent cylinder in models M07S and M09S. At the equator a zone of rapid rotation at the base of the convection zone is present in these two models. At high latitudes, the differential rotation exhibits fast flows akin to polar vortices already described in the literature (Brown et al., 2008; Featherstone & Miesch, 2015). Case M05S rotates so slowly that the longitudinal average is not well defined, hence the asymmetric profile observed in the upper left panel of Figure 11. This is due to a global dipolar mode of the convective flows which imprint itself on the overall dynamics.

As discussed in the introduction, the behavior of the differential rotation can be understood to be a result of the change in the amplitude of the Rossby number of the models. These three main states: anti-solar $(R_{of}>1)$, solar-like $(0.3<R_{of}<0.9)$, jupiter-like or cylindrical-banded $(R_{of}\lesssim 0.3)$ are encountered in our series of models. For those fastest rotating cases, the cylindrical differential rotation expected from Taylor-Proudman contraints transits from a monotonic behaviour from equator to pole into a banded structure of alternating jets. Those jets are commonly seen in planets like Jupiter and Saturn. Their spacing can be related to the compressible Rhines scale $\lambda_{g}$ as discussed in details in Gastine et al. (2014a). In models M05R3 and M05R5, $\lambda_{g}/2$ is found to be of the order of $15\%$ to $20\%$ of the stellar radius which is in qualitative agreement with the banded structure seen in Figure 11 for these two cases.

These results are also compatible with global 3D MHD simulations performed by other authors to model differential rotation and stellar magnetism in the convection zone (Miesch et al., 2006; Ghizaru et al., 2010; Racine et al., 2011; Käpylä et al., 2011; Augustson et al., 2015; Karak et al., 2015), particularly for solar like stars (Brun et al., 2004; Brown et al., 2010, 2011; Brun et al., 2011; Varela et al., 2016). These studies pointed out the large magnetic temporal variability and the critical effect of stellar rotation and mass on magnetic field generation through dynamo mechanism, which for some parameter regimes leads to cyclic activity (Gilman & Miller, 1981; Gilman, 1983; Nelson et al., 2013; Käpylä et al., 2013; Augustson et al., 2013; Guerrero et al., 2016; Augustson et al., 2015). The definition of these three main states allows a fast and straightforward identification of the expected magnetic temporal variability of the solar like stars, for instance the stars with anti-solar-like rotation profiles should exhibit smaller magnetic temporal variability than stars with solar like rotation profile, because the magnetic field regeneration via convective motions dominates the regeneration via differential rotation, leading to non oscillatory $\alpha^{2}$ dynamos instead of oscillatory $\alpha-\Omega$ dynamos (see, e.g., Varela et al., 2016).

5.2. Tachoclines

By considering the coupling between the convective envelope and the stable radiative interior for each modelled stars, we observe the natural development of a transition layer between the two zones for various physical quantities such as velocity, entropy and temperature fluctuations, and the overall dynamics. Among these variations, that of the rotation profile is crucial. In the Sun, this transition has been named tachocline (Spiegel & Zahn, 1992), and it is thought to play an important role in the organization of the eleven years cycle (Charbonneau, 2005; Brun et al., 2015a). How such tachoclines evolve in other solar-like stars with different global parameters is largely unknown, so we will assess that dependence here.

The overall shape of the tachoclines achieved in our simulations is shown in Figure 12, with radial cuts at several latitudes of the temporal and longitudinal average of the angular velocity for two representative cases. To help quantify these shapes, we fit the rotation profile in the tachocline with

[TABLE]

where parameter $A$ represents the amplitude of the differential rotation in the tachocline, $x_{c}$ its inward/outward drift, and $w$ its thickness. The fitted profile for models M09s and M09R3 are indicated in red in Figure 12.

Applying this methodology to all models, we obtain the following trends, summarized in Figure 13 for models M09. The amplitude $A$ of the differential rotation is found to increase with the rotation rate. This effect is predominantly seen at high latitude beyond the tangent cylinder. It may be due to the smaller lever arm in this region with respect to the equatorial one. The location $x_{c}$ of the tachocline shows a drastic change of behaviour between anti-solar (red circles) and solar-like cases. Indeed one notices it is closer to the surface at low latitude and deeper at mid-to-high latitudes whereas it is the opposite in the solar-like cases. The width $w$ of the tachocline, in part controlled by our choice of diffusivities, still exhibit a similar clear trend. The anti-solar cases possess a thicker tachocline at low latitudes compared to higher latitudes. Again, this trend reverses for the solar-like cases. The overall shape of the tachocline is summarized in the lower panels of Figure 13.

In order to disentangle the effect of viscous stresses from dynamical effects linked to rotation, we have further analyzed (not shown) a series of models with approximately the same rms Reynolds number ($\sim$ 100) spanning a fluid Rossby number from 0.28 to 0.84. By considering a constant Reynolds number for all cases we assume that the same degree of turbulence and convective advection with respect to viscouss effect are realized in the simulations and that we can even better isolate the role of rotation. The observed trends in shape, location and amplitude are confirmed with those models, confirming they originate from dynamical effects rather than being viscously controlled.

5.3. Scaling laws

We define the mean latitudinal contrast of differential rotation $\Delta\Omega$ as the difference taken at the top of the domain of the azimuthally and temporally averaged profile of $\Omega$ between the equator and latitude $60^{\circ}$. As a result, a positive $\Delta\Omega$ denotes a solar-like differential rotation with an equator rotating faster than the higher latitudes. The values of $\Delta\Omega$ are reported in Table 5 for all our models. We show them in the upper panel of Figure 14 as function of the fluid Rossby number $R_{of}$, coloured by mass and labeled by rotation rate. The latitudinal differential rotation generally drops with the fluid Rossby number and increases with the mass of star, and it possibly undergoes a transition in the anti-solar cases when $R_{of}>1$. The normalized differential rotation (middle panel) $\Delta\Omega/\Omega_{*}$ tends to show saturation at low Rossby number, and the anti-solar cases exhibit a transition similar to Featherstone & Miesch (2015). Note that due to a slightly different definition of the Rossby number, in Featherstone & Miesch (2015) the transition occurs around $R_{of}\sim 0.1$ while in our case it occurs around $R_{of}\sim 1$. Finally, this transition is less clear in the differential rotation kinetic energy (lower panel), which is defined by

[TABLE]

Indeed, the kinetic energy of the differential rotation decreases smoothly with the fluid Rossby number, more simulations with a larger Rossby number would be needed to confirm this lack of transition in the kinetic energy.

We use our set of models to fit its dependency upon mass and fluid Rossy number (as well as mass and rotation rate) and obtain

[TABLE]

If we do not retain the anti-solar cases in our regression fit, we have:

[TABLE]

The observational trends for the differential rotation rate exhibit either a similar dependency with rotation rate ($\Delta\Omega\propto\Omega_{\star}^{0.7}$, see Donahue et al., 1996; Saar, 2009) or a significantly lower one ($\Delta\Omega\propto\Omega_{\star}^{0.15}$, see Barnes et al., 2005; Reinhold et al., 2013), making a direct comparison difficult. One explaination could be the lack of dynamo generated magnetic fields in this series of hydrodynamic models. Such magnetic field should feedback on the global balance establishing the large scale differential rotation pattern (see § 7.1). Indeed, a preliminary study of the cases presented in this work that also include a dynamo field finds a lower exponent for the dependency upon $\Omega$ (see Varela et al., 2016), improving the comparison of our theoretical models with observational trends.

The trend with mass while retaining the same positive variation with $M_{*}$ also differs somewhat from the observational trends ($\Delta\Omega\propto M_{\star}^{5.6}$ in Barnes et al. 2005; Reinhold et al. 2013, $\Delta\Omega\propto M_{\star}^{5.4}$ in Collier Cameron 2007) . We attribute these differences in part to the non perfect relationship between mass and effective temperature in the observational data and also to a lower upper bound in stellar mass in our study compared to the observational data. Indeed in Augustson et al. (2012) we have simulated more massive F stars up to 1.4 $M_{\odot}$ and found a larger dependency of the differential rotation with stellar mass. Hence it is likely that the differential rotation contrast increases more steadily with mass between 1.1 and 1.4 $M_{\odot}$ given the shallow convective envelop in these stars than it does for lower masses (see also the mean field computations of Küker & Rüdiger, 2007; Küker et al., 2011, that confirms the trends found in 3-D simulations). More recent observations tends to also show a change of slope around the stellar spectral type F (see for instance Reinhold & Reiners, 2013).

6. Meridional Circulation

Meridional circulation results from the slight imbalance between the terms acting predominantly in geostrophic balance. These are the horizontal pressure gradients, Reynolds stresses, curvature terms and Coriolis force, in a purely hydrodynamical setup. As we change the aspect ratio and the rotation rates of our convective shells, we expect the relative amplitude of these terms to change as well, resulting in different meridional circulation profiles.

6.1. Profiles and Amplitudes

In solar-like star simulations, we expect the meridional circulation in convective envelope to be weak because geostrophy is mostly satisfied (Pedlosky, 1987; Ballot et al., 2007; Brown et al., 2008; Augustson et al., 2012). In Figure 15 we display the meridional circulations realized in the fifteen models as contours of the meridional streamfunction $\Psi$, defined as in Miesch et al. (2000):

[TABLE]

We note two main trends: solar-like models have multi-cellular flow structures, and the anti-solar ones possess mostly unicellular meridional circulations per meridional quadrant. As the Rossby number is decreased (from right to left), the number of cells increases, in particular near the polar cap. The amplitude of the meridional circulation in the convective envelope is of order of meters per second (see Table 3), being weaker for the low mass stars compared to the massive ones (as already discussed in § 4). In the radiative interior, this flow is extremely weak, the radial velocity dropping by several orders of magnitude. This results in a penetration of the meridional circulation of less than $3$ to $5\%$ of the stellar radius. The flow in the anti-solar cases is directed poleward in both hemispheres at the surface, with a return flow at the base of the convection zone. This is also true for the upper meridional cells in the faster rotating cases.

A more direct way to understand the maintenance of the meridional circulation is to consider the angular momentum balance. By splitting the right hand-side term between a global net torque coming from the difference between Reynolds and viscous stresses and from the advection of angular momentum by the meridional circulation (and by assuming stationarity), we get:

[TABLE]

where

[TABLE]

with ${\cal T}$ the global net torque (whose expression will be made more explicit in section 7.1), and $\varpi=r\sin\theta$ the moment arm.

If the global net torque ${\cal T}$ is zero, then there is no meridional circulation. In stars, we do not expect the viscous stresses to play a major role and hence the meridional circulation arises in order to compensate for the angular momentum transport due to convection (Reynolds stresses), provided that the magnetic effects are negligible. Since we are considering purely hydrodynamical cases and the viscous stresses necessarily contribute to the angular momentum transport in our simulations, the meridional circulation develops as a response to the net torque exerted by both the Reynolds stresses and the viscous diffusion of the differential rotation. The responding circulation to an applied torque is due to a physical mechanisms called gyroscopic pumping that generalizes Ekman pumping (McIntyre, 2007; Garaud & Bodenheimer, 2010; Brun et al., 2011; Miesch & Hindman, 2011; Featherstone & Miesch, 2015).

6.2. Scaling laws

We quantify the strength of the meridional circulation by calculating its associated kinetic energy defined as

[TABLE]

where $\left\langle\right\rangle$ stands for the azimuthal average and MCKE has been further averaged in time over 100 days towards the end of each simulations. We display in Figure 16 the trends of MCKE as a function of fluid Rossby number. We notice that the energy contained in the meridional circulation increases with Rossby number, and decreases with mass. The meridional circulation energy also increases with the Reynolds number (not shown here), which is naturally expected as the meridional circulation results from the imbalance between the turbulent Reynolds stresses and the way it advects angular momentum (see Featherstone & Miesch 2015 and §6.1). Finally, we use our set of models to fit the dependency of MCKE upon mass and fluid Rossby number and obtain

[TABLE]

7. Analyzing the Dynamics

In order to better understand the various dynamical states analyzed in the previous sections, we now turn to a more quantitative analysis of the dynamics. We will first look at the angular momentum redistribution in convective shells achieved in our simulations (§ 7.1), then on the thermal wind balance (§ 7.2), and finally of the energy exchange and maintenance of the differential rotation (§ 7.3).

As a preliminary analysis, we turn to Table 5 where various global quantities characterizing the dynamics achieved in our simulations are listed. We note that the kinetic energy contained in the convective shell increases with the rotation rate due to a global increase in the energy contained in the differential rotation (DRKE). By contrast, the energy contained in the convective motion (CKE=KE-DRKE-MCKE) is found to play a lesser role and decreases in overall amplitude and significantly in proportion. This is likely due to the propensity of rapidly rotating convective sphere to predominantly inject energy into the longitudinal flows. The kinetic energy contained in the meridional circulation is always small, and it globally decreases in amplitude to become almost negligible at the fastest rotation rates. In agreement with the DRKE, and as discussed in § 5, the absolute angular velocity contrast in latitude from 0 to 60∘ ($|\Delta\Omega|$) is found to increase with rotation rate, whereas the relative contrast decreases (see Figure 14). We will come back to the temperature and entropy contrasts in § 7.2.

7.1. Angular momentum balance

The differential rotation and meridional circulation profiles in our simulations are established and maintained through the transport of momentum and energy by convective motions that are influenced by the rotation, stratification, and spherical shell geometry.

Following Elliott et al. (2000) and Brun & Toomre (2002, hereafter BT02), an equation for the angular momentum transport can be deduced from the $\phi$-component of the momentum equation:

[TABLE]

with $L=r\sin\theta v_{\phi}$ the specific angular momentum and $\tau$ the net torque applied to the convective envelope. Assuming a statistically stationary state, applying longitudinal and temporal averages and writing the net torque as a divergence of a flux one gets (see also Pedlosky, 1987):

[TABLE]

where ${\cal F}_{r}(r,\theta)$ and ${\cal F}_{\theta}(r,\theta)$ represent the mean radial and latitudinal angular momentum fluxes whose expressions are given by:

[TABLE]

${\cal F}_{r,V}$ (resp. ${\cal F}_{\theta,V}$) is the flux associated to viscous transport, ${\cal F}_{r,R}$ (resp. ${\cal F}_{\theta,R}$) that related to Reynolds stresses and ${\cal F}_{r,MC}$ (resp. ${\cal F}_{\theta,MC}$) represents the angular momentum flux due to meridional circulation. Their detailed expression can be found in BT02 and in subsequent publications. As was done in BT02 we then integrate respectively each flux over colatitude and radius to assess the net flux through a sphere of varying radius and through cones of varying inclination:

[TABLE]

These integrated fluxes are presented in Figure 17 for our cases M07R1 and M11R1, and have been averaged over 5 rotation periods. For cylindrical cases such as M05R3,R5 the balance is very similar to that shown for M07R1. For simplicity we drop the letter $I$ when discussing the individual contribution of the flux. Since a statistically stationary state is realized in our simulations, the sum of all fluxes must be close to zero as there are no net torques left.

We start by discussing the angular momentum balance realized in case M07R1. We note that in the radial direction the prograde Reynolds stresses act to accelerate the equator opposed mainly by the viscous stresses and to a lesser extent by the meridional circulation. As discussed in Brun et al. (2011) and § 6, the meridional circulation is a response to the net torque applied by the sum of the Reynolds and viscous stresses, it helps reaching a stationary state. We see that it is both positive and negative, reflecting the presence of multiple cells. Convection is thus carrying angular momentum such as to accelerate the equator and slow down the deep layers. Turning to the latitudinal flux balance, here too the Reynolds stresses are found to carry angular momentum towards the equator (positive/negative in northern/southern hemisphere respectively), with a peak value near $15^{\circ}$ and a secondary peak near the latitude of the tangent cylinder. The viscous stresses are poleward as they tend to erase the differential rotation in the convective envelope. The role of the meridional circulation is mostly to transport angular momentum poleward at low latitude, with a small counter cell higher up. Overall the balance is well realized as shown by the sum being nearly equal to zero in both radial and latitudinal balances.

Turning now to case M11R1, we see that in the radial direction the Reynolds stresses are now transporting angular momentum inward, hence slowing down the surface and speeding up the deep layers. The viscous stresses and meridional circulation have very similar amplitudes and profiles and transport angular momentum outward. In this case, the flux associated with meridional circulation is positive at all depths and relatively large, reflecting the mostly unicellular profile of the meridional flow (illustrated in Figure 15). In the radial direction, the angular momentum balances in cases M07R1 and M11R1 are thus of an opposite sign, in agreement with their differential rotation being respectively solar-like and anti-solar-like. In the latitudinal direction, the situation is not as dissimilar as one could have anticipated. Indeed, the Reynolds stresses have not changed sign with respect to case M07R1. What has changed is the amplitude and profile of the meridional circulation, with the viscous stresses now aiding the Reynolds stresses instead of opposing them. This is again due to the sign of the differential rotation. With a fast pole and a slow equator, the viscous stresses tend to slow down the pole and speed up the equator. Through gyroscopic pumping, this results from a meridional circulation profile that transports angular momentum poleward such that a stationary state is established. These subtle differences between solar-like and anti-solar-like cases are in agreement with the results from Featherstone & Miesch (2015), which were obtained without an underlying stable zone. As can be seen by the solid curve, the latitudinal angular momentum fluxes are not strictly balancing one another for case M11R1. This is in contrast to the three other panels, and it is due to the slow equilibration of this case.

7.2. Thermal wind balance

We now turn to analyze the role played by thermal effects and heat redistribution in the dynamical balance realized in our 3-D stellar models. As published in Brun et al. (2010), a general meridional force balance equation can be derived that reveals the subtle role of all processes in maintaining a non-cylindrical rotation profile that differs from the “classical” thermal wind (TW) balance (Durney, 1999; Brun & Toomre, 2002; Balbus et al., 2009). It is straightforward to use our numerical simulation to evaluate what are the dominant terms and how this meridional force balance comes about. Its derivation is summarized in Appendix  C and can be compactly written as

[TABLE]

where the various right-hand side terms are:

• $\bullet$

$\mathcal{S}$ describes the stretching and tilting of the vorticity due to velocity gradients;

• $\bullet$

$\mathcal{A}$ describes the advection of vorticity by the flow;

• $\bullet$

$\mathcal{C}$ describes the change of vorticity due to compression;

• $\bullet$

$\mathcal{B}$ is the so-called baroclinic contribution from the latitudinal entropy gradient and the product of the radial background entropy gradient with the latitudinal pressure gradient. The former dominates when the stratification is nearly adiabatic;

• $\bullet$

$\mathcal{V}$ accounts for the viscous diffusion of vorticity.

Under the assumption that the convection zone is nearly adiabatic and hydrostatic, that the fluid Rossby number $R_{of}$ is small, and that viscous stresses can be neglected, Equation (25) simplifies to:

[TABLE]

This is the “classical” thermal wind equation. It states that baroclinicity can break the Taylor-Proudman constraint of $\partial\langle v_{\phi}\rangle/{\partial z}=0$, implying a cylindrical rotation profile (Zahn, 1992). This is due to the fact that baroclinic torques suppress the Coriolis-induced meridional circulation that would otherwise tend to establish a cylindrical state of rotation. It is instructive to use our numerical simulations to evaluate the role played by all the terms of the zonal vorticity equation identified above and to discuss the nature of the meridional force balance.

In Figure 18, we display the full thermal wind balance achieved in two models having respectively a prograde and retrograde differential rotation profile, e.g. M07R1 and M11R1. For each case the lhs and the various terms composing the rhs are shown, using the same color table and minimum and maximum bounds. The first point to notice is that for both models rhs $=$ lhs to a high degree of fidelity, meaning that they have achieved an equilibrium and a well-relaxed state. Turning to M07R1 (top row of Figure 18), we see that the dominant term is the baroclinic one. It is mostly negative (positive) in the northern (southern) hemisphere as is the lhs and possesses elongated island features in the radiative zone. In that later zone, the “classical” thermal wind balance is realized. In the convective envelope, in particular near the surface, this is less the case. Stretching ($\mathcal{S}$) and advection ($\mathcal{A}$) terms, and to a lesser extent the viscous term ($\mathcal{V}$), contribute to the overall balance. The first two of those terms have the same sign in each hemisphere as the baroclinicity ($\mathcal{B}$), whereas $\mathcal{V}$ has the opposite sign. This confirms the role played by convective motions in the meridional balance, and its departure from a strict classical thermal wind balance. This can be easily understood by the fact that $R_{of}$ is not very small in that model (see Table 4).

This departure from a strict TW balance is even more apparent in model M11R1 shown in the bottom row of Figure 18. With $R_{of}>1$, the full TW balance must be considered as most terms are expected to contribute to the rhs. The baroclinic term shows the largest and most systematic contribution to the overall balance. In that case too, it dominates the balance in the radiative interior and possesses a predominantly positive contribution from the low to mid-latitudes up to the poles. Case M07R1 has the opposite response, showing one important difference between prograde and retrograde cases. This is directly linked to the different entropy fluctuation profiles realized in the models, as we will discuss just after in commenting Figure 19. Nevertheless, $\mathcal{S}$, $\mathcal{A}$, $\mathcal{C}$ and $\mathcal{V}$ terms now contribute everywhere in the convective envelope and not only near the surface of the simulation domain. Note that in that slowly rotating case, $\mathcal{S}$ and $\mathcal{A}$ have an opposite sign with respect to $\mathcal{B}$ and contribute a little in the radiative zone due to the development of turbulence in the overshooting region by the more vigorous convective downdrafts. In both models (this is also true for the other cases), $\mathcal{S}$ and $\mathcal{A}$ have in the large the same sign. This is consistent with the analysis of the latitudinal angular momentum transport performed in the previous section, where their sign was found to remain the same when transiting from $R_{of}<1$ to $R_{of}>1$. This is an important property of prograde and retrograde models, and it holds for all the cases considered here. Such an invariance indicates that the turbulent latitudinal Reynolds stresses do not change sign as the rotation rate is varied, whereas the entropy and temperature fluctuations do.

We conclude from this TW balance analysis that for low $R_{of}$ the baroclinic term remains the key player in tilting the iso-contours of $\Omega$ with some contribution from turbulence. However this becomes inefficient if $R_{of}$ becomes too small $<0.1$ as is the case for models M05R3 and M05R5 that are mostly cylindrical. This is in part due to the increased influence of rotation that cannot be fully compensated by thermal effects as we will see when commenting upon the trends shown in Figure 19. For high $R_{of}$ cases the contribution of all terms composing the rhs is more balanced, but their cumulative action also tends to slightly bend the iso-contours of $\Omega$. However, in such cases, the Taylor-Proudman constraint does not play as important a role as in the very low Rossby number cases.

In order to understand further how the baroclinic term arises, it is useful to look at meridional cuts of the entropy and temperature fluctuations averaged over time and longitudes as shown in Figure 19 (top row). We note first that the fluctuations of $S$ and $T$ are symmetric with respect to the equator for both cases (this is also true for all models except M05s) and that their latitudinal variations are of opposite sense. In M07R1, the poles are warm and the equatorial regions cool, whereas it is the reverse for M11R1. Such a latitudinal gradient of entropy (and temperature) is consistent with the baroclinic term $\mathcal{B}$ shown in Figure 18 and discussed earlier. The amplitude of the variations is also quite different, with case M11R1 having fluctuations about a factor of 10 larger. In M07R1, $\Delta T\sim 10K$ whereas it is around $100K$ in M11R1. Given the value of the background temperature at the base of the convective envelope of each cases ($\sim 3-6MK$), these variations can be seen as tiny but there are large enough to drive large baroclinic torques due to the large heat capacity of stellar plasmas. The larger variations seen in the M11R1 case is in agreement with the more vigorous convective flows realized in the simulation due to the significantly larger stellar luminosity that needs to be carried out in this star with respect to a 0.7 $M_{\odot}$ star (see MLT discussion in §2). In the radiative interior we further see that the variations are the largest where the rotation profile transits from differential to uniform, e.g. in the tachoclines that possess the largest radial shear as discussed in Brun et al. (2011). We can understand how such states are established in our models by looking at the latitudinal heat balance as we will discuss soon, after commenting on Figure 20.

However, we first discuss how the latitudinal entropy and temperature fluctuations contrast vary with mass and rotation rate by considering the 15 models of this study. In the bottom row of Figure 19, we show how $\Delta T$ and $\Delta S$ vary with $R_{of}$ for the four masses considered here. We chose to compute $\Delta T$ and $\Delta S$ at the surface and between 0 and 60∘ of latitude and to plot their absolute value to avoid artificial behavior due to their change of sign as $R_{of}$ is varied. We note that for faster rotation rates both entropy and temperature contrasts increase, reaching amplitudes of several hundred $K$ for the temperature. Likewise, we find that the fluctuations grow with increasing stellar mass, which is in concordance with the more increasingly intense convection achieved in those more massive stars. As seen in (Brown et al., 2008; Augustson et al., 2012), the latitudinal entropy and temperature contrasts have a strong mass dependence and a weaker one on rotation (or Rossby number). These trends can be summarized by deriving scaling relationships for $\Delta T$ and $\Delta S$, which are obtained through multi-parameter regression fits to the data shown in Figure 19. We show the resulting scalings and their uncertainties for two cases: one including the anti-solar cases and the other not. We see that those cases influence the exponents of the scalings, but not the overall trends.

With the anti-solar cases:

[TABLE]

Without the anti-solar models:

[TABLE]

The latitudinal variation of the thermodynamics variables $S$ and $T$ is established by anisotropic heat transport in the convective envelope. This is an expected behavior in a rotating convective envelope, as the Coriolis force varies with respect to latitude and so does its influence on the convective motions (see for instance Durney, 1989, 1999; Elliott et al., 2000; Miesch, 2005; Brun & Rempel, 2009; Käpylä et al., 2011). Hence, along with the radial heat transport realized in the convective shell of our models, as discussed in §4, there is a latitudinal heat transport that is worth studying in greater detail. As done in Elliott et al. (2000) and Brun & Palacios (2009), we can decompose the latitudinal heat transport into various terms involving diffusion and advection processes. In Figure 20, we display, using the same color table and minimum and maximum bounds, the five terms contributing to the latitudinal transport of heat after having converted and normalized them to the adequate stellar luminosity. We note that the main terms are the latitudinal enthalpy and entropy fluxes, with a minor contribution from the kinetic energy flux. In this context, the viscous and radiative fluxes are negligible. To see their structure, they have been multiplied by very large factors (e.g., $10^{14}$ and $10^{5}$ respectively) in Figure 20. Hence, the balance is mainly between the entropy flux and the enthalpy flux, which are mostly of opposite sign. The leading contribution comes from the correlation of the fluctuating latitudinal velocity and the temperature fluctuations to yield an enthalpy flux that efficiently transports the heat, hence establishing the fluctuations seen in Figure 19. Once these profiles have been established, the latitudinal entropy flux develops to establish a balance. Turning to case M07R1, we see that the entropy flux tends to cool down the polar regions and heat the equatorial region, and it is maximum at the base of the convection zone. For case M11R1 (lower panels), the balance is such that the entropy flux mostly warms up the poles. In a small equatorial region, the entropy flux changes sign and tends to warm up this zone in agreement with Figure 19. In that case, the role of the enthalpy flux is less clear, but it tends to also balance the entropy flux even though the system has not yet reached a complete equilibrium. As a summary, we see that the heat transfer and the associated thermodynamic perturbations in the solar and anti-solar cases differ, yielding hot poles and cool equator for the solar-like case, and the reverse for the anti-solar ones.

7.3. Energy exchange and maintenance of differential rotation

It is also useful to understand the energy exchanges maintaining the differential rotation. Following Rempel (2005), the kinetic energy balance for the differential rotation can be written

[TABLE]

where the various terms represent the work of the Coriolis force, Reynolds stresses, ${\rm DRKE}$ advection and viscous forces. Their expression is given by (Rempel, 2005) and defined in Appendix  D.

The Coriolis force relates to energy transfers between the differential rotation and meridional circulation. The non-linear advection is decomposed into Reynolds stresses associated with the radial ($Q_{r}^{r}$) and latitudinal ($Q_{r}^{l}$) profiles of the differential rotation, and the differential rotation advection $Q_{A}$. Finally, we regroup all the contributions from the viscous stress tensor into $Q_{V}$ and the geometric terms into $Q_{\rm curv}$.

We display in Figure 21 the various contributions to the kinetic energy balance of the differential rotation for our 15 models (the numerical values can be found in Table 6). In all models, the differential rotation is mainly maintained by a balance between the non-linear Reynolds stresses and the viscous stress. The Reynolds stress contribution itself is strongly dominated by the radial component. It only moderately increases with rotation except for models with $M_{\star}=0.9\,M_{\odot}$. In less massive stars, the convective flows maintaining the differential rotation correspond to a higher fraction of the stellar luminosity than in more massive stars. This behaviour can be naively expected from the simple scaling laws derived in Section 2, namely $\rho_{\rm bcz}v^{2}\sim M_{\star}^{-0.9}$. Finally, the advection of the differential rotation $Q_{A}$ is completely negligible in all cases.

It is also instructive to note that for the 4 models possessing an anti-solar differential rotation (M05s, M07s, M09s and M11R1), the latitudinal transport of energy by Reynolds stresses is negative and the differential rotation is mainly sustained by the radial transport of angular momentum by Reynolds stresses. In these models the Rossby number is higher than 1 (see Table 4), the convective motion are much less vortical than in the other cases and this directly impacts how the Reynolds stresses transfer kinetic energy to the differential rotation.

Following Rempel (2005, 2006), we can finally identify the amount of energy transfered to the large-scale differential rotation with $|Q_{V}|$. We note in Table 6 that $|Q_{V}|$ increases with rotation rate and decreases with mass. For instance the model with $M_{\star}=0.5\,R_{\star}$ and $\Omega_{\star}=5\Omega_{\odot}$ converts about $23\%$ of the stellar energy flux into the large-scale differential rotation. Note that this model is twenty times less luminous than the Sun, which means that only $1\%$ of a solar luminosity is converted into differential rotation (second column in Table 6). Model M11 R5 converts the most energy (in absolute value) into differential rotation with $|Q_{V}|\sim 21\%L_{\odot}$. We recall that $|Q_{V}|$ is only a proxy of the converted energy, as the differential rotation is maintained due to the angular momentum transfer by the turbulent Reynolds stress (see § 7.1).

8. Discussion and Conclusions

In this study we have focused our analysis on the characterization of the angular velocity profiles in the convective envelope of solar-like stars of various masses and rotation rates. Starting from a mixing length argument, we have shown that a Rossby number can be derived based on fundamental stellar parameters such as mass and rotation rate that gives an important insight into the expected differential rotation profile. Three main categories have been identified: anti-solar-like, solar-like and cyclindrical/Jupiter-like. These categories depend strongly on the Rossby number of the simulations, being large for anti-solar cases and smaller than $R_{of}<0.1$ for Jupiter-like profiles.

Solar-like rotation profiles are found for intermediate $R_{of}$ values that are less than unity. Such differential rotation states with conically-tilted isocontours of $\Omega$ are achieved only if the baroclinic term plays a dominant role in the thermal wind balance, as discussed in §7. As the rotation rate is increased and the Rossby number decreases toward 0.1, the angular velocity profile tends to become more cylindrical, retaining its monotonic behavior with fast equator and slow poles. Then for even lower Rossby numbers, the alternating prograde and retrograde jets become increasingly apparent. Eventually, they tend toward Jupiter and Saturn’s surface angular velocity profiles. A good guess of the number of jets in such cases can be obtained by computing the compressible Rhines scale (see e.g. Gastine et al., 2014a).

Of course the values quoted for the Rossby number are somewhat dependent on the definition used. Observers tend to favor the stellar Rossby number, comparing the rotation period to the convective overturning time deduced from stellar evolution models (Landin et al., 2010). Even in that case, there is some freedom in choosing the convective over turning time, where one can consider either the base of the convective envelope, the mid-depth value, or an averaged value. Another definition of the Rossby number can be used, the so-called convective Rossby number as discussed in Glatzmaier & Gilman (1982). This is relevant in analyzing rotating convection simulations such as the ones presented in this study, because it uses key nondimensional numbers such as the Rayleigh, Taylor and Prandtl numbers. It evaluates the ratio between the driving buoyancy force and the rotation constraints, getting rid of diffusion effects. A third definition that is not employed in this study, but that is found to be relevant for the study of stellar dynamos, is the scale dependent Rossby number. It is evaluated by comparing the kinetic energy contained in the first convective modes to the energy contained in all the scales. Finally, there is the fluid Rossby number, which is used in most of this work. It relates the turbulent vorticity to the planetary vorticity. It is an a postiori measurement of the effect of the Coriolis force on the turbulence. We find it useful because it is quite straightforward to compute, and it gives a robust assessment of the rotating dynamics. The transition values given above and throughout the paper are based on the fluid Rossby number.

However, for the sake of completeness, we have also quoted in table 4 the convective and stellar Rossby numbers. What is clear from this table is that those numbers are indeed different but their relative trends with respect to global stellar parameters follow a simple quasi-linear relationship. Hence, what matters more than the magnitude of the Rossby number at the point of the the transition is to distinguish the rotational regimes. For large Rossby numbers, we expect to have anti-solar behavior, for intermediate values to be solar-like, and for small values to have a cylindrical rotation profile that can be either monotonic with respect to latitude or exhibit alternating zonal jets when the rotational constraint is very large.

Intrinsically, mixing length theory cannot predict differential rotation states. It can only illustrate how the Rossby number changes with respect to global stellar parameters as demonstrated in §2 and Figure 1. In contrast, the natural outcome of our 3-D numerical simulations of rotating stellar convection is to provide, among other dynamical properties, the state of rotation achieved in a model for a given set of stellar parameters. We can thus assess, with a consistent definition of the Rossby number, what are the resulting states of differential rotation for the simulations. In doing so, we can calibrate the nondimensional constant $c_{1}$ used in plotting Figure 1. Note that no attempts at fine tuning have been made to precisely find the transitions, instead a systematic scan of the stellar parameters has been performed so as to obtain the three differential rotation states identified. With a multi-parameter regression fit to our set of 15 numerical simulations, we obtain the following scaling relation between the fluid Rossby, stellar rotation and mass:

[TABLE]

We note that it is not far from the one we derived in §2 from back of the envelope arguments using MLT. The constant $c_{1}\sim 0.89$ and the exponents are very close as well. To further illustrate this important result, we display this scaling relation based upon our set of 3-D simulations of rotating stellar convection in Figure 22 using color contours. What our study confirms and Figure 22 summarizes is that for a given observation of the stellar rotation rate or $v\sin i$, the rotation state depends on the stellar mass of the solar-like star. In particular, for a given rotation rate, the more massive stars may be more likely to achieve an anti-solar-like rotational state. This trend could be observable and a systematic search for anti-solar stars should be undertaken, as it will greatly aid in the constraints upon our models. And indeed some attempts using Kepler data have already been started (Reinhold & Arlt, 2015; Varela et al., 2016).

It is interesting to note that the transition of the state of differential rotation with Rossby number that we have indentified could have a direct correspondance with trends found in stellar X-ray luminosity studies (see, e.g. Pizzolato et al., 2003; Wright et al., 2011, and references therein). Jupiter-like profile could correspond to the saturated X-ray regime, solar-like profiles to the linear regime and the anti-solar state to the well less defined regime found in these studies. Also this modification of the large differential rotation profile could lead to different dynamo regimes and magnetic field topology, impacting directly wind braking and stellar spin-down (Brown, 2014; Vidotto et al., 2016; Metcalfe et al., 2016). We intend to verify these relations and changes of regime with dedicated dynamo simulations.

We also find that $\Delta\Omega$ varies significantly in amplitude with stellar mass, its amplitude being larger for more luminous stars as observed by Barnes et al. (2005). It is also possible that the dependency (exponent) vs stellar mass becomes steeper for F-type stars than for G and K stars (see, e.g. Augustson et al., 2012). The trend with rotation rate is also to have a larger contrast for faster rotation rate, but the relative differential rotation $\Delta\Omega/\Omega_{*}$ is found to decrease with $\Omega_{*}$. The scaling relationships have a larger dependency than is advocated by Collier Cameron (2007), but they are in close agreement with the studies of (Donahue et al., 1996; Saar, 2009). Furthermore, we know that in MHD dynamo models the angular velocity contrast has a weaker dependence on the rotation rate, as demonstrated in Varela et al. (2016) which utilizes MHD versions of the 3-D simulations discussed here. So, the scaling relationships given in these paper with respect to the rotation rate should be considered as upper limits.

All our simulations include a stable radiative interior below the convective envelope. The seamless nonlinear coupling of these regions greatly improves the realism of the bottom boundary condition of the convective layer relative to an impenetrable one. The convective downdrafts can plumet through the domain without hitting a solid wall, and so they can be buoyantly braked. A careful analysis of the overshooting layer that results from the pummeling of the convective motions reveals interesting trends. We find that the amount of overshooting decreases with mass, when it is characterized by latitudinally averaging the radial extent over which enthalpy flux is negative at the base of the convective envelope. We further find that stars with prograde rotation have a smaller normalized overshooting extent with increasing Rossby number, whereas it is the reverse for anti-solar-like stars. The latitudinal variations are interesting too. We find that solar-like stars have prolate overshooting layers and that in anti-solar-like stars they are oblate.

As for the overshooting layer, characterizing the shape and amplitude of stellar tachoclines at the base of convective envelopes is of importance for stellar magnetism and chemical mixing. We find the following trends in our study: the prolateness of the tachocline changes with rotation rate. Anti-solar-like star have oblate tachoclines, thicker at the equator, whereas solar-like star have prolate tachoclines. We also find that for faster rotation rates the overall shear is larger, in agreement with $\Delta\Omega$ being larger in the convective envelope. However, our choices for the radial profile of the thermal and viscous diffusivities influence the thickness and location of the tachoclines. Indeed, since the pioneering work of Spiegel & Zahn (1992), it has been well known that the viscous and radiative spreading of a tachocline in the radiative interior of a solar-like star is at work, unless there is some yet-to-be-identified physical process that acts against it. Several scenarios have been proposed: anisotropic turbulence, gravity waves, primordial or cyclic dynamo-generated magnetic field. However, a specific numerical setup will likely be required to be able to disentangle their various impacts, which is a task far beyond the scope of this study. Nevertheless, we believe that the trends found for the shape of the tachoclines is robust as all models have been built with the same physical ingredients and numerical accuracy.

Meridional circulation is found to change significantly with rotational influence. Its shape changes from a monolithic unicellular poleward flow in the anti-solar cases, to multi-cellular flows, with the number of cells increasing both in radius and latitude, for solar-like and Jupiter like states of differential rotation. This can be understood by the ability or not of the meridional cells to extend beyond the latitude corresponding to the tangent cylinder of each models. With stronger rotational constraint, the meridional cells are more and more aligned with the rotation axis and are confined to lower latitudes as discussed in Featherstone & Miesch (2015). Our study also confirms previous findings that the meridional circulation weakens as the rotation rate is increased. This is linked to the fact that more kinetic energy is being channeled to the longitudinal motions. However, it could be the case that this decrease in the amplitude of the meridional circulation with respect to the rotation rate is due to a lower level of supercriticality of the simulations, since it is difficult to maintain it. Indeed, we know that this large-scales flow is a direct response to any net longitudinal torque applied to the envelope through the effect of gyroscopic pumping. As a consequence, if the amplitude of the Reynolds stresses weakens because of a less intense degree of turbulence, then the meridional circulation will be weaker as well. Even taking into account the decrease of amplitude of the Reynolds stresses, the global trend of weaker meridional circulation for faster rotation rate is confirmed. As we clearly see in Table 4, the Reynolds number does not decrease. On the contrary, it can increase.

Two cases exhibit active nests of convection (M05R3, M05R5), but they do not impact the conclusions derived in this study. Their co-existence with an underlying stably-stratified region in our models warrants further investigation regarding their formation, which we intend to explore in the near future with a specific parameter study.

To summarize, we have seen that many properties of stellar convection, with its associated mean flows and transport mechanisms, are influenced by rotation and stellar mass. We have been able to find useful trends and to anticipate the rotation regime of candidate stars thanks to MLT and 3-D numerical simulations. This study was done with purely hydrodynamical models and lacks a consideration for the nonlinear feed backs related to a dynamo generated magnetic field. Preliminary studies by our group, and published in Varela et al. (2016), seem to confirm the main trends for the differential rotation states but with a weaker sensitivity to global parameters. Work by Käpylä et al. (2014); Gastine et al. (2014b); Karak et al. (2015); Guerrero et al. (2016) also find some differences arising from magnetic fields, but they do not change the global trends presented in this study. For instance we find that the anti-solar cases presented in this hydrodynamical study remain anti-solar when magnetic field is taken into account (see Figure 2 of Varela et al., 2016, for a preliminary study of our MHD antisolar cases). This is likely due to the fact that their Rossby ($R_{os}$ or $R_{of}$) remain larger than 1.0 even after the Lorentz force has influenced the convective flows. We expect to publish a detailed analysis of the dynamo counterpart of the 15 hydrodynamical models discussed in this study in the near future.

Of course the turbulence degree used in the simulations discussed in this work are still limited by the current computer resources and one must be extremely careful in comparing directly numerical results with observations. Still we find these systematic parameters study useful to delineate the main trends and identify the key physical mechanisms and it is reconforting to see that the observational tendencies are recovered qualitatively. We are also convinced that anti-solar like differential rotation states are worth searching for observationaly by selecting stars with large fluid or stellar Rossby numbers.

A.S. Brun and A. Strugarek dedicate this paper to Professor Jean-Paul Zahn, whose continuous support and advice over the years and great expertise in stellar fluid dynamics have been of invaluable help in analyzing non linear 3-D simulations and in improving our understanding of stellar structure, evolution, and dynamics in general. We miss him dearly. We also thank Paul Charbonneau, Steve Saar, Thomas Gastine, Michio Yamada, Shin-Ichi Takehiro and Rafael Garcia for useful discussions. We acknowledge funding by ERC STARS2 207430 grant, ANR Blanc Toupies SIMI5-6 020 01, INSU/PNST, CNES SolarOrbiter, PLATO and GOLF grants, FP7 SpaceInn 312844 grant, and NASA grants NNX11AJ36G, NNX13AG18G and NNX16AC92G. K. C. Augustson is funded through the ERC SPIRE 647383 grant. A. Strugarek acknowledges support from the Canadian Institute of Theoretical Astrophysics (National Fellow), from Canada s Natural Sciences and Engineering Research Council and from CNES postdoctoral fellowship. Simulations have been performed on GENCI and PRACE supercomputer infrastructures under grant 1623 and RA1964. A.S. Brun wishes to thank the University of Colorado and JILA as well as the University of Kyoto and RIMS for their hospitality.

Appendix A Model ingredient parameters

In Table 7 we list some of the parameters used in the simulations.

Appendix B Rossby numbers

There are multiple definitions of the Rossby number in the literature that quantify the influence of rotation on the dynamics of the system. We chose to use the fluid Rossby number $R_{of}$, which is a direct comparison of the advection term and the Coriolis force in the Navier-Stokes equation. It is defined by

[TABLE]

where $\tilde{\omega}$ is the rms vorticity at mid-depth in the convection zone. We choose to evaluate the Rossby number at mid-depth of the convection zone as it is close to the location of the maximum angular momentum transfers shown in Figure 17. We also tested a different definition of the Rossby number for which $\tilde{v}$ is averaged over the all convection zone rather evaluated at mid-depth and did not find any significant difference in the analysis reported in this work.

The stellar Rossby number, which is often used in the literature, corresponds to the ratio of the rotation period of the star ($P_{\rm rot}=2\pi/\Omega_{*}$), to the convective over-turning time of convection ($\tau_{\rm conv}=d_{\rm CZ}/\tilde{v}_{r}$, with $d_{\rm CZ}$ the thickness of the convective envelope). With the definition of the convective over-turning time being different for various authors, it is generally deduced from stellar structure models using mixing-length at various locations in the convective envelope, leading to some confusion in its exact definition (Landin et al., 2010). In our simulations, we directly use the mid-depth value of the radial velocity. It is defined by

[TABLE]

The convective Rossby number $R_{\rm oc}$, first introduced by Gilman & Glatzmaier (1981), is a combination of the Taylor, Rayleigh and Prandtl numbers and is defined by

[TABLE]

Finally, a modified Rossby number $R_{\rm ol}$ was introduced by Christensen & Aubert (2006) to take into account the characteristic length scale of the flow rather than the shell thickness $d_{\rm CZ}$. It is defined by

[TABLE]

where $U$ is the rms velocity at mid-depth, $L$ the size of the convective enveloppe, and the characteristic length scale $\bar{l}_{u}$ is defined as

[TABLE]

where ${\bf v}_{l}$ is the velocity field at scale $l$ in the spherical harmonics spectral space.

Appendix C Thermal wind balance

The complete thermal wind balance equation can be derived from the vorticity equation

[TABLE]

with \mbox{\boldmath\omega}_{a}=\mbox{\boldmath\nabla\times}{\bf v}+2\mbox{\boldmath\Omega_{*}} the absolute vorticity and \mbox{\boldmath\omega}=\mbox{\boldmath\nabla\times}{\bf v} the vorticity in the rotating frame. Averaging the zonal component of this vorticity equation over longitude and time and assuming a statistically stationary state yields the general equation for force balance in the meridional plane:

[TABLE]

where $\displaystyle\frac{\partial}{\partial z}=\cos\theta\frac{\partial}{\partial r}-\frac{\sin\theta}{r}\frac{\partial}{\partial\theta}$ and

[TABLE]

Appendix D Kinetic energy balance for the differential rotation

We recall the evolution equation of the kinetic energy associated with the differential rotation:

[TABLE]

The various terms of Equation D1 are defined by

[TABLE]

where $\left\langle.\right\rangle_{r,\theta}$ stands for the average over the meridional plane, and the overbar stands for the azimutal average.