{'fraction_non_alphanumeric': 0.03202862132118063, 'fraction_numerical': 0.01482175561139742, 'mean_word_length': 4.636658264314008, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Differential rotation in solar-like stars from global simulations

To explore the physics of large-scale flows in solar-like stars, we perform 3D anelastic simulations of rotating convection for global models with stratification resembling the solar interior. The numerical method is based on an implicit large-eddy simulation approach designed to capture effects from non-resolved small scales. We obtain two regimes of differential rotation, with equatorial zonal flows accelerated either in the direction of rotation (solar-like) or in the opposite direction (anti-solar). While the models with the solar-like differential rotation tend to produce multiple cells of meridional circulation, the models with anti-solar differential rotation result in only one or two meridional cells. Our simulations indicate that the rotation and large-scale flow patterns critically depend on the ratio between buoyancy and Coriolis forces. By including a subadiabatic layer at the bottom of the domain, corresponding to the stratification of a radiative zone, we reproduce a layer of strong radial shear similar to the solar tachocline. Similarly, enhanced superadiabaticity at the top results in a near-surface shear layer located mainly at lower latitudes. The models reveal a latitudinal entropy gradient localized at the base of the convection zone and in the stable region, which however does not propagate across the convection zone. In consequence, baroclinicity effects remain small and the rotation iso-contours align in cylinders along the rotation axis. Our results confirm the alignment of large convective cells along the rotation axis in the deep convection zone, and suggest that such"banana-cell"pattern can be hidden beneath the supergranulation layer.


Introduction
Stellar magnetic activity is thought to be a result of inductive effects of large-scale shearing flows as well as complicated collective effects of small-scale turbulent motions in convective envelopes. Cyclic behaviour has been detected in several late-type stars in the emission cores of the Calcium H and K line profiles, with periods between 7 and 22 years (Baliunas & Vaughan 1985;Hall 2008). The best known and studied example is, however, the 11-years sunspot cycle, where observed activity is caused by yet poorly understood dynamo processes in the convection zone occupying the upper 30% of the solar radius. A similar mechanism is expected to operate in other solar-like stars.
Several different dynamo models have been proposed over the years in order to describe the solar cycle using a MHD mean-field (Steenbeck et al. 1966) formalism (see Brandenburg & Subramanian 2005;Charbonneau 2010, for reviews on the topic). However, since the dynamics of the magnetic field in the solar interior remains hidden to the observations, there are strong model ambiguities.
The large-scale plasma flows in the solar interior, such as differential rotation and meridional circulation have been probed thanks to helioseismology. We know that the Sun rotates differentially in latitude through the whole convection zone with a profile that forms conical iso-rotation contours (Schou et al. 1998). Below the bottom of the convection zone the rotation becomes rigid creating a thin radial shear layer (Kosovichev 1996) called the tachocline (Spiegel & Zahn 1992). Besides, at the upper 50Mm of the convection zone the rotation rate decreases apparently at all latitudes forming another thin layer of negative radial shear so-called near-surface shear layer (Thompson et al. 1996;Corbard & Thompson 2002).
The meridional circulation has been observed at the solar surface through different techniques (e.g., Giles et al. 1997;Zhao & Kosovichev 2004;Ulrich 2010;Hathaway & Rightmire 2010;Gizon & Rempel 2008). This flow is directed poleward with amplitudes ∼ 20 m/s peaking at mid latitudes. By tracking supergranules Hathaway (2012) was recently able to infer a return flow at 70 Mm below the photosphere. Besides, helioseismic inversions by Zhao et al. (2012Zhao et al. ( , 2013 revealed that the return flow is located in the middle of the convection zone and suggested the existence of a second circulation cell located deeper.
More recently asteroseismology opened the door to the exciting possibility of measuring the large scale flows in different classes of stars. Of special interest are stars at the same stage of evolution and structurally similar to the Sun. The study of these stars could provide a better understanding of the stellar dynamo mechanism and the generation of cosmic magnetic fields in general.
The observational results from helioseismology have imposed some constrains to the mean-field dynamo models but are still insufficient to provide a complete scenario.
Thus, global numerical simulations seem to be the most promising approach towards the understanding of the solar dynamics and dynamo, filling the gap left by observations. Steady dynamo solutions have been found in numerical models at the solar rotation rate.
These models showed that the magnetic field in the convection zone is organized in the form of toroidal wreaths (Brown et al. 2010). Oscillatory dynamo solutions have also been found in simulations of stars rotating at 3 or more times faster than the Sun (Käpylä et al. 2012;Nelson et al. 2013). However, none of the current numerical models has been able to reproduce sufficiently well the solar differential rotation and meridional circulation, in particular, the tachocline and the subsurface shear layer, and still are not able to explain the solar activity cycles.
In spite of the fast developments of computer capabilities, all global simulations are far from realistically modeling the Sun or solar-like stars. The problem arises from the low values of the dissipative coefficients in the stellar interior which impose a large separation among the scales of energy injection and diffusion. Most of the simulations fall in the group of the so-called Large-Eddy Simulation (LES) models. These models resolve explicitly the evolution of large-scale motions while the contribution of unresolved, sub-grid scale (SGS), motions is represented by ad-hoc turbulence models (e.g. Smagorinsky 1963). For many years the small-scale contribution used in solar global simulations was the most simple, e.g., in the form of enhanced eddy viscosity (e.g., Gilman 1976;Brun et al. 2004). Recently, more sophisticated SGS models, like the dynamic Smagorinsky viscosity model (Germano et al. 1991), have been implemented leading to a more turbulent regime in the global models (Nelson et al. 2013).
An interesting alternative to such explicit SGS modeling consists in developing numerical schemes which accurately model the inviscid evolution of the large scales whilst the SGS contribution is intrinsically captured by numerical dissipation. This class of modeling is called Implicit Large Eddy Simulations (ILES), and is implemented in the EULAG code (Smolarkiewicz 2006;Prusa et al. 2008), used primarily in atmospheric and climate research.
Recently, oscillatory dynamo solutions were studied in global simulations performed with the EULAG code for the case of the solar rotation (Ghizaru et al. 2010). Although the period of the simulated activity cycles is of ∼ 40 years (instead of 11), and also there is no clear latitudinal migration of the toroidal field towards the equator, these results are a promising step forward in the understanding of the solar dynamo. These simulations with EULAG were mostly focused on the dynamo rather than on the differential rotation problem. The present work aims to explore the hydrodynamical case of rotating convection with the EULAG code, giving attention to the physics behind the development of large-scale flows: meridional circulation and differential rotation.
The main challenge is to reproduce the solar differential rotation profile inferred by helioseismology with the characteristics described above. Starting with the seminal work by Gilman (1976) this problem has been studied numerically during the last few decades without obtaining a completely satisfactory result. This is a complex problem, and its solution must fulfill several observed properties: i reproduce the latitudinal differential rotation with a difference in the rotation rate between the equator and 60 • latitude of ∼ 90nHz; ii obtain iso-rotation contours of Ω(r, θ) aligned along conical lines in the bulk of the convection zone; iii reproduce the sharp transition region between the differentially rotating convection zone and the solid body rotating radiative zone, i.e., the tachocline, iv reproduce a sharp decrease in the solar rotation rate in the upper 50Mm of the solar convection zone, in other words, a near-surface rotational shear layer with a negative radial gradient of the angular velocity.
The item i) has been explored since Gilman (1976), and is understood as a combined effect of the buoyancy and Coriolis forces. Rotation leads to correlations between the three components of convective motions, reflected in anisotropy of the Reynolds stress tensor, which determines the direction of the angular momentum flux. In global simulations with the anelastic code ASH, Brun & Toomre (2002) scanned the parameter space, particularly in terms of the Rayleigh (Ra) and Prandtl (Pr) numbers, to determine a trend at which the latitudinal gradient of differential rotation could be sustained at higher levels of turbulence.
Item ii) is perhaps the most cumbersome of all the items above. Most of the global simulations of rotating convection obtain the iso-lines of angular velocity aligned along the rotation axis, following the so-called Taylor-Prouman state, which can be explained from the φ component of the vorticity equation: where ω φ is the azimuthal component of the vorticity, Ω is the angular velocity, g is the gravity acceleration, s is the entropy and c p is the specific heat at constant pressure. The derivative along the cylindrical radius is defined as ∂ z = cos θ∂ r − r −1 sin θ∂ θ .
In the steady state (∂ω φ /∂t = 0), the cylindrical iso-rotation could be attributed to insufficient baroclinicity in the system, i.e., if entropy does not exhibit a pronounced latitudinal gradient then ∂Ω 2 /∂z = 0 . In hydrodynamic mean-field models, the Taylor-Proudman state is broken by considering anisotropic heat conductivity which could result from rotational influence on the turbulent eddies (Kitchatinov & Ruediger 1995). Similarly, Rempel (2005) studied how a subadiabatic stratification beneath the convection zone can produce a gradient in the entropy perturbations that could be transported into the convection zone. Miesch et al. (2006) implemented the Rempel's idea in global convection simulations by imposing a latitudinal gradient of entropy as a bottom boundary condition.
They found with such model that the iso-contours of angular velocity form conical lines such as determined by helioseismology.
As for the item iii), the tachocline has been included in dynamo simulations by imposing a forcing term in the momentum equation (Browning et al. 2006 (Ghizaru et al. 2010;Racine et al. 2011;Guerrero et al. 2013, and the ones presented in this paper) include a strongly sub-adiabatic radiative zone, and have been able to reproduce the rotational shear layer in agreement with the solar observations.
The Furthermore, recent studies of the amplitude of convective velocities at the solar interior have pointed out the discrepancy between the turbulent velocities obtained by helioseismology observations (Hanasoge et al. 2012), theoretical estimations (Miesch et al. 2012) and numerical models (e.g., Miesch et al. 2008). This has emphasised the need of a better understanding of the multiscale turbulent velocities in the solar convection zone.
Both issues, the meridional circulation and turbulent velocities deserve further investigation.
Having a solar model as our reference, in this paper we investigate the items (i) to (iv) for the Sun and also study the generation of differential rotation and meridional flows in solar-like stars. Specifically, we study turbulent convection in spherical shells, the stratification of which resembles the structure of the solar interior. Considering simulations for different rotation rates we first explore the phenomena of angular momentum transport and its relationship to different configurations of mean flows within the convective layer. As a test for the implicit SGS method employed in this paper, we next explore the convergence of the numerical scheme with increasing resolution. Finally, considering angular momentum transport arguments we explore the possibility of forming a near-surface rotational shear layer such as observed in the upper layers of the solar convection zone.
This paper is organized as follows: in §2 we describe the numerical model, in §3 we present our results, including the properties of convective flows and the resulting differential rotation and meridional circulation. We also describe the angular momentum mixing in terms of the Reynolds stresses. Finally we discuss our attempt to reproduce the near-surface shear layer. We conclude and discuss our results in §4.

Model
Similarly to the previous dynamo simulations with the EULAG code (Ghizaru et al. 2010;Racine et al. 2011), our model covers a full spherical shell, i.e., 0 ≤ θ ≤ π, and 0 ≤ φ ≤ 2π in latitude and longitude, respectively. In radius our domain spans from 0.61R ⊙ to 0.96R ⊙ (although in section 3.3 we extend our domain up to r = 0.985R ⊙ ).
We solve the anelastic set of hydrodynamic equations following the formulation of where D/Dt = ∂/∂t + u · ∇ is the total time derivative, u, is the velocity field in a rotating frame with constant angular velocity Ω, p ′ and Θ ′ are the pressure and potential temperature fluctuations with respect to an ambient state (discussed below), respectively, H(θ ′ ) represents the radiative and heat diffusion terms, ρ s and Θ s are the density and potential temperature of the reference state chosen to be isentropic (i.e., Θ s = const) and in hydrostatic equilibrium, g = GM/r 2 is the gravity acceleration. The potential temperature, Θ, is related to the specific entropy through s = c p ln Θ + const.
Finally, the term αΘ ′ represents the balancing action of the turbulent heat flux responsible for maintaining the steady, axi-symmetric, solution of the stellar structure (see next section). It is characterised by the ambient state Θ e (together with the corresponding pressure p e ). In the reported simulations, the αΘ ′ term forces the system towards the ambient state on a time scale τ = α −1 = 1.55 × 10 8 s (∼ 5 years). This timescale is much shorter than the timescale of radiative losses and heat diffusion, meaning that the turbulent heat flux is the primary driver of convection in the system (see also Cossette et al. 2013).
After verifying that the terms in H(θ ′ ) do not have any effect on the model solutions, we did not included these terms in the simulations below. However, the choice of τ can affect the results since the final potential temperature is the total of the ambient state and fluctuations. The longer (shorter) τ , the larger (smaller) the value of Θ ′ . We have verified (in simulations not presented here) that models with shorter τ suppress convection and reproduce a profile close to Θ e . In models with longer τ , on the other hand, the large fluctuations tend to homogenize the final potential temperature creating a flat, more adiabatic, profile. Consequently, the selected value of τ has been determined empirically as a compromise between a weak forcing of the energy equation and a reasonable relaxation time to a statistically steady state of the system, that still leads to the desired results while maintaining convection and the stellar structure over a long observation time.

Ambient state
The ambient state in the simulations is approximated by a polytropic model. The polytropic stratification for an ideal gas, p = ρRT , is given by the hydrostatic equation: where g(r) = g b (r b /r) 2 is the gravity acceleration, r b is the radius of the bottom boundary, g b is the gravity acceleration at r b , R is the gas constant, and m is the polytropic index.
The solution of Eq. (5) integrated from r b to r and assuming m = const is: are the values of the temperature, density and pressure at the bottom of the domain, r b . These equations are solved recursively from bottom to top for a polytropic index that varies with radius as follows: with ∆m = m cz − m r , m r = 2.5 and m cz = 1.49995. The potential temperature -12associated with this stratification is computed through: where γ = 5/3 defines the adiabatic exponent of an ideal gas. This setup defines a strongly sub-adiabatic radiative layer (stable to convection) below r tac = 0.7R ⊙ , followed by a slightly super-adiabatic, convectively unstable layer. An erf function couples the two shells with a width of transition w t = 0.01R ⊙ (see Fig. 1). shows a zoom of the convection zone from r = 0.7R ⊙ to r = 0.96R ⊙ .

Numerical method and boundary conditions
Eqs.
(3) to (4) are solved numerically using EULAG code (Smolarkiewicz et al. 2001;Prusa et al. 2008). EULAG is an HD and MHD code for incompressible or anelastic fluid dynamics developed originally for research of atmospheric processes and climate and generalized over years for a range of geophysical and astrophysical problems. The time evolution is calculated using a unique semi-implicit method based on a high-resolution  (Elliott & Smolarkiewicz 2002;Domaradzki et al. 2003;Margolin et al. 2006). Thus, the results of MPDATA are often interpreted as Implicit LES (ILES) models (Smolarkiewicz & Margolin 2007).
The boundary conditions are specified as follows. In the latitudinal direction, discrete differentiation extends across the poles, while flipping sign of the longitudinal and meridional components of differentiated vector fields. In the radial direction we use stress free, impermeable, boundary conditions for the velocity field, and asume that the normal derivative of Θ ′ is zero. The simulations are initialized applying a small-amplitude random perturbation to the velocity field; the results presented in the next section are computed once the simulations have reached a statistically steady state. It is characterized for the random fluctuations of the volume averaged fluid quantities around a well defined mean.

Results
The convective properties of the non-rotating solution of Eqs. (2-4) exhibit convection cells roughly of the same spatial extent with broad and slow upflows and sharp and fast downflows (see Fig. 2a). The horizontal scale of this motion depends on the resolution of the model. Initially we consider the same number of grid points as in Ghizaru et al. (2010), n φ = 128, n θ = 64, n r = 47. Since the viscous terms are not explicitly considered in the equations, the standard Reynolds or Taylor numbers, cannot be computed. However, following Christensen & Aubert (2006), we consider a modified Rayleigh number which does not depend on the dissipative terms, Ra * = 1 where s e is the specific entropy of the ambient state. From the solution we compute the Mach number, Ma = u rms /c * s , where u rms is the volume averaged RMS turbulent velocity and c * s = γRT * s , T * s is the temperature of the background state at the middle of the unstable layer (r = 0.85R ⊙ ); and the Rossby number, which compares the relative importance of turbulent convection over rotation: where D is a typical length scale, given here by the thickness of the convective unstable layer, 0.25R ⊙ . In addition we compute a parameter characterizing the latitudinal differential rotation; where Ω eq = Ω(R ⊙ , 90 • colatitude) and Ω p = Ω(R ⊙ , 30 • colatitude). Overline represents azimuthal averaging. In the Sun χ Ω ≃ 0.2. For models with different Rossby numbers different values of χ Ω are expected due to different velocity correlations. In the next section we study these correlations and the development of large scale motions in terms of the balance of angular momentum fluxes. g ∂se ∂r , u rms is the volume average rms velocity (in m/s) in the unstable layer, Ma = u rms /c * s is the Mach number, with c * s = γRT * s | r=0.85R ⊙ being the sound speed at the middle of the unstable layer, Ro = urms   Table 1). Yellow/dark color scale represent upward/downward motions.

Forces balance and angular momentum fluxes
The development and sustenance of large-scale flows depend on correlations between the components of the turbulent velocities, the so-called Reynolds stresses. The relative importance of these quantities defines the redistribution of angular momentum. Among these correlations, of particular importance are the vertical velocity and the rotation rate of the system since they define the formation of fast rotating zonal flows either at higher or lower latitudes. To understand this behaviour in our first set of numerical experiments we change the rotation rate while keeping the ambient state constant (i.e., the same gradient of potential temperature, Eq. 8).
Models L1 to L7 in Table 1 summarize the parameters and results for this set of simulations. The structure of the convective pattern for some representative cases is shown in Fig. 2b- where the over-line corresponds to the time and longitude average, u m and u ′ m are the mean and turbulent meridional (r and θ) components of the velocity, respectively. In steady state the left hand side of the equation vanishes remaining The first term inside the divergence corresponds to the angular momentum flux due to meridional circulation and the second one is the flux due to small-scale correlations, the Reynolds stresses. Each of these terms has radial and latitudinal component, A more illustrative representation of these quantities is in the cylindrical coordinates ̟ and -18z thus: All these quantities were computed in the statistically steady state attained after a long time interval depending on the rotation rate, and averaged in time for 3.4 years.
In the case of the slower rotation rate, Ω 0 = Ω ⊙ /4 (run L1), the convective pattern exhibits big convection cells which are slightly elongated in the direction opposite to the rotation. They span the range ∼ ±45 • of latitude. As for the mean flows, the azimuthal motion is slower at the equator than at the higher latitudes (left panel of Fig. 3a). The transport of angular momentum due to Reynold stress in the direction of ̟ is negative at all latitudes (third panel of Fig. 3a), indicating that the rotation increases inwards. The most important contribution to the angular momentum flux is, however, the meridional circulation (rightmost panel of Fig. 3a). Single large scale meridional flow cells develops at each hemisphere. The circulation is counterclockwise (clockwise) in the northern (southern) hemisphere (second panel of Fig. 3a). This flow transports the angular momentum from the lower to higher latitudes so that the zonal flow is accelerated in the direction of rotation at latitudes between 30 • and 80 • .
Increasing the rotation rate to Ω 0 = Ω ⊙ /2 (Run L2) or Ω 0 = 0.57Ω ⊙ (Run L3) we find that the convection cells are concentrated in a band around ∼ ±30 • latitude. In this region the cells are mostly elongated in the longitudinal direction (Fig. 2b). Similarly to the previous model, the latitudinal differential rotation is anti-solar (Fig. 3b). However, the radial negative shear at the equator starts not at the base of the convection zone, like in case L1, but at r ≃ 0.73R ⊙ . Between 0.7R ⊙ < r < 0.73R ⊙ the plasma rotates faster than the rotating frame. In this case the largest contribution to the mixing of angular momentum comes again from the meridional circulation. However, this time two strong meridional flow cells in radius are formed in each hemisphere. In Run L2 the latitudinal velocities are so high (∼ 40 m/s) that the flow crosses from one hemisphere to another (second panel of Fig.   3b). Thus, the meridional flow is not fully symmetric across the equator even after a long temporal average.
In runs L4 and L5 the contributions of meridional circulation and Reynold stresses to the flux of angular momentum are similar. For these cases the plasma at lower latitudes rotates roughly at the same rate as the rotating frame. Above ∼ 60 • latitude the zonal flow decelerates and has the minimal rotation rate at the poles. The meridional flow shows a multicellular pattern with three or more cells in radius per hemisphere, with latitudinal velocities of ∼ 10 m/s. Run L6 (Ω 0 = Ω ⊙ ) shows a rotation profile accelerated at the equator (Fig. 3c).
Similar to the Sun, the rotation decreases monotonically towards the poles, and is in iso-rotation with the stable layer of the domain at middle latitudes. In this case the Reynolds stress components of the angular momentum flux dominate over the meridional circulation components. Although at lower latitudes, at the bottom of the convection zone F RS ̟ is negative, it changes sign above r ∼ 0.85R ⊙ . However, F RS z , which transports the angular momentum towards the equator is more important. Something similar happens in the fastest rotating case, L7 (Ω 0 = 2Ω ⊙ , Fig. 3d) where the Reynolds stress drive the differential rotation. In this case, however, the transport along the rotation axis, F RS ̟ , is the most important flux component. The rotation is faster at the equator and slower at the poles, but it exhibits jets of slow rotation in a mid latitude range which corresponds to a cylinder tangent to the base of the convection zone. The term F M C z has larger values but varies quickly in radius and latitude. Thus, there are no organized meridional circulation flows capable to effectively advect the angular momentum. Thus, in general, the differential rotation establishes as a competition between the buoyant and Coriolis forces. When buoyancy dominates then the differential rotation is of the anti-solar type, and when the Coriolis force rotation dominates the rotation is solar-like.
To demonstrate this point we have performed the simulations labeled as L8 and L9 in Table 1, for which the rotation has the solar value, but the profile of Θ e was modified to decrease and increase the buoyancy term by making the stratification less and more superadiabatic, respectively (see Fig. 4). The model with less superadiabatic stratification, L8, develops a faster rotating equator as a consequence of the dominant role of the Corioles force. The distribution of the angular momentum fluxes seems to be an intermediate case between the models L6 and L7 (Fig. 6a). On the other hand, a more superadiabatic Θ e results in antisolar rotation because the buoyancy force dominates in this case. The results of this model are compatible with the case L2 or L3 (Fig. 6b). Another way to modify the adiabaticity of the system is considering different relaxation times, τ (1/α in Eq. 4). We have verified in simulations not presented here that models with shorter τ develop small fluctuations of Θ ′ , so that the final profile of Θ is more superadiabatic. Also, models with longer τ develop large fluctuations that tend to homogenize the potential temperature in the bulk of the convection zone, leading to a less superadiabatic system. Changes in τ , however, are not as important as directly modifying the profile of Θ e . It is important to notice that qualitatively the shape of the iso-rotation contours does not depend on the dominance of the buoyancy or Coriolis forces; in all cases the rotation contours are aligned along the rotation axis. To study the role of baroclinicity in the resulting Taylor-Proudman balance, in Fig. 5 we plot the azimuthal average of the potential temperature fluctuation, Θ ′ (r, θ), for three representative models: L2, L6 and L7. Slow rotation models (L1-L3) exhibit small variations of potential temperature predominantly in the radial direction. In rapid rotating models (L4-L7), below the convection zone, Θ ′ aquires larger values, it is positive at the poles and negative at the equator, indicating warmer poles. However, within the convection zone, this latitudinal gradient does not propagates and Θ ′ varies mainly radially. The lack of a latitudinal gradient in temperature explains why the contours or differential rotation remain cylindrical (see Eq. 1). The reasons and implications of this behavior will be explored in a different paper.
The transition between the solar and anti-solar rotation regimes (faster or slower equator) has been obtained since the early models of Gilman (1976) where he compared the Rayleigh with the Taylor numbers. In his Fig. 6 he depicts the regimes for transition between solid rotation, equatorial acceleration and high latitude acceleration. Similar results for the differential rotation and meridional circulation in the context of rotation of gas planets were recently reported by Gastine et al. (2013). Performing simulations for different rotation rates they found a sharp transition from the prograde to retrograde rotation of equatorial regions. They compared the Rossby number at the equator with the Rayleigh number at a mid depth, Ra * , and found that this transition occurs at Ra * ∼ 1.
In this paper, we compare the differential rotation parameter, defined in Eq. (9) with the Rossby number as it is depicted in Fig. 7. For the models with the same ambient state we find that the sharp transition between the prograde and retrograde regimes occurs at Ro ≃ 0.063 (Ra * ≃ 3 in Table 1), indicated by the vertical dashed line in Fig. 7. In models with the same rotation rate but with higher (lower) superadiabaticity (models L8 and L9), the transition can happen at lower (higher) values of Ro. It is noteworthy to mention that the model L6 is in a good agreement with the solar value of χ Ω , indicated by the horizontal dashed-dotted line. For the model with faster rotation, L7, the differential rotation parameter, χ Ω , decreases. Here it happens due to the formation of a cylinder of slow rotation at intermediate latitudes. However, similar global simulations of fast rotating stars also show a decrease in the latitudinal differential rotation (Brown et al. 2008;Dubé & Charbonneau 2013). Käpylä et al. (2011b) also explored the transition from anti-solar to solar-like differential rotation in compressible simulation of convection in -23spherical wedges. In their Fig. 17 they compare the rotation parameter with the Coriolis number (inverse to Ro). Unlike our results and those in Gastine et al. (2013), in their case this transition seems to be smoother. However, due to all these differences in the model and experimental design this comparison remains inconclusive and warrants further studies.

Convergence test
As a test for the implicit SGS model, we performed a set of simulations presented for the same stratification (Θ e ) and the same rotation rate (Ω ⊙ ) but for different mesh resolutions. Run L6, presented in the previous section was for a relatively coarse grid model with 128 × 64 × 47 mesh points. Models M1 and L1 have 2 and 4 times higher resolution, respectively. Snapshots of the vertical velocity for the cases L6, M1 and H1 are shown, from top to bottom panels, in Fig. 8. In all three cases convective "banana" cells elongated along the rotation axis appear in a belt of ±30 • . However, the scale of these banana cells depends on the resolution, i.e., a larger number of azimuthal modes seems to be excited for higher resolutions. The higher latitudes are populated with smaller and more symmetrical convection cells whose horizontal scale also decreases with the increasing of the resolution.  (Table 1)  To study the energetics of the interacting large and small scales for each of these simulations we have computed the power spectrum of the kinetic energy at a given radius, using the Fourier transform in the azimuthal direction: where m is the azimuthal wave number, and the power spectrum definition (Mitra et al. 2009): where the brackets denote average over latitude. Note that E u = u 2 . The spectra for models L6, M1 and H1 are plotted in Fig. 9. For the coarser resolution (continuous black line) the spectrum shows an energy decay from the smaller wavenumber up to 7 m 15.
There are peaks of power at m ≃ 5 and m ≃ 7, which probably correspond to the scales of the banana cells. From m ≃ 15 to the smallest resolved scale, m = 64, there are two different sub-ranges, one of which is compatible with the inertial turbulent decay (although not exactly following the Kolmogorov k −5/3 law, see dot-dashed line) up to m ≃ 30, followed by a faster decay. In total, the energy spectra spans ∼ 3 orders of magnitude. The spectrum for the run M1 (red line) is similar to L6, but in this case the excess of power is at m ≃ 6.
There is a plateau of energy from 10 m 20 followed by a decay up to the smallest wave number m = 128. The model H1 with higher resolution (blue line) does not have a plateau, and the energy seems to decay continuously, with a seemingly inertial sub-range between 12 m 30. The peak at m ≃ 12 contains the energy of thin banana cells observed in this simulation. In total this model spans ∼ 5 orders of magnitude in energy. For comparison we have plotted also the power spectrum of the models without rotation for low, L0, and high, H0, resolutions (black and blue dotted lines). Evidently the convective energy is higher in the models without rotation. It is interesting how this simulations show a turbulent decay consistent with k −5/3 which spans several orders of magnitude in energy and for a long range of scales.  Table 1). Continuous (dotted) lines correspond to the models with (without, models L0 and H0) rotation. The dot-dashed line indicates the −5/3 power decay.
Similarly to Fig. 3, Fig. 10 shows the profiles of differential rotation and meridional circulation obtained for models M1 and L1, as well as the angular momentum flux components. In model M1 the mean rotation and meridional flows slightly resemble those of model L6. In this case, however, the rotation shows a band of rotation with the speed of the stable layer, which spans over a large latitudinal extent. Both quantities, the Rossby number, Ro, and the rotation parameter, χ Ω are similar to those of model L6, so that this case appears almost at the same position in Fig. 7. Regarding the angular momentum fluxes, the meridional distribution of F RS ̟ and F RS z is qualitatively similar to L6. In the case M1, the Reynolds stress contribution is about one half of that in case L6 (see Fig.   10a). On the other hand, F M C ̟ and F M C z increase by a factor of two. Given the resemblance in the Reynolds stresses profiles, the differential rotation does not differ abruptly between cases L6 and M1. The meridional motions are, however, more pronounced in this case.
Like in case L6, there are several cells in the meridional flow profile elongated in the z direction.
This contribution of meridional circulation to the angular momentum flux appears further enhanced in the high resolution model, H1 (Fig. 10b). Although there is an equatorial acceleration due to a positive F RS ̟ flux, both F RS ̟ and F RS z are decreased in amplitude. Most importantly, there is no morphological similarity to the cases L6 and M1.
This leads to a different profile of Ω. Instead of a continuous decrease with latitude, there is a retrograde jet which appears at the surface between 20 • and 40 • . This seems to be the consequence of meridional motions aligned along the tangent cylinder at r ≃ 0.71R ⊙ .
Another prograde jet appears at mid latitudes, from there up to the poles the azimuthal velocity equals the rotation rate of the frame, with a slight decrease towards the poles. In this case the meridional motions enhance F M C ̟ and F M C z by another factor of 2.
It is surprising that although the Rossby number, Ro, is similar for the models with three different resolutions, indicating that the balance between the Coriolis and buoyancy forces is roughly the same, in case H1 the value of χ Ω (≃ 0.09) is smaller than in the lower resolution models L6 and M1, while F M C is larger. It might be the case that when a more complex and turbulent flow develops, small-scale motions contributing to the angular momentum transport diffuse on shorter timescales so that the Coriolis force does not affect them. The appearance of these new scales of motion might also change the role of the sub-grid scale transport, relegating its action to scales not resolved. The situation could be that the smallest resolved scales are already modifying the system. The question that arises, and deserves further investigation, is why these effects are not captured by the sub-grid scale transport in the coarser resolution cases? It is puzzling why our model L6 with the rather coarse resolution reproduces the solar rotation closer than the higher resolution models, M1 and L1. Perhaps it may be necessary to include explicit eddy transport in the higher resolution models (i.e., include in the equation an explicit term for the turbulent heat diffusion or for the viscosity) in order to increase the efficiency of the unresolved turbulent transport. In other words, this will modify the effective Prandtl number of the simulation.
Another possibility could be to readjust the parametrized turbulent heat flux profile, Θ e .

The solar near-surface shear layer
In the upper part of the solar convection zone the velocities are such that the turnover time is of the order of minutes (for granulation) or hours (super-granulation). These time scales are smaller when compared with the solar rotation period, 28 days, indicating that there the buoyancy force is dominating over the Coriolis counterpart. Observations indicate that in this region the radial shear is negative, i.e., the angular velocity decreases with radius. In the simulations described in the previous sections, an inward angular momentum flux is clearly observed in the buoyant dominated regime. It corresponds to an angular velocity which decreases radially outwards. In order to simulate the buoyancy dominated regime, we can take advantage of our formulation of the energy equation and impose a rapid decrease of the potential temperature in the upper part of the convection zone. Although radiative cooling is not considered in the simulations, this sharp decline of the potential temperature (or entropy) occurs indeed at the upper part of the convection zone due to hydrogen ionization.
In the model labeled as N1 we extend the domain in the radial direction up to r = 0.985R ⊙ and modify the ambient state (Eq. 7) by considering the polytropic index m(r) = m r + ∆m 1 2 1 + erf r − r tac w t (16) where ∆m 2 = m nssl − m cz , with m nssl = 1.4996 and r nssl = 0.96R ⊙ (see black dotted line in Fig. 11).  Like in the Sun, the rotation profile shows acceleration at the equator. Above r = 0.96R ⊙ the rotation decreases with radius at lower latitudes (see blue line in Fig. 11).
This negative shear does not remain for latitudes above 15 • , and the contrast of differential rotation in this case is χ Ω = 0.087. Our interpretation is that the near-surface shear layer is associated with an inwards transfer of the angular momentum as a consequence

Discussion
Large scale flows in rotating turbulent convection, like in giant planets and stars, are developed due to collective turbulent effects. These flows are important in astrophysics since they could explain the observed patterns of surface differential rotation and, in addition, because they could be responsible for the generation of large-scale magnetic fields via the dynamo process.
In this paper we have explored the development of these large scale flows in global simulations of a convective envelope whose stratification resembles the solar interior.
Although our main goal is study the pattern of differential rotation observed in the Sun, our results with different rotation rates are applicable also to solar-like stars.
The model considered here has proven to work well in simulating stellar interiors. The inviscid numerical method allows higher turbulent levels even with a coarse resolution.
Besides, the energy formulation allows convection development in an slightly super-adiabatic state, which is expected to happen in an environment where the radiative flux is small.
In our numerical experiments we have obtained convective patterns, rotation and meridional circulation profiles for solar-like stars for a broad range of rotation rates. The basic properties of the flow structure are compatible with previous results of rotating convection obtained with different numerical codes (e.g., Käpylä et al. 2011b;Gastine et al. 2013). The Reynolds stress components computed from our simulations also agree with previous results in spherical geometry (e.g., Käpylä et al. 2011b) and partially agree with the mean-field theory of Lambda-effect (Rüdiger & Kitchatinov 2007).
We have explained the develpment of differential rotation as a consequence of the competition between buoyancy and Coriolis forces. This balance results in different Reynolds stress profiles which define the angular momentum distribution. For slow rotating models where buoyancy dominates, the most important contribution to the angular momentum flux comes from the meridional circulation, which results from the turbulent correlation u ′ r u ′ θ . For large Rossby numbers the results show either one (model L1) or two (models L2 and L3) meridional cells per hemisphere.
For the rapidly rotating models the correlations u ′ φ u ′ θ and u ′ φ u r dominate, and u ′ r u ′ θ might have only a marginal contribution. The differential rotation forms due to these components of the Reynolds stress and the meridional flow is formed due to inertial forces due to the rotation (Miesch & Hindman 2011). The meridional circulation exhibit several cells forming at each hemisphere. This result was found in other global convection simulations (e.g., Brun & Toomre 2002;Käpylä et al. 2011b;Ghizaru et al. 2010). However, the formation of several circulation cells is at odds with hydrodynamic mean-field models.
In order to obtain a solar-like differential rotation profile, our simulations require convective velocities of the order of m/s in the bulk of the convection zone. The Rossby number associated with this rotation rate and u rms is of the order of 10 −2 which defines the Sun as a fast rotator, see Fig. 7. Furthermore, the models close to the transition between anti-solar and solar-like differential rotation also show multiple circulation cells. Recent observational results also indicate that the Sun has more than one cell (Hathaway 2012;Zhao et al. 2012Zhao et al. , 2013. Further studies are required, however, to clarify this issue. We have also shown that the transition between anti-solar to solar-like rotation could also be obtained by changing the forcing of the system, i.e., making convection more or less vigorous (models L3, L8 and L9). However, as mentioned before, for the solar rotation rate, only the ambient state defined in Eqs. (6 to 8) results in a fair agreement with the observed properties described in Introduction as items (i) to (iv). In the model L6 we have been able to satisfy the properties (i) and (iii). An important result is the formation of the tachocline at the interface between stable and unstable layers.
In all the models presented here the rotation contours have cylindrical shape. Although at just below the tachocline a latitudinal gradient in Θ ′ is observed in the rapid rotating simulations, it is not transported to the inner convection zone, thus being insufficient to break the columnar rotation. The solar rotation property (ii) is not reproduced for our model. A different choice of the ambient state Θ e has proven to be able to create radial contours of rotation at middle an high latitudes (Racine et al. 2011). On the other hand, Warnecke et al. (2013) have obtained conical rotation profiles at lower latitudes by considering an extended domain where the solar convection zone is surrounded by a simplified coronal model. We believe that a more fundamental property (including the contribution of the magnetic field) of the plasma could be the responsible for the radial rotation profiles. A detailed study of this will be the subject of a forthcoming paper.
To tackle the formation of the near-surface shear layer, property (iv), we have constructed a model with an extended radial domain up to r = 0.985R ⊙ . In the outer 6% of the domain the entropy of the ambient state decreases rapidly with radius. The results of the model N1 show that the Rossby number increases in this thin layer to values above Ro = 0.1. This indicates that convection is vigorous in this region, evolves on a short timescale, and is only marginally affected by rotation. In this case, above r = 0.96R ⊙ the angular velocity decreases with radius at near equatorial latitudes 15 • . Due to the thin width of the layer no poleward meridional circulation (like the one observed in slow rotating models) is formed due to this negative shear. Similar results have been obtained by Käpylä et al. (2011a) and Gastine & Wicht (2012) with the use of a strong density stratification. Although we keep the adiabatic background stratification constant in the energy equation, the density constrast associated with the ambient state, ρ e , increases.
Thus, both methods seem to be equivalent. In order to reproduce a near-surface shear layer spaning over all latitudes a further extended radial domain is perhaps required. Tilted contours of rotation in the bulk of the convection zone might also facilitate the formation of this layer such as it is observed.
We have tested the convergence of our models by increasing the resolution of the model L6 by a factor of 2 and 4 in models M1 and H1, respectively. The results show that for higher resolutions the solar-like differential rotation profile disappears in the deep convection zone and is reduced at the surface. We noticed that that meanwhile the Rossby number remains roughly constant, the contribution of the Reynolds stress flux decreases and the meridional circulation flux increases. This indicates that the smaller scales resolved when the turbulence of the system is increased are modifying the mixing of angular momentum in an unexpected way. This issue has been found and studied in previous works Elliott et al. 2000;Brun & Toomre 2002), by exploring different parameter regimes and boundary conditions. Brun & Toomre (2002) found that a solar-like latitudinal difference in Ω could be achieved at higher Reynolds number by decreasing the Prandtl number. This implies a more homogeneous and efficient exchange of heat and filtering certain scales. This kind of analysis involves the use of explicit dissipative