NUMERICAL SIMULATIONS OF A ROTATING RED GIANT STAR. I. THREE-DIMENSIONAL MODELS OF TURBULENT CONVECTION AND ASSOCIATED MEAN FLOWS

Almost all, if not all stars rotate, and the distribution of their angular momentum appears to change along the course of their evolution as can be seen from the υ sin i data collected during the last decades for stars all over the HR diagram (e.g., see for instance Pilachowski & Milkey 1984, 1987; de Medeiros & Mayor 1999; Glebocki & Stawikowski 2000; Royer et al. 2002a, 2002b; Barnes et al. 2005; Karl et al. 2005; Carney et al. 2008).

To constrain the internal angular momentum profile at each evolutionary phase, the easiest approach would be to probe stellar interiors in order to derive the angular velocity profile. Helioseismology allows us to invert the inner solar rotation profile down to r = 0.25R_☉ (Schou et al. 1998; Antia & Basu 2000; Antia et al. 2008; García et al. 2008). Asteroseismology should soon offer similar opportunities for other stellar spectral types thanks to the CoROT (Goupil et al. 2006; Baglin et al. 2007) and Kepler satellites. In the meantime, however, we are left with indirect probes of the internal angular momentum distribution in all stars: (1) the surface velocity measurements for stars of similar initial mass at different evolutionary stages; (2) the surface abundance anomalies resulting from the action of internal transport processes, in part due to rotation, that connect the stellar envelopes to the nuclearly processed internal regions. These indirect probes can be compared to the results of one-dimensional rotating stellar evolution models (Pinsonneault et al. 1989; Fliegner & Langer 1994; Meynet & Maeder 1997; Talon et al. 1997; Maeder & Meynet 2000; Heger et al. 2000, 2005; Palacios et al. 2003, 2006; Chanamé et al. 2005; Suijs et al. 2008) in order to guess the angular velocity profile evolution. The introduction of rotation and associated transport processes in one-dimensional stellar evolution models leads indeed to a clear improvement of the comparison between theoretical predictions and observations, in particular from the point of view of chemical abundances. However, a number of points, in particular concerning the evolution of the surface rotation velocities, remain to be elucidated. It is the case for the differences between predicted and observed spin rates for white dwarfs and neutron stars (Heger et al. 2005; Suijs et al. 2008), or the interpretation of the surface velocities distribution of horizontal branch stars (Sills & Pinsonneault 2000; Recio-Blanco et al. 2002, 2004).

While important theoretical developments have been devoted to the description of the one-dimensional angular momentum transport in the radiative stellar interiors (Zahn 1992; Maeder & Zahn 1998; Mathis & Zahn 2004), the transport of angular momentum in convective regions, and in particular in extended convective envelopes, is poorly described in one-dimensional models. This is of particular importance for the case of giant stars that possess very extended convective envelopes, occupying up to 80% and more of the total stellar radius. These stars are expected to have very slow surface rotation due to their large radius, which is confirmed by observations (de Medeiros & Mayor 1999; Carney et al. 2008), but may also have a rapidly rotating core (Sills & Pinsonneault 2000).

Understanding the angular momentum distribution during the giant evolutionary phases is crucial to understand the subsequent evolutionary phases, and the spin rates of compact stellar remnants. It relies on the understanding of the interplay between convection and rotation in these very large convective envelopes. Few observational constraints are available though, due to the difficulty to measure rotation velocities for cool and slowly rotating stars exhibiting broadened spectral lines. For the case of active giants exhibiting spots at their surface, we may mention the possibility to reconstruct their latitudinal surface rotation from spectropolarimetric data (Petit et al. 2002). These type of data are however limited, and three-dimensional numerical experiments appear to be the best start to ascertain the behavior of rotating extended convective envelopes of red giant stars.

The study of turbulent convection in rotating media requires global three-dimensional simulations of convection in spherical geometry, which is both computationally very demanding and difficult to achieve due to the large number of scales that have to be represented. It is only in the last two decades that such an approach was made possible, in particular thanks to the development of massively parallel computer architectures. Most of the simulations of convection in astrophysics have been obtained in the context of the Sun (e.g., Stein & Nordlund 1998; Brun & Toomre 2002; Miesch et al. 2006, 2008), and led to a great improvement of our understanding of the solar convection. Concerning the convection realized in extended envelopes of red giant stars, three-dimensional hydrodynamical simulations of the convection in a red supergiant star (Freytag et al. 2002, hereafter FSD02; "star-in-a-box" experiment) and in asymptotic giant branch (AGB) stars (e.g., FSD02; Woodward et al. 2003, hereafter WPJ03) have been performed using fully compressible codes. These simulations do not include rotation and present as a common feature of the development of pulsations without any κ-mechanism. In the AGB star by WPJ02, the convective pattern appears to have a dipolar nature, which is also found by Kuhlen et al. (2006) in their simulation of nonrotating fully convective spheres in carbon–oxygen white dwarfs, and by Steffen & Freytag (2007, hereafter SF07) in their nonrotating "star-in-a-box" mini-sun experiment. In the FSD02 simulations of both the red supergaint and the AGB star, the convective pattern achieved consists in a small number of large cells covering the surface, the number of which increases with resolution. The surface patterns obtained by WPJ03 in their higher resolution simulations are more complex, yet resembling that found by FSD02: large warm upgoing flows surrounded by narrow cooler downgoing structures.

Recently, SF07 presented rotating "star-in-a-box" experiments. They find no striking difference in the convective pattern compared to their nonrotating case. Concerning rotation, they find a strong differential rotation, which is anti-solar in latitude, and very strong meridional circulation flows, comparable to typical convective velocities.

Other three-dimensional simulations of the convection in the upper envelope and the atmosphere of red giants have been performed in a more local approach, within the framework of realistic simulation of spectral lines formation in these stars (e.g., Chiavassa et al. 2006; Collet et al. 2007; Robinson et al. 2004). In these works, the hydrodynamical simulation of convection itself is not discussed in details. Let us finally mention the work by Herwig et al. (2007), who investigated the convection in the thermal pulses of an AGB star.

In order to further study the heat, energy, and angular momentum redistribution occurring in the extended convective envelope of slowly rotating red giant branch (RGB) stars, we have started a series of three-dimensional hydrodynamical simulations with the Anelastic Spherical Harmonics (ASH) code (Palacios & Brun 2007). In the present paper, we explore the parameter space by varying the rotation rate and the turbulence level of our RGB models and discuss in length the dynamics of such complex systems. More specifically, in Section 2 we describe our numerical approach and the three-dimensional ASH code as well as specify the one-dimensional stellar evolution model used to set our background reference state. Then, in Section 3 we discuss the properties of the turbulent convection achieved in our simulations. Large-scale flows are discussed in Section 4, with a particular emphasis on the internal rotation profile. We finally analyze and describe in Section 5 the transport of heat and angular momentum in our simulations, before summarizing the overall results and concluding in Section 6.

The simulations described here were performed using our ASH code. ASH solves the three-dimensional anelastic equations of motion in a rotating spherical geometry using a pseudospectral semi-implicit approach (e.g., Clune et al. 1999; Miesch et al. 2000; Brun et al. 2004). These equations are fully nonlinear in velocity variables and linearized in thermodynamic variables with respect to a spherically symmetric mean state. This mean state is taken to have density $\bar{\rho }$, pressure $\bar{P}$, temperature $\bar{T}$, specific entropy $\bar{S}$; perturbations about this mean state are written as ρ, P, T, and S. Conservation of mass, momentum, and energy in the rotating reference frame are therefore written as

$$\begin{equation} \mbox{\boldmath $\nabla $}\cdot (\bar{\rho }{\bf v}) = 0, \end{equation} \tag{ 1 }$$

$$\begin{eqnarray} &&\bar{\rho }\left(\!\frac{\partial {\bf v}}{\partial t}+({\bf v}\cdot \mbox{\boldmath $\nabla $}){\bf v}+2\Omega_{\rm o}\times {\bf v}\!\right)\\ &&\quad =-\mbox{\boldmath $\nabla $}P \nonumber + \rho {\bf g} - \mbox{\boldmath $\nabla $}\cdot \mbox{\boldmath $\cal D$}-[\mbox{\boldmath $\nabla $}\bar{P}-\bar{\rho }{\bf g}], \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \bar{\rho }\bar{T}\frac{\partial S}{\partial t}&=&\mbox{\boldmath $\nabla $}\cdot [\kappa _r \bar{\rho }c_p \mbox{\boldmath $\nabla $}(\bar{T}+T)+\kappa \bar{\rho }\bar{T}\mbox{\boldmath $\nabla $}(\bar{S}+S)] \\ &&\quad-\bar{\rho }\bar{T}{\bf v}\cdot \mbox{\boldmath $\nabla $}(\bar{S}+S)+2\bar{\rho }\nu \left[e_{ij}e_{ij}-1/3(\mbox{\boldmath $\nabla $}\cdot {\bf v})^2\right], \nonumber \end{eqnarray} \tag{ 3 }$$

where c_p is the specific heat at constant pressure, v = (v_r, v_θ, v_ϕ) is the local velocity in spherical geometry in the rotating frame of constant angular velocity ${\Omega _{\rm o}} = \Omega _0 {\bf \hat{e}_z}$, g is the gravitational acceleration, κ_r is the radiative diffusivity, and $\mbox{\boldmath ${\cal D}$}$ is the viscous stress tensor, with components

$$\begin{eqnarray} &&{\cal D}_{ij}=-2\bar{\rho }\nu [e_{ij}-1/3(\mbox{\boldmath $\nabla $}\cdot {\bf v})\delta _{ij}], \end{eqnarray} \tag{ 4 }$$

where e_ij is the strain rate tensor. Here, ν and κ are effective eddy diffusivities for vorticity and entropy. To close the set of equations, linearized relations for the thermodynamic fluctuations are taken as

$$\begin{equation} \frac{\rho }{\bar{\rho }}=\frac{P}{\bar{P}}-\frac{T}{\bar{T}}=\frac{P}{\gamma \bar{P}} -\frac{S}{c_p}, \end{equation} \tag{ 5 }$$

assuming the ideal gas law

$$\begin{equation} \bar{P}={\cal R} \bar{\rho }\bar{T}, \end{equation} \tag{ 6 }$$

where ${\cal R}$ is the gas constant. The effects of compressibility on the convection are taken into account by means of the anelastic approximation, which filters out sound waves that would otherwise severely limit the time steps allowed by the simulation.

For boundary conditions at the top and bottom of the domain, we impose

• 1.  
  impenetrable walls:

  $$v_r=0|_{r=r_{\rm bot},r_{\rm top}},$$
• 2.  
  stress free conditions:

  $$\frac{\partial }{\partial r}\left(\frac{v_{\theta }}{r}\right)=\frac{\partial }{\partial r}\left(\frac{v_{\phi }}{r}\right)=0|_{r=r_{\rm bot},r_{\rm top}},$$
• 3.  
  and constant entropy gradient:

  $$\frac{\partial \bar{S}}{\partial r}=cst|_{r=r_{\rm bot},r_{\rm top}}.$$

Convection in stellar environments occurs over a large range of scales. Numerical simulations cannot, with present computing technology, consider all these scales simultaneously. We therefore seek to resolve the largest scales of the nonlinear flow, which we think are likely to be the dominant players in establishing differential rotation and other mean properties in convective envelopes. We do so within a large-eddy simulation (LES) formulation, which explicitly follows larger scale flows while employing subgrid-scale (SGS) descriptions for the effects of the unresolved motions. Here, those unresolved motions are treated as enhancements to the viscosity and thermal diffusivity (ν and κ), which are thus effective eddy viscosities and diffusivities. For simplicity, we have taken these to be functions of radius alone, and to scale as the inverse of the square root of the mean density. We emphasize that currently tractable simulations are still many orders of magnitude away in parameter space from the highly turbulent conditions likely to be found in real stellar convection zones. These LESs should therefore be viewed only as indicators of the properties of the real flows. We are however encouraged by the success that similar simulations (e.g., Miesch et al. 2000; Elliott et al. 2000; Brun & Toomre 2002; Miesch et al. 2006) have enjoyed in matching the detailed observational constraints for the differential rotation within the solar convection zone provided by helioseismology.

Thermodynamic variables within ASH are expanded in spherical harmonics Y^m_ℓ(θ, ϕ) in the horizontal directions and in Chebyshev polynomials T_n(r) in the radial direction. Spatial resolution is thus uniform everywhere on a sphere when a complete set of spherical harmonics of degree ℓ is used, retaining all azimuthal orders m. We truncate our expansion at degree ℓ = ℓ_max, which is related to the number of latitudinal mesh points N_θ (here ℓ_max = (2N_θ − 1)/3), take N_ϕ = 2N_θ azimuthal mesh points, and utilize N_r collocation points for the projection onto the Chebyshev polynomials. We have considered various grid resolutions N_r × N_θ × N_ϕ depending on the degree of turbulence of the model (see Table 1). A semi-implicit, second-order Crank-Nicolson scheme is used in determining the time evolution of the linear terms, whereas an explicit second-order Adams-Bashforth scheme is employed for the advective and Coriolis terms. The ASH code has been optimized to run efficiently on massively parallel supercomputers such as the IBM SP-6 or SGI Altix, and has demonstrated excellent scalability on such machines (see Clune et al. 1999 and Brun et al. 2004 for more details).

Table 1. Parameters for the Four Simulations

Case	RG1	RG1t	RG2	RG2t
Nr, Nθ, Nϕ	257, 256, 512	257, 512, 1024	257, 256, 512	257, 768, 1536
Ra ...	8.25 × 105	3.18 × 106	7.45 ×105	2.69 × 106
Ta ...	1.32 × 106	4.89 × 106	5.44 × 104	1.96 × 105
Pr ...	1	1	1	1
Roc ...	0.791	0.806	3.703	3.703
νtop ...	1.2 × 1015	5 × 1014	1.2 × 1015	5 × 1014
κtop ...	1.2 × 1015	5 ×1014	1.2 × 1015	5 × 1014
$\tilde{Re^{\prime }}$ ...	256	623	372	882
$\tilde{Pe^{\prime }}$ ...	256	623	372	882
$\tilde{Ro^{\prime }}$ ...	0.281	0.284	2.03	2.01

Notes. All simulations have an inner radius r_bot = 1.36 × 10¹¹ cm, an outer radius r_top = 1.36 × 10¹² cm, with L = 1.22 × 10¹² cm (≃17.6 R_☉). The number of radial, latitudinal, and longitudinal mesh points are N_r, N_θ, and N_ϕ, respectively. The higher degree of turbulence in cases RG1t and RG2t was obtained by maintaining the Prandtl number at 1 and lowering both the eddy viscosity ν and eddy diffusivity κ at the top edge of the domain. ν and κ are given in units of cm s⁻². The characteristic numbers evaluated at mid-layer depth are the Rayleigh number R_a = (∂ρ/∂S)ΔSgL³/ρνκ, the Taylor number T_a = 4Ω²₀L⁴/ν², the Prandtl number P_r = ν/κ, the convective Rossby number $R_{oc} = \sqrt{R_a / T_a P_r}$, the rms Reynolds number $\tilde{R_e^{\prime }} = \tilde{\upsilon }^{\prime } L / \nu$, the rms Péclet nulber $\tilde{P_e^{\prime }} = \tilde{R_e^{\prime }} P_r$, and the rms Rossby number $\tilde{R_o^{\prime }} = \tilde{\omega }^{\prime }/2 \Omega _O \approx \tilde{\upsilon }^{\prime }/2 \Omega _0 L$. We use the rms convective velocity $\tilde{\upsilon }^{\prime }$ at mid-depth given in Table 2 to compute these quantities and the following mid-depth viscosity of 8.4810¹⁴ (3.5310¹⁴) cm s⁻² for, respectively, the laminar (turbulent) cases.

Download table as:  ASCIITypeset image

The models considered here are intended to be simplified descriptions of the bulk of the deep extended convective envelope of an evolved (RGB phase) 0.8 M_☉ star. We do not consider in this study the presence of an overshooting layer (compare Miesch et al. 2000; Browning et al. 2006; Brun 2009). Contact is made with a one-dimensional stellar model (at an age of 11 Gyr) for the initial conditions, which adopts realistic values for the radiative opacity, density, and temperature Palacios et al. (2006). This one-dimensional model was computed with the STAREVOL stellar evolution code (see Siess et al. (2000) and Siess (2006) for more details on the numerical methods and physical ingredients adopted). Convection is computed within STAREVOL using a classical mixing length treatment with a parameter α_MLT = 1.75, and the convective limits are determined according to the Schwarzschild criterion for the convective instability. The luminosity of the modeled star L_* is 425 L_☉, its radius R_* is about 40 R_☉. Our simplified three-dimensional simulations were initialized using the radial profiles of gravity g, radiative diffusivity κ_rad, and the mean density ${\bar \rho }$ of the one-dimensional model along with a prescribed mean entropy gradient $d\bar{S}/dr$ as the starting points for an iterative Newton–Raphson solution for the hydrostatic balance and for the gradients of the thermodynamic variables. The mean temperature $\bar{T}$ is then deduced from Equation (6). This technique yields background reference profiles in reasonable agreement with the thermally relaxed one-dimensional stellar model as can be seen in Figure 1, and we are thus confident that the background state of our simulations is close to their final relaxed state. We have built four different models that all share this one-dimensional structure, and list their most important parameters in Table 1. Cases RG1 and RG1t are computed with an initial bulk rotation rate Ω₀ = Ω_☉/10 = 2.6 × 10⁻⁷ rad s⁻¹ (or a rotation period P_rot of ∼280 days), while we have used Ω₀ = Ω_☉/50 = 5.2 × 10⁻⁸ rad s⁻¹ (or P_rot∼1400 days) for cases RG2 and RG2t. Cases RG1t and RG2t are more turbulent versions of cases RG1 and RG2, respectively (that have been published in Palacios & Brun 2007). Case RG2 has been evolved from scratch, but we have also checked that reducing the rotation rate of case RG1 down to the value adopted for case RG2 leads to the same configuration in terms of convection and rotation. For this reason, cases RG1t and RG2t have been, respectively, evolved from cases RG1 and RG2 by progressively reducing the effective values of κ_top and ν_top.

Zoom In Zoom Out Reset image size

Although the large rotation adopted for cases RG1 (and RG1t) corresponds to a surface equatorial velocity much larger than that inferred from observations, these two cases could be viewed as less evolved counterparts of cases RG2 and RG2t, respectively, since the surface rotation will decrease as the radius increases during the ascent of the RGB due to angular momentum conservation.⁴

All the models presented here have been computed over several rotation periods and about one hundred convective overturning times (or about 20–30 years of physical temporal evolution). This long temporal evolution is required in order to reach a statistically stationary state and equilibrated balances of energy, heat, and angular momentum within the convective shell. Even though the Kelvin–Helmotz time τ_KH ∼ GM²_*/R_*L_* is about 1200 years for our RGB models, we believe that the properties of the models discussed in the subsequent sections correspond to a mature and well equilibrated state. Of course, we cannot rule out a very slow evolution of the model when integrated over "very long period of time" (as done in Chan 2007), but given the realistic one-dimensional structure used for the background state and the small Mach number (i.e., small temperature and density fluctuations with respect to the mean background values) realized in our simulations we are convinced that our results are robust and a fair description of the nonlinear dynamics present inside a rotating convective envelope. In order to reach such an equilibrated state, we have use a large amount of CPU-time on massively parallel computers (of order 600,000 hr and more per simulation), in particular with cases RG2 and RG2t, due to their long rotation period.

Given reasonable computing resources, the ASH code does not allow at present to compute the full convective envelope of a RGB star, in which the density varies by more than 4 orders of magnitude, but only a portion of it (see Figure 1). We will thus assume in all the cases presented in this paper a density variation of 2 orders of magnitude. This choice translates into computing the inner half (e.g., from r_bot = 0.05 R_* up to r_top = 0.5 R_*) of the extended convective envelope of a 0.8 M_☉ RGB star. We concentrate on the inner part of the convective envelope, since near the surface the Mach number would be too large for the anelastic formalism to remain valid. Moreover, this deeper part of the extended convective envelope is particularly interesting from the point of view of one-dimensional stellar evolution models, because it is the region that may gradually be incorporated into the underlying radiative zone, and thus establish the link between nuclearly processed regions and the chemically homogeneous convective envelope. Indeed, during the RGB ascent, the angular momentum of the innermost part of the convective envelope will be transferred to the underlying radiative interior as the convective envelope retreats in mass after the first dredge-up Clayton (1968). This will influence the rotation-induced transport of both angular momentum and chemicals in the radiative zone.

The strong radial enthalpy flux found in our simulations results from the correlation between radial velocity and temperature fluctuations shown for RG1t in Figure 4 (left and middle panels), where we display horizontal maps of v_r and T (and enstrophy).

Zoom In Zoom Out Reset image size

The convective patterns realized in our four models vary significantly depending on the adopted reference rotation frame Ω₀. Cases RG1 and RG1t rotating at one tenth the solar rate possess about 10 downflow lanes over the uppermost layers. These relatively cool downflow lanes are interconnected, forming an intricate network surrounding warm broad upflows. They appear respectively in blue (cool downflows) and reddish (warm upflows) in Figure 5, where we present a three-dimensional rendering of the temperature and velocity in the model RG1t. They present some systematic alignment with the north–south direction, such property being quite obvious for the least turbulent case RG1. Large velocity and temperature fluctuations are observed throughout the domain (as also evident in Figure 8 for the velocity where we plot the radial profile of the rms velocity components), a direct consequence of the high luminosity of the modeled star. Typical values for these quantities are ±300 K and ±1000 m s⁻¹ in most of the computational domain, in keeping with the necessity of convection to carry outward the large imposed radiative flux. Comparatively to the simulations of the solar convection envelope, these values are at least 1 order of magnitude larger (Brun & Toomre 2002). The patterns described for the velocity fluctuations are also visible in the temperature fluctuations. The contrast between hotter and cooler parts is however weaker than that existing between upflows and downflows in the velocity maps. Nonetheless, due to our choice of a Prandtl number Pr = 1, they exhibit the same fine structures. Those structures are particularly obvious in the enstrophy map (right panel of Figure 4), with downflow lanes being part of an intricate network and possessing the highest enstrophy concentration near the interstices. This is a direct consequence of the fact that as flows converge toward the downdraft they acquire intense cyclonicity, whereas upflows expand as they rise and lose their cyclonic character. The vortex tubes seen in the enstrophy map actually come by pair of opposite gyres. They are the weaker and broader equivalent to the intense "cold spinners," seen in solar convection simulations (Miesch et al. 2008). In Figure 6, we plot the radial vorticity map at the same radius as in Figure 4 for cases RG1t and RG2t. For case RG1t, the radial vorticity presents near the top of the domain a clear dominant negative (positive) value in the northern (southern) hemisphere, respectively, a direct consequence of the action of the stronger Coriolis force in this more rapidly rotating model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Turning to the slowly rotating cases RG2 and RG2t in Figure 4, we clearly note a pattern change. This is mostly due to a weaker influence of rotation on convection (as illustrated in Table 1 by their larger convective Rossby number with respect to cases RG1 and RG1t) leading to more vigorous and isotropic flow. The downflow lanes show much less systematic alignment with respect to the north–south direction. The weaker influence of rotation also leads to a shift toward a smaller number of downflow lanes over the surface, with here only between 4 and 6 lanes covering the top of the computational domain. Further a strong dipole in temperature is observed, with two separate large regions of the surface being respectively predominantly cool and warm (cf. Section 5.1 for more details). As a consequence of this thermal structuration of the convective flow, a large coherent circulation sets in between these two regions, yielding two zones that are respectively filled mostly with converging narrow downflows on the one hand and diverging broad upflows on the other hand. This dipole in temperature and velocity persists in depth as can be seen from the equatorial slices presented in Figure 7 for case RG2t. We clearly note the large inflow in the top part of the figure with the corresponding outflow on the other side of the central region. This large-scale flow goes through the central part without any difficulty connecting efficiently regions of the domain that are far apart. The temperature patterns are more finely structured but overall there is good correlation between respectively the cold (warm) and inflow (upflow) regions. This bipolar structure of convection is a very interesting property of our slowly rotating case and points toward the results previously obtained by Porter & Woodward (2000), WPJ03, and SF07 in nonrotating hydrodynamical simulations. Indeed these authors report in their nonrotating giant stars the presence of dominant large-scale convective cells over the star's top domain embedded in a dipolar flow crossing the whole convection domain. It is interesting to note that Chandrasekhar (1961) in his linear axisymmetric study of the onset of convection in a spherical shell, has demonstrated that the ℓ = 1 mode is for small aspect ratio, one of the most unstable modes to excite. This is mostly due to a geometrical effect with convection rolls having more space to develop for low-aspect ratio (or full sphere) than for large-aspect ratio (thin shell), thus leading to larger physical structures and favoring the low ℓ modes. We have also run a nonrotating simulation for the very same initial structure as that used for the four cases discussed in the present paper, and this dipolar structure also shows up clearly in that simulation (which is not presented here). This dipole can only be guessed in the three-dimensional equatorial cut presented in Figure 5 for model RG1t.

Zoom In Zoom Out Reset image size

As for cases RG1 and RG1t, the more slowly rotating models possess some enstrophy, with the downflows having the largest values due to their cyclonic character. The radial vorticity map for case RG2t shown in Figure 6 does not possess a dominant sign per hemisphere, contrary to case RG1t even though there seems to be some weak tendency for a negative/positive orientation (sign) to dominate near respectively the warm/cool regions. The less systematic orientation in each hemisphere of the radial vorticity, again, is a consequence of the larger convective Rossby number of this slowly rotating model, and of the weaker hemispherical anisotropy imposed by the Coriolis force of the flows.

We display on Figure 8 the radial profile of the total rms velocity and its radial, latitudinal and longitudinal components for cases RG1, RG1t, RG2, and RG2t. Typical values are of 1–3 km s⁻¹ to be compared with a sound speed evolving in the range 25–80 km s⁻¹ in the models, thus a posteriori confirming the validity of the anelastic approximation (see also Figure 1). We note that the radial velocity dominates in the bulk of the computational domain roughly between r = 0.07 R_* and r = 0.3 R_*. This is due to the continuous acceleration of the convective plumes over the whole domain depth, with $\tilde{v}_r$ becoming maximal near the bottom. In the more slowly rotating cases RG2 and RG2t, $\tilde{v}_r$ is even faster due again to the reduced influence of rotation on the vigor of the convective motions and possess a pronounced bell-like profile. It is indeed known since (Chandrasekhar 1961) that rotation tends to inhibit convection by increasing the critical Rayleigh number, thus leading to a decrease of the supercriticality of the models for a given Ra. Due to our rigid (impenetrable) boundary conditions, the radial velocity is forced to vanish at the domain's edges explaining its peculiar shape. Independent of the adopted Reynolds number, the horizontal rms velocity components have very similar amplitudes. $\tilde{v}_{\phi }$ profiles slightly differ between the slow and faster models, with the slower models having a steeper profile in the outer parts of the computational domain. This is due to the shellular rotation existing in the slowly rotating models as will be discussed in the following sections. In all our four cases, the horizontal components dominate toward the top of domain. At the bottom boundary, the horizontal flow amplitudes increase significantly possibly due to local angular momentum conservation as the flow converge toward the center. This is even more evident for cases RG2 and RG2t, where the inner shells have much faster horizontal velocities due to the shellular rotation.

Zoom In Zoom Out Reset image size

The temporal evolution of the convective patterns for case RG1t is shown in Figure 9, which displays four sequences of images of the radial velocity at three different depths over a period of about 96 days. The first row shows the equatorial plane at each timeshot, and the different depths represented in Mollweide projections are indicated with dotted, dashed, and dashed-dotted black circles for 19 R_☉ (second row), 11 R_☉ (third row), and 3 R_☉ (forth row), respectively. As can be seen from the equatorial cuts and the first row of Mollweide projections at 19 R_☉, the convective structures at the top of the computational domain move from east to west (retrograde, i.e., clockwise). Deeper in the convective region, structures have a prograde (counterclockwise) displacement. The equatorial cuts indicate that the up and downflows maintain their coherence over most of the computational domain. This is far less evident from the Mollweide projections due to the large geometrical factor between both edges of the domain, as well as to the intrinsic tilt of the structures.

Zoom In Zoom Out Reset image size

The evolution presented here corresponds to about a third of the rotation period of case RG1t (see Table 1). The global convective overturn given by

$$\begin{equation*} \tau _c \approx \int _{R_b}^{R_t} \frac{dr}{\upsilon _{rms}} \end{equation*}$$

is of 75 days in this case, which is much shorter than the rotation period. For this reason, we completely lose track of the structures when considering a complete rotation. This is also the case for our model RG1. For the slowly rotating cases RG2 and RG2t, this is even more striking, since the convective overturn is again of about 80 days, while the rotation period is 1398 days.

As discussed in the introduction, the rotation profile established in the convective envelope of RGB stars is particularly important to the understanding of their global evolution. We believe that our models can give us a good hint of such profiles as a function of depth and latitude.

In Figure 10, we display contours of the longitudinal and temporal average of the angular velocity realized in our four models over a rotation period (cf. Table 1) along with radial cuts at indicated latitudes. For the cases rotating at a tenth of the solar rate, we note a large angular velocity contrast both in radius and latitude, with similar amplitudes in both directions. The central regions are found to rotate extremely fast in a prograde sense, whereas the uppermost layers have a retrograde rotation (but remain prograde when considering the bulk rotation Ω₀). The differential rotation is slower at the equator than in the polar regions, thus yielding an anti-solar rotation profile. Such a profile has also been found by SF07 in their simulation of rotating convection as mentioned in Section 1. Let us also note that anti-solar differential rotation has also been reported at the surface of giant stars (e.g., Weber et al. 2005; Weber 2007). The most striking properties of the differential rotation achieved in these models is that it is almost invariant along the z-axis (i.e., which coincides with the rotation axis). This is a well-known dynamical property of rotating fluids, in which the fluid velocity tends to be uniform along lines parallel to the rotation axis (Pedlosky 1987). The Taylor–Proudman constrain on rotating flows and the way of potentially braking it have been extensively studied in the context of the conical rotation profile of the Sun (e.g., Kitchatinov & Rüdiger 1995; Durney 1999; Brun & Toomre 2002; Rempel 2005; Miesch et al. 2006). In these papers, it has been shown that along with the Reynolds stresses, the baroclinic term involving latitudinal gradient of the temperature and entropy fluctuations plays a significant role in shaping the solar differential rotation. The anisotropic transport of heat by convection or the influence of a tachocline at the base of the convection zone both seem to contribute to break the cylindricity (Brun & Rempel 2008). We delay to Section 5.3 the quantitative discussion of all these effects in our RGB models.

Zoom In Zoom Out Reset image size

In Figure 11, we show latitudinal profiles of Ω for cases RG1 and RG2 at four different instants along with the average formed over the whole temporal interval sampled (solid line). The various curves correspond to temporal averages formed over either successive one rotation periods for RG1 or over temporal intervals of hundreds of days compatible with the period of time for which RG2 possesses one or multiple meridional cells (see below Section 4.2). We note the relatively good stability of the profiles found for RG1 over several convective overturning times and rotation periods (except near the polar regions due to the small level arm there making it difficult to form meaningful averages), thus confirming the well equilibrated state reached by this case. For case RG2 the situation is less obvious, with clear departure from the mean at given instant mostly because this case undergoes a large dynamical oscillation due to its slow rotation rate that we discuss in length in Section 4.2. Note however the much smaller range of variation as a function of latitude found in that case (30 nHz) with respect to case RG1 (200 nHz). This is due to the peculiar angular velocity profile realized in this model.

Zoom In Zoom Out Reset image size

Returning indeed to Figure 10 (bottom rows) showing the angular velocity in the two cases rotating at a fiftieth of the solar rate, we note a large angular velocity contrast in radius solely due to the striking shellular state (i.e., angular velocity uniform on spherical surfaces) achieved by the simulations. Further, the rotation near the top edge of the computational domain is found to be retrograde in an absolute way, i.e., even taking into account the reference frame rotation. Given the large differences between cases RG1 (RG1t) and RG2 (RG2t), we have run verification models to be sure of the robustness of the results. For instance, we have progressively decelerated case RG1 to 1/50th the solar rate to check if we recover the shellular state achieved in case RG2. We find that when reaching a value of the convective Rossby number larger than about 1 (i.e., for a rotation around 1/15th solar), the influence of rotation on convection is reduced (as expected from earlier studies; Glatzmaier & Gilman 1982; Browning et al. 2004; Ballot et al. 2007; Brown et al. 2008). Thus, the cylindrical profile of the differential rotation observed in cases RG1 and RG1t is lost, the flow being less constrained to be quasi two-dimensional along the rotation axis. This has for direct consequence that such slowly rotating models do not show any angular velocity contrast in latitude.

Figure 12 presents the mean radial profile for the specific angular momentum obtained for each of our four models after averaging the angular velocity over longitudes and latitudes, and in time over one rotation period (except for case RG2t, averaged over all the iterations available). Also plotted on this view graph are the expected specific angular momentum profiles associated with uniform angular velocity profiles with Ω₀ = 0.1Ω_☉ (long-dashed triple dotted line) and Ω₀ = 0.02Ω_☉ (long-dashed line). The profiles achieved by models RG1 and RG1t are very similar both in shape and amplitude. Similarly, the profiles obtained for models RG2 and RG2t are also almost identical. This indicates that the distribution of specific angular momentum (and of the mean radial angular velocity) achieved in the simulations is not influenced by the turbulent state (e.g., Reynolds number) of the simulations up to the values that we have been able to compute. This is a quite different behavior compared to the solar convection, where lowering the diffusivities while keeping P_r constant in order to increase R_e by a factor of 2 leads to some modification of the angular velocity profile (see Figure 4 of Brun & Toomre 2002). On the other hand, the specific angular momentum distribution achieved in our convective shell strongly depends on the bulk rotation considered. For models RG1 and RG1t with a bulk rotation a tenth of the solar value, the profile presents a positive slope throughout the computational domain. In these models, due to the cylindrical profile of Ω (see Figure 10), the radial contrast at high latitudes is smaller than at the equator so that these latitudes contribute to soften the mean radial profile presented in Figure 12.

Zoom In Zoom Out Reset image size

The specific angular momentum profile for the slower cases RG2 and RG2t is radically different, with a change of slope below r ≈ 0.2 R_* and a negative slope in the outer part of the computational domain. This is directly related to the shellular rotation existing in these simulations. At small radii, the angular velocity is large, and although its radial mean profile has a negative slope, the contribution of the increasing radius (i.e., the specific angular momentum scales as Ωr²) maintains a positive slope for the specific angular momentum j. Moving toward the top of the domain, the mean radial angular velocity dramatically drops at all latitudes in these slowly rotating simulations (the total variation of Ω(r) at all latitudes is of more than 2 orders of magnitude; see Figure 10), decreasing faster than the increase of the radius, and leading to the observed change of slope. Near the top of the convective shell the absolute retrograde rotation seen in Figure 10 leads to a second change of slope in the mean radial profile of specific angular momentum.

Let us finally note that none of our simulations approach the extreme cases of uniform mean radial specific angular momentum or uniform mean radial angular velocity, as have been assumed in the modeling of angular momentum transport in the convective envelopes in one-dimensional stellar evolution models.

Another important large-scale flow established in rotating convective envelopes is the mean (axisymmetric) meridional circulation (i.e., the flow in the r–θ plane). This flow is maintained by small imbalance between latitudinal pressure gradient, Coriolis force acting on the differential rotation, Reynolds stresses, and buoyancy forces in the purely hydrodynamical case (e.g., Miesch 2005; Brun & Rempel 2008). In the Sun, this relatively weak flow (with respect to the solar differential rotation) is thought to play an important role in setting the solar cycle (Jouve & Brun 2007). It also plays important role for the redistribution of angular momentum, even though it only contains about 0.5% of the total kinetic energy (Brun & Toomre 2002). By contrast, in our RGB simulations we find that the kinetic energy contained in the meridional circulation (MCKE) accounts for about 8% of the total kinetic energy (see Table 2 and Figure 2) and thus cannot be neglected in the energy balance of our convective envelope. We display in Figure 13 the axisymmetric streamlines along with the radial and latitudinal components of the mean (axisymmetric) velocity for the four cases analyzed in the present paper. Again the rotation plays a significant role in shaping the flow properties, leading to a clear difference in the realized meridional circulation profiles. For cases RG1 and RG1t, we find a very stable profile consisting of one large single cell per hemisphere. Each cell is directed toward the pole at the top of the domain and toward the equator deeper down. Near the star's rotation axis the flow is mostly radial (vertical), and near the equator it is directed perpendicular to that axis. The contour maps reveal that the axisymmetric 〈v_r〉 and 〈v_θ〉 components (with 〈 〉 denoting longitudinal average) are of the same order of magnitude, with 〈v_θ〉 being antisymmetric with respect to the equator and 〈v_r〉 symmetric. For cases RG2 and RG2t, the picture is completely different. Due to the large dipole in the temperature maps (cf. Figure 4), the meridional circulation consists of one large single cell spanning the whole shell. The direction of this flow is different for models RG2 and RG2t, and depends on the latitudinal fluxes (see below). In both cases, there is no radial transport in the equatorial region. The amplitude of 〈v_r〉 and 〈v_θ〉 are larger than in the faster rotating cases (1000 m s⁻¹ versus 700 m s⁻¹). Due to its peculiar topology, the radial and latitudinal components of the mean velocity have opposite symmetries with respect to the equator and cases RG1 and RG1t.

Zoom In Zoom Out Reset image size

Table 2. Representative Velocities, Energy Contents, and Differential Rotation Contrast

Model	$\tilde{v}_r$	$\tilde{v}_{\theta }$	$\tilde{v}_{\phi }$	$\tilde{v}_{\phi }^{\prime }$	$\tilde{v}$	$\tilde{v}^{\prime }$	KE	DRKE	CKE	MCKE	ΔΩ
RG1	1348	1015	1406	848	2202	1780	1.3 × 106	4.7 × 105 (36%)	7.1 × 105 (55%)	1.2 × 105 (9%)	78.9
RG1t	1370	1040	1456	914	2263	1803	1.3 × 106	4.6 × 105 (36%)	7.2 × 105 (56%)	1.1 × 105 (8%)	64.7
RG2	1905	1097	1511	1455	2670	2586	2.37 × 106	1.5 × 105 (6%)	2.11 × 106 (90%)	1.0 × 105 (4%)	4.2
RG2t	1939	1188	1380	1343	2666	2552	2.15 × 106	1.4 × 105(6.5%)	1.85 × 106 (86%)	1.6 × 105 (7.5%)	4.4

Notes. Temporal averages of the rms velocity $\tilde{v}$, of the rms components $\tilde{v}_r$, $\tilde{v}_{\theta }$, $\tilde{v}_{\phi }$, and of fluctuating velocities $\tilde{v}^{\prime }$ and $\tilde{v}_{\phi }^{\prime }$ (axisymmetric components are removed) are estimated at mid-layer depth, all expressed in units of m s⁻¹. For cases RG2 and RG2t that undergo a long-term oscillation of KE, we have chosen to form the temporal average over a period representative of the simulations by avoiding peak (min/max) values. Also listed are the time average over the complete volume of the total kinetic energy KE and that associated with the (axisymmetric) differential rotation DRKE, the (axisymmetric) meridional circulation MCKE, and the nonaxisymmetric convection CKE, all in units of erg cm⁻³. The latitudinal contrasts between 0° and 60° of angular frequencies ΔΩ, quoted in nHz, are computed at the top of the domain.

Download table as:  ASCIITypeset image

We also find that depending on the period over which the time average is performed, the pattern of the meridional circulation will present one or two cells in models RG2 and RG2t.

In order to illustrate the rich dynamics operating in the slowly rotating case RG2, we display in Figure 14 the temporal evolution of the meridional streamlines along with the temporal traces of the kinetic energy densities (KE, DRKE, CKE, and MCKE) over an interval of 2500 day near the end of the total evolution of the simulation. This allows us to identify the phases in which model RG2 possesses more than one large single cell. It turns out that the appearance of a second cell in the southern hemisphere is linked to a weakening of MCKE with respect to DRKE. When MCKE becomes smaller than DRKE we find a multicellular meridional circulation profile similar to that achieved in cases RG1 and RG1t, for which DRKE is always larger than MCKE. On the contrary when MCKE is larger than DRKE, situation that occurs often and lasts over a much longer periods of time as illustrated in the figure (over the ranges a and c), the meridional flow possesses one big clockwise cell. We believe that the reduced rotational constraint in the slowly rotating case RG2 is at the origin of such an intriguing behavior by allowing the meridional circulation to evolve freely. Further the oscillating meridional circulation patterns are likely a consequence of the convective structures that are clearly influenced by the strength and the location of the self-established temperature dipole (Figure 14).

Zoom In Zoom Out Reset image size

Analyzing the kinetic energy budget contained in the turbulent convective motions is very useful to understand key features of our simulations, such as convection and large-scale flows, and how they vary with the model's parameters. As we have seen in Figure 3, the large stratification present in our models leads to an inward kinetic energy flux, with fast concentrated downflows overwhelming the slow broad upflows in transporting the kinetic energy. Indeed, the kinetic energy flux F_k scales as v_rv² so it is sensitive to the sign of the radial velocity component v_r, and as v_r in downflows is larger than in upflows, it yields when averaging over horizontal surfaces this asymmetric (downward) transport of the kinetic energy. The slowly rotating cases possess a more vigorous convection, with faster rms velocities (about a factor 1.3 higher) as indicated in Table 2 (and also Figure 8). Associated with these faster flows we also find that cases RG2 and RG2t possess a larger inward directed kinetic energy flux than their faster rotating counterparts (RG1 and RG1t). Converted to luminosity the kinetic energy flux reaches −120% (compared to −50%) of the total luminosity leading to an extremely large convective luminosity of 220% (compared to 140%; see Figure 3). Here, again the stabilizing effect of rotation on convection plays a central role in limiting the amplitude of the flows and thus the strength of F_k in cases RG1 and RG1t. Even if we compare case RG1 to case RG2 which have the same Reynolds number, the asymmetry between upflows and downflows is stronger in the slowly rotating case RG2.

To get a better insight on how the kinetic energy is redistributed among the various axisymmetric and nonaxisymmetric motions, we again split the KE into it components DRKE, MCKE, and CKE. We have already partly discussed the different behavior of these components in Section 4.2 via Figure 13. We report in Table 2 the volume and time-averaged values of KE and its components for the four cases discussed in the paper. These values reveal a strong difference between cases RG1 (RG1t) and RG2 (RG2t). The kinetic energy associated with the differential rotation DRKE represents a third of the total kinetic energy KE in cases RG1 and RG1t but only 6% in models RG2 and RG2t, where the kinetic energy density is essentially concentrated in the (nonaxisymmetric) convection itself. This is again due to the mild influence of rotation on convection in these simulations with slow bulk rotation. As the bulk rotation rate increases, more kinetic energy is diverted to longitudinal flows (i.e., differential rotation), and much less so in the meridional flows. It is striking to note that MCKE in slowly rotating case is close to DRKE and leads to fascinating modulation of the dynamics as seen in Section 4.2 and Figure 14. In cases RG2 and RG2t, the contribution of the mean flows to the total kinetic energy budget is noticeably weak and the convective motions are mostly nonaxisymmetric. Further these two cases exhibit an oscillation of their total kinetic energy that lasts about a rotation period. The difference between evaluating the rms velocities over a period of time for which KE is minimum or maximum is or order 3%–5%. Since case RG2t is slowly emerging from a pronounced deep of KE this explains why in Table 2, some of the rms velocities for this case are smaller than in RG2. Nevertheless given the smaller kinetic viscosity used in RG2t, overall this case has a larger Reynolds number (cf. Table 1) and is effectively more turbulent than RG2.

Another useful information on the turbulent properties of convection can be retrieved by plotting the kinetic energy spectra as a function of the spherical harmonic degree ℓ as was done in Figure 15. In this figure, we clearly see that the difference of amplitude between the large (low ℓ degree) and the small scales is important. It involves many orders of magnitude, thus confirming the well-resolved character of our simulations. For the fast rotating cases, the profiles present a plateau at low ℓ degree and the maximum of the energy distribution is located around ℓ = 10. None of these features exist for cases RG2 and RG2t. Quite interestingly these more slowly rotating cases display a conspicuous peak at ℓ = 1 that dominates all the other scales. This points toward the dipole in temperature and velocity seen in the Mollweide maps of Figure 4, in the three-dimensional rendering of Figure 5 and in the equatorial slices of Figure 7. As previously mentioned, this is expected since in his linear study of the onset of convection in spherical shells of various thickness, Chandrasekhar (1961) has demonstrated that the ℓ = 1 mode is among the easiest mode to excite when the topology is close to being a full sphere with an aspect ratio r_min/r_max equal or close to zero. It is important to realize that our simulations are fully nonlinear as demonstrated by the energy spectra, but that nevertheless the ℓ = 1 mode remains a dominant convection mode since our choice in our models of a small aspect ratio (i.e., r_min/r_max = 0.1) favors low ℓ modes.

Zoom In Zoom Out Reset image size

Convection not only transports heat but also redistributes angular momentum. We seek here to understand the mechanisms responsible for the transport of angular momentum within our rotating convective shells yielding the differential rotation discussed in Section 4.1. With our choice of boundary conditions (see Section 2.1), no net torque is applied to the shell and the total angular momentum must be conserved. We indeed find that in our four simulations the total angular momentum averaged over the whole volume is conserved to within 10⁻⁷.

The transport of angular momentum may be assessed by considering the mean radial and latitudinal angular momentum fluxes ${\cal F}_r$ and ${\cal F}_\theta$. Following Elliott et al. (2000) and Brun & Toomre (2002), the expression of these fluxes is extracted from the ϕ component of the momentum equation after averaging in time and longitude (denoted by the symbol $\hat{\ }$; see also Palacios & Brun 2007):

$$\begin{equation} {\cal F}_r=\hat{\rho }r\sin \theta \left[ \underbrace{-\nu r\frac{\partial }{\partial r}\left(\frac{\hat{v}_{\phi }}{r}\right)}_{{\cal F}_{r,V}}+\underbrace{\widehat{v_{r}^{^{\prime }} v_{\phi }^{^{\prime }}}}_{{\cal F}_{r,R}}+\underbrace{\hat{v}_r(\hat{v}_{\phi }+\Omega _0 r\sin \theta)}_{{\cal F}_{r,MC}} \right] \end{equation} \tag{ 7 }$$

$$\begin{eqnarray} {\cal F}_{\theta }&=&\hat{\rho }r\sin \theta \left[\underbrace{-\nu \frac{\sin \theta }{r}\frac{\partial }{\partial \theta }\left(\frac{\hat{v}_{\phi }}{\sin \theta }\right)}_{{\cal F}_{\theta,V}}+\underbrace{\widehat{v_{\theta }^{^{\prime }} v_{\phi }^{^{\prime }}}}_{{\cal F}_{\theta,R}}\right.\nonumber\\ &&\left.+\underbrace{\hat{v}_{\theta }(\hat{v}_{\phi }+\Omega _0 r\sin \theta)}_{{\cal F}_{\theta,MC}} \right]. \end{eqnarray} \tag{ 8 }$$

${\cal F}_{r,V}$ (resp. ${\cal F}_{\theta,V}$) is the flux associated with 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. As was done in Brun & Toomre (2002), 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:

$$\begin{eqnarray} I_r(r)&=&\int _0^{\pi } {\cal F}_r(r,\theta) \, r^2 \sin \theta \, d\theta, \nonumber\\ I_{\theta }(\theta)&=&\int _{r_{\rm bot}}^{r_{\rm top}} {\cal F}_{\theta }(r,\theta) \, r \sin \theta \, dr \,. \end{eqnarray} \tag{ 9 }$$

These integrated fluxes are presented in Figure 16 for our cases RG1 and RG2, and have been averaged over 8.5 rotation periods. For simplicity we drop the letter I when discussing the individual contribution of the flux.

Zoom In Zoom Out Reset image size

For the faster case RG1, the radial flux ${\cal F}_{r,MC}$ is strong due to the two cells pattern, yielding a strong radial transport in the equatorial regions. Despite some viscous transport associated with the radial shear, the angular momentum fluxes are mainly the result of a balance between the inward flux due to Reynolds stresses and the outward flux due to meridional circulation. The Reynolds stresses thus extract angular momentum from the surface and speed up the inner region whereas the meridional circulation flow, mostly directed outward at the equator, does the opposite job and extracts momentum from the deep inner parts.

Turning to the latitudinal fluxes, angular momentum is transported equatorward by Reynolds stresses, that act against the poleward transport by meridional circulation. This is consistent with the poleward streamlines existing in each hemisphere that corresponds to a counterclockwise circulation in the northern hemisphere. A small equatorward viscous transport occurs in both hemispheres due to the direction of latitudinal variation of Ω.

For case RG2, the radial transport of angular momentum by meridional circulation is weak. Indeed, for this slowly rotating simulation, the meridional circulation becomes unicellular, with one big cell covering the whole convective shell, and no radial transport is achieved in the equatorial regions, which are the regions that contribute the most in the integration of the radial fluxes. On the other hand, the viscous transport becomes important, related to the strong radial shear existing in this model in absence of radial transport by meridional circulation. The mechanisms yielding important transport of angular momentum are thus viscous terms balancing the Reynolds stresses, a configuration which is radically different from that achieved in model RG1 due to the very peculiar form of the meridional circulation appearing in the slowly rotating models. Let us note however that the net radial flux of angular momentum achieved in case RG2 is similar in amplitude to that achieved in case RG1. The transport due to Reynolds stresses is weaker in case RG2 compared to case RG1 because in that simulation, the convective pattern is less dominated by north–south aligned banana cells, and develops a more isotropic flow in which correlations are weaker.

The latitudinal fluxes are less symmetric with respect to the equator in case RG2, due to the large single cell of meridional circulation. The temporal average encompasses periods during which a second cell appears in the southern hemisphere (e.g., Figure 14), which translates into some weak latitudinal transport of angular momentum in this hemisphere. Similar to what is obtained for case RG1, the latitudinal flux balance is mostly between Reynolds stresses and meridional circulation because shellular rotation leads to almost zero viscous transport.

In both directions and for both cases, we note that the different processes balance each other so that the simulations tend to reach a statistical equilibrium in which almost no net flux is found in the radial or latitudinal direction. Due to the nonlinearity present in our turbulent simulations, it is to be expected that the total net flux of angular momentum (represented by the solid lines) is not strictly zero everywhere. However, averaging over several rotation periods leads already to a fairly good balance, confirming the mature state of the analyzed solutions.

We now turn to discussing the transport of heat with latitude achieved in our RGB models. We find that there are large fluctuations of temperature and entropy in our simulations. In Figure 17, we present temporal and longitudinal averages of the temperature fluctuations for cases RG1 and RG2. In model RG1, we find an organized symmetric pattern with respect to the equator, consisting of a warm equator and cooler poles. On the other hand, the northern hemisphere in RG2 is mostly hot while the southern hemisphere is cold, leading to an anti-symmetric temperature fluctuations profile consistent with the presence in this model of a dominant ℓ = 1 mode. Again there is a clear dichotomy between the two series of models, the cases rotating faster having smaller fluctuations and a more structured (banded) thermodynamic background.

Zoom In Zoom Out Reset image size

Convective motions under the influence of rotation lead to efficient transport of heat in latitude or what is sometime referred to as anisotropic heat transport (Rüdiger et al. 2005). As with the transport in radial direction (see Figure 3), several physical processes contribute to the transport of heat in latitude. In order to determine which physical processes dominate in our models we have derived a new diagnosis in the ASH code allowing us to assess the latitudinal energy flux balance. By taking the scalar product between the velocity $\bf v$ and the Navier–Stokes equation (see Equation (2)) one can derive the equation for the temporal evolution of the kinetic energy (Heyvaerts 1998; Miesch 2005). By adding the heat equation (see Equation (3)) to the kinetic energy equation and after some algebra consisting in putting under a flux conservative form some of the terms (Heyvaerts 1998), one gets the anelastic equation for the total density energy (Miesch 2005):

$$\begin{equation} \frac{\partial (U + {\rm KE})}{\partial t} = - \mbox{\boldmath $\nabla $}\cdot {\bf F} + {\cal B,} \end{equation} \tag{ 10 }$$

where $U = \bar{\rho }\bar{T}S$ is the internal density energy, ${\rm KE} = \frac{1}{2} \bar{\rho }v^2$ is the kinetic density energy, ${\cal B} = -v_r\bar{\rho }\bar{T}d\bar{S}/dr$ is the work against the stratification, and F is the total energy flux. The ${\cal B}$ term comes about because of anelastic formulation (see Bannon (1996) for a modified anelastic formulation in which the equation of state is altered to make that term disappear). F can be decomposed into several contributions:

$$\begin{equation} {\bf F} = - {\bf v} \cdot {\cal D} - \kappa \bar{\rho }\bar{T}\mbox{\boldmath $\nabla $}(S + \bar{S}) - \kappa \bar{\rho }C_p \mbox{\boldmath $\nabla $}(T + \bar{T}) + \frac{1}{2}\bar{\rho }v^2 {\bf v} + \bar{\rho }C_p {\bf v} T, \end{equation} \tag{ 11 }$$

the various terms being in turn heat transport by viscous stresses, by entropy eddy diffusion, by radiative diffusion, along with the kinetic energy and enthalpy fluxes. By assuming stationarity and that there is no net mass flux through horizontal surfaces, i.e., $\langle \bar{\rho }v_r \rangle _{\theta \phi } = 0$ (with 〈 〉_θϕ denoting an horizontal average), such that ${\cal B}$ vanishes, one is left with the simple flux conservation equation $\mbox{\boldmath $\nabla $}\cdot {\bf F} = 0$.

In the latitudinal direction, F_θ has the following expression (Elliott et al. 2000; Rüdiger et al. 2005; Miesch 2005):

$$\begin{equation} F_{\theta } = - ({\bf v} \cdot {\cal D})_{\theta } - \frac{\kappa \bar{\rho }\bar{T}}{r} \frac{\partial S}{\partial \theta } - \frac{\kappa \bar{\rho }C_p}{r} \frac{\partial T}{\partial \theta } + v_{\theta }\frac{1}{2}\bar{\rho }v^2 + \bar{\rho }C_p v_{\theta } T. \end{equation} \tag{ 12 }$$

One may further decompose the kinetic and enthalpy fluxes into their mean and fluctuating components in order to isolate the contribution coming from the mean meridional circulation. Assuming that v can be decomposed into an axisymmetric (mean) part denoted 〈v〉 and a fluctuating contribution denoted v', it can easily be shown that:

$$\begin{equation} F_{\theta,\, e} = \bar{\rho }C_p \left[\underbrace{\langle v_{\theta }^{\prime } T \rangle }_{\rm turbulent} + \underbrace{\langle v_{\theta }\rangle \langle T \rangle }_{\rm mean} \right], \end{equation} \tag{ 13 }$$

where we recall that the temperature fluctuations T with respect to the background reference temperature $\bar{T}$ possesses a non-vanishing axisymmetric mean (e.g., 〈T〉 ≠ 0) as already discussed in Figure 17. For the kinetic energy flux, we similarly get:

$$\begin{equation} F_{\theta,\, k} = 0.5 \bar{\rho }\left[\underbrace{\langle v_{\theta }^{\prime } v^{\prime 2} \rangle + 2 \langle v \rangle \langle v_{\theta }^{\prime } v^{\prime } \rangle + \langle v_{\theta }\rangle \langle v^{\prime 2} \rangle }_{\rm turbulent} + \underbrace{\langle v_{\theta }\rangle \langle v \rangle ^2}_{\rm mean} \right]. \end{equation} \tag{ 14 }$$

We will assume that the last term of the two above expressions corresponds to the latitudinal transport of mean enthalpy and kinetic energy by the mean meridional (latitudinal) flow. We then integrate over the radius r for each θ (as we did in Section 5.2), to get the net flux through cones of various colatitudes and deduce the latitudinal heat flux balance. Figure 18 displays the balance obtained for cases RG1 and RG2. For the faster model RG1, we see that the turbulent F_e (we omit the θ index for simplicity) is transporting heat poleward while both the turbulent F_k and its mean part (cf. Equation (14)) are warming up the equator balancing the turbulent F_e. The mean enthalpy flux F_e,m=0 (cf. Equation (13)) is constant with poleward heat transport, helping its turbulent counterpart. All the other terms contribute little to this latitudinal transport. For case RG2, the overall picture is messier and significantly different. The mean enthalpy flux plays a more important role certainly because of the large one cell meridional circulation established in that model, while strikingly enough the fluctuating one is almost negligible (i.e., has almost no net latitudinal effect). The turbulent F_k is now poleward while its mean contribution remains equatorward thus in the same direction as the mean enthalpy flux. The balance achieved in this model is not as good as in RG1 mostly due to the fact that the model possesses larger temporal fluctuations and that several terms of comparable amplitude contribute to the latitudinal heat transport. Similar results for the more turbulent cases RG1t and RG2t are found. In summary, we note that in all models a subtle latitudinal heat balance between the various processes is established, involving mostly the enthalpy and kinetic fluxes and that their relative importance including their mean and turbulent contributions is clearly influenced by the rotation rate.

Zoom In Zoom Out Reset image size

In rotating convection, both radial and latitudinal heat transport occur, as we have just discussed, establishing latitudinal gradients in temperature and entropy within the convective zone (see Figures 17 and 18). A direct consequence of these gradients is that the surfaces of pressure and density fluctuations will not coincide anymore, thereby yielding baroclinic effects. We can turn to the vorticity equations (Pedlosky 1987; Zahn 1992) to analyze the role of the turbulence and baroclinic effects in setting the large-scale flows discussed in Section 4. Briefly summarized, the baroclinic terms drive meridional flows that under the influence of Coriolis force yield longitudinal flows also known as the thermal wind. These baroclinic terms may also lead to a departure from cylindricity of the angular velocity profiles by braking Taylor–Proudman constraint (Pedlosky 1987).

The thermal wind balance equation can be derived from the vorticity equation as discussed in detail by Brun & Toomre (2002) and Miesch et al. (2006). The equation for the vorticity ω = ∇ × v in the purely hydrodynamical case can be derived under the anelastic approximation by taking the curl of the momentum equation (Equation (2); see also Fearn 1998 and Brun 2005 for its magnetohydrodynamic (MHD) generalization and the notion of magnetic (thermal) wind), it reads

$$\begin{eqnarray} \frac{\partial \mbox{\boldmath $\omega $}}{\partial t}&= &\!\!(\mbox{\boldmath $\omega $}_a\cdot \mbox{\boldmath $\nabla $}){\bf v} - ({\bf v}\cdot \mbox{\boldmath $\nabla $})\mbox{\boldmath $\omega $}_a - \mbox{\boldmath $\omega $}_a(\mbox{\boldmath $\nabla $}\cdot {\bf v})\\ \nonumber &\quad + &\frac{1}{\bar{\rho }^2}\mbox{\boldmath $\nabla $}\bar{\rho }\times \mbox{\boldmath $\nabla $}P - \mbox{\boldmath $\nabla \times $}\left(\frac{\rho g}{\bar{\rho }}\mbox{\boldmath $\hat{\bf e}_r$}\right) - \mbox{\boldmath $\nabla \times $}(\frac{1}{\bar{\rho }}\mbox{\boldmath $\nabla $}\cdot \mbox{\boldmath $\cal D$}), \end{eqnarray} \tag{ 15 }$$

with $\mbox{\boldmath $\omega $}_a=\mbox{\boldmath $\nabla \times $}{\bf v} + 2\mbox{\boldmath $\Omega _0$}$ the absolute vorticity, $\mbox{\boldmath $\omega $}=\mbox{\boldmath $\nabla \times $}{\bf v}$ the vorticity in the rotating frame, and $\cal D$ the viscous tensor (Equation (4)).

The ϕ component of the vorticity equation helps understanding the relative importance of the different processes acting in the meridional planes (Brun & Toomre 2002). In the stationary case ($\frac{\partial \mbox{\boldmath $\omega $}}{\partial t} = 0$), and assuming an axisymmetric configuration ($\frac{\partial }{\partial \varphi } = 0$) the azimuthal component of Equation (15) is

$$\begin{eqnarray} 2\Omega _0\frac{\partial v_{\phi }}{\partial z}&=&\underbrace{-(\mbox{\boldmath $\omega $}\cdot \mbox{\boldmath $\nabla $})v_{\phi } - \frac{\omega _{\phi }v_r}{r} - \frac{\omega _{\phi }v_{\theta }\cot \theta }{r}}_{\rm Stretching}\nonumber \\ &&+\, \underbrace{({\bf v}\cdot \mbox{\boldmath $\nabla $})\omega _{\phi } +\frac{v_{\phi }\omega _r}{r} + \frac{v_{\phi }\omega _{\theta }\cot \theta }{r}}_{\rm Advection} \nonumber\\ &&-\, \underbrace{\omega _{\phi }v_r\frac{d\ln \bar{\rho }}{d r}}_{\rm Compressibility} +\underbrace{\frac{g}{r c_p}\frac{\partial S}{\partial \theta }}_{\rm Baroclinicity} + \underbrace{\frac{1}{r \bar{\rho }c_p}\frac{d\bar{S}}{dr}\frac{\partial P}{\partial \theta }}_{\rm Departure\, from\, adiabaticity}\\ &&+\,\underbrace{\frac{1}{r}\left[\frac{\partial }{\partial r}(r {\cal A}_{\theta }) - \frac{\partial }{\partial \theta }{\cal A}_r \right],}_{\rm Viscous\; stresses} \nonumber \end{eqnarray} \tag{ 16 }$$

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

$$\begin{eqnarray} {\cal A}_r&=& \frac{1}{\bar{\rho }}\left[\frac{1}{r^2}\frac{\partial (r^2{\cal D}_{rr})}{\partial r}+\frac{1}{r\sin \theta }\frac{\partial (\sin \theta {\cal D}_{\theta r})}{\partial \theta } - \frac{ {\cal D}_{\theta \theta } + {\cal D}_{\phi \phi }}{r} \right] \nonumber \\ {\cal A}_{\theta }&=& \frac{1}{\bar{\rho }}\left[\frac{1}{r^2}\frac{\partial (r^2{\cal D}_{r\theta })}{\partial r}+\frac{1}{r\sin \theta }\frac{\partial (\sin \theta {\cal D}_{\theta \theta })}{\partial \theta } + \frac{ {\cal D}_{\theta r} - cot \theta {\cal D}_{\phi \phi }}{r} \right] \nonumber \end{eqnarray}$$

In the above expressions, we have omitted for the sake of clarity the 〈 〉 denoting the axisymmetric average.

Here, we have identified several contributions to characterize the shape of the angular velocity along the z-axis that we display (only for the dominant term) in Figure 19 for case RG1 and in Figure 20 for case RG2: Stretching describes the stretching/tilting of the vorticity due to velocity gradients; Advection describes the advection of vorticity by the flow; compressibility describes the stretching of vorticity due to the flow Compressibility; ${\frac{g}{r c_p}\frac{\partial S}{\partial \theta }}$ is the baroclinic term; $\frac{1}{r \bar{\rho }c_p}\frac{d\bar{S}}{dr}\frac{\partial P}{\partial \theta }$ results from departure from adiabaticity; Viscous accounts for the diffusion of vorticity due to viscous stresses. Figure 19 displays for RG1 the left-hand side of Equation (16), along with the dominant terms of the right-hand side, their sum and the difference RHS − LHS. The first point to notice is that the thermal wind balance, which would be achieved in the stationary case by having a strict equality between the LHS and the baroclinic term, is not achieved in our simulation (contrary to what we find in most of the solar convective envelope Brun & Toomre (2002)). We note that LHS ∼ RHS which is a hint of the mature dynamical state of the model. The dominant terms for the evolution of vorticity also involve both turbulent transport by advection thus emphasizing the key role played by turbulent convection in setting the angular velocity profiles as discussed in Section 5.2. The viscous stresses contribute in key regions (mostly at the top of the domain) where the isocontours of Ω are distorted. Overall all these large terms cancel each other thus resulting on isocontours of Ω in RG1 that are close to cylinders as observed. A similar analysis of RG2 reveals again that there is no strict thermal wind balance. Large advection terms oppose the baroclinic contribution with compressibility also playing a role. The quality of the achieved balance is not as good in part due to the fact that the meridional flows possess large fluctuations that make temporal average less meaningful, with the direct consequence that the temporal derivative of ω_ϕ cannot be neglected as easily. The shellular shape of the differential rotation in the slowly rotating cases RG2 and RG2t is clearly neither cylindrical nor conical, and it is unclear given their large convective Rossby number (cf. Table 1) and the weak influence of rotation on their dynamics whether their angular velocity profile can be usefully analyzed with the thermal wind balance contrary to RG1 and RG1t, for which the diagnosis is more adequate. We thus conclude that even though the amplitude of the baroclinic effect is large, it does not play a role in shaping the angular velocity profile as important as in solar simulations, at least at the rotation rates that we have explored in this study.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size