Adam S Jermyn

Shashikumar M Chitre

Pierre Lesaffre

Christopher A Tout

Convective Differential Rotation in Stars and Planets II: Observational and Numerical Tests

Keywords: convection -Sun, rotation -stars, rotation -stars, evolution -stars, interiors

Differential rotation is central to a great many mysteries in stars and planets. In Part I we predicted the order of magnitude and scaling of the differential rotation in both hydrodynamic and magnetohydrodynamic convection zones. Our results apply to both slowly-and rapidly-rotating systems, and provide a general picture of differential rotation in stars and fluid planets. We further calculated the scalings of the meridional circulation, entropy gradient and baroclinicity. In this companion paper we compare these predictions with a variety of observations and numerical simulations. With a few exceptions we find that these are consistent in both the slowly-rotating and rapidlyrotating limits. Our results help to localize core-envelope shear in red giant stars, suggest a rotation-dependent frequency shift in the internal gravity waves of massive stars and potentially explain observed deviations from von Zeipel's gravity darkening in late-type stars.

INTRODUCTION

Driven largely by high-cadence precision photometry from the CoRoT (Roxburgh & COROT Team 1999), Kepler (Gilliland et al. 2010) and TESS (Ricker et al. 2015) missions, the ability of observations to reveal stellar rotation has increased dramatically in recent years. Asteroseismology now enables strong constraints to be placed on radial rotation profiles (Aerts et al. 2019), revealing large differential rotation in both red giant (Beck et al. 2012) and Sun-like stars (Ouazzani et al. 2019). Similarly, star spot timing has provided measurements of latitudinal shear at the surfaces of stars (Donahue et al. 1996;Reinhold & Gizon 2015;Lurie et al. 2017).

Spectral deconvolution (Donati & Collier Cameron 1997) and spectropolarimetry (Reiners & Schmitt 2003) also provide a handle on latitudinal differential rotation and, because these do not require the presence of spots, they can be applied to a wider range of stars. Importantly, the use of spectra rather than spots also removes uncertainty in the spot latitude.

The breadth of these observations of other stars complements the depth of those of the Sun, which reveal its detailed E-mail: adamjermyn@gmail.com rotation profile (Thompson 1991;Schou et al. 1998). In addition, some compact objects now provide strong constraints on the rotation profiles of their progenitors. For instance, gravitational waves allow the spins of merging black hole binaries to be measured (Kushnir et al. 2016;The LIGO Scientific Collaboration et al. 2018;Zackay et al. 2019). Assuming no significant spin changes owing to accretion or supernovae, such measurements then indicate the angular momenta of the cores of the progenitor stars. Similarly, photometry of white dwarfs provides rotation rates and hence constrains the spins of cores of lower-mass stars (Hermes et al. 2017).

The growing diversity and depth of observations makes the distribution of angular momentum a key theoretical question which has driven the development of mean field turbulence theories (Kitchatinov & Ruediger 1995), thermal wind balance arguments (Brun et al. 2010;Balbus et al. 2012) and numerical simulations (Käpylä 2011;Miesch & Toomre 2009). Importantly, rotation and differential rotation play an active role in the structure and evolution of stars. For instance by inducing mixing (Eddington 1929;Maeder 1998;Chaboyer & Zahn 1992) and generating magnetic fields (Spruit 2002) and activity (Wright & Drake 2016).

Previously (Jermyn et al. 2020, hereinafter Paper I) we provided predictions for the magnitude of differential rota- tion in the convection zones of stars and gaseous planets. These predictions are order-of-magnitude scaling relations based on considerations of the asymptotic scalings of different physical processes. Along the way we also predicted the scaling of other quantities such as the baroclinicity and meridional circulation. In this companion paper we show that our predictions are generally in good agreemnt with a variety of different observations and numerical simulations. The greatest disagreements are with the simulations, which arise primarily when these are highly diffusive, highlighting the importance of developed turbulence in angular momentum transport.

r θ z R Ω ∂ z Ω ∂ R Ω
In the next section we define our notation and review the key assumptions and results of Paper I. We then proceed to compare our predictions to observations of both radial and latitudinal shear in low mass, solar-type and red giant stars, as well as Juno measurements of latitudinal shear in Jupiter (section 3). We also compare our predictions of baroclinicity with measurements of the solar latitudinal temperature gradient. Following this we turn to differential rotation in hydrodynamic and MHD simulations in section 4, where we also examine related quantities such as the baroclinicity, convection speed and, for MHD simulations, the magnetic energy density. In section 5 we describe further tests which could be done given more observations and numerical simulations. We conclude with a discussion of the results and their astrophysical implications in section 7.

OVERVIEW

Rotation breaks spherical symmetry. This means that rotating stars may be very different from non-rotating stars. In Paper I we studied the response of convection zones to both slow and rapid rotation. We focused on systems with efficient convection zones, meaning those in which most of the heat is transported by advection and very little by diffusion. In such systems we expect the microscopic thermal diffusivity, composition diffusivity and viscosity to all be irrelevant, which allows us to specify the problem fully by giving just the geometry of the convection zone, the total angular momentum and the profile of the Brünt-Väisälä frequency

N 2 = - γ -1 γ g • ∇s, ( 1
)
where g is the acceleration due to gravity, which need not be radial, s is the dimensionless entropy

1 s = 1 γ -1 (ln P -γ ln ρ) , ( 2
)
where P is the pressure, ρ is the density, and γ is the adiabatic index.

Unfortunately, turbulence is difficult to analyze. Because of this we just attempted to understand the orders of magnitude and scaling relations involved in rotating convection. With this, the problem simplifies to one specified just by the geometry, the ratio of the angular velocity Ω to the magnitude of the Brünt-Väisälä frequency |N |, and the ratio of the pressure scale height h ≡ |∇ ln P| -1

(3) to the spherical radius r. In the same spirit we further assume that the precise geometry does not matter, and treat spherical shells and fully spherical convection zones alike. Our analysis began in the non-rotating limit, which on average must be spherically symmetric. In that limit the convection zone is parameterized by the scale height h and the Brünt-Väisälä frequency |N |. From these we then expect to determine the convection speed c , which is the root-meansquare of the difference between the velocity of a fluid element and its mean velocity, as well as the Alfvén velocity

A ≡ B √ 4πρ , ( 4
)
where B is the root-mean-squared magnetic field strength. These quantities are shown in the upper-left and lower-right panels of Fig. 2. In ordinary mixing length theory for efficient convection c ≈ h|N |, and A is usually expected to be of the same order (Cantiello & Braithwaite 2019), so we assume both of these hold in the non-rotating limit. We next turned to the limit of slow rotation (Ω |N |). In this limit rotation breaks spherical symmetry. This allows there to be steady meridional circulation currents, differential rotation, and misaligned entropy and pressure gradients, shown in Fig. 1 and the upper-right and lower-left panels of Fig. 2. For convenience in working with these quantities we defined the baroclinicity

ξ ≡ e φ • (∇ ln P × ∇s) |∇ ln P||∇s| , (5)
which is directly proportional to the thermal wind term in the vorticity equation (Part 1 equation 3) and, when small, is approximately the angle between the pressure and entropy gradients.

To determine how each of these quantities scales with Ω we used symmetry arguments to constrain the possible scalings of the different components of the turbulent Reynolds stress. So, for instance, because mapping Ω → -Ω is equivalent to mapping the azimuthal coordinate φ → -φ, we know that the Reynolds stress component T r θ must be an even function of Ω because it is even under φ → -φ. If the stress is analytic in a region around Ω = 0 this means that T r θ scales at least as fast as Ω 2 . We then assumed that each quantity scales at the lowest allowed order and obtained scalings for each of the quantities shown in Figs. 1 and2, including the shear and baroclinicity. These scalings suggest, that in steady state, the thermal wind term, mechanical forcing by turbulent stress, and losses owing to turbulent viscosity are all of the same order and serve to set the shear.

In the limit of rapid rotation we then assumed that there are no remaining symmetries, so that all components of the stress and turbulent diffusivity and so on are of the same order. We supplemented this with the scaling laws of Stevenson (1979) for the convection speed and Christensen & Aubert (2006) for the Alfvén velocity. Importantly these scalings imply that the Brünt-Väisälä frequency depends on the rotation rate in this limit, so we denote the non-rotating frequency as |N | 0 .

We were then able to determine the lowest possible order of each quantity in Ω -1 , and assumed that each declines as slowly as possible given these constraints. The resulting scaling laws indicate that the Taylor-Proudman term ∂Ω 2 ∂z balances azimuthal turbulent stresses and thereby sets the differential rotation, and that advection of angular momentum balances meridional turbulent stresses to set the meridional circulation.

To summarize, our predictions were derived under a few assumptions.

(i) Dimensionless factors arising from geometry are of order unity unless symmetries require them to be otherwise.

(ii) All external perturbing forces, such as tides or external heating, are negligible in the regions of interest.

(iii) The material is non-degenerate, compressible and not radiation-dominated.

(iv) All microscopic (i.e. non-turbulent) diffusivities are negligible, such that (a) convection is efficient, so the gas is nearly isentropic, (b) the Reynolds and Rayleigh numbers are much larger than critical, and (c) magnetohydrodynamical processes are ideal.

(v) The system is axisymmetric in a time-averaged sense.

(vi) Convection is subsonic.

(vii) The system is chemically homogeneous.

(viii) The pressure scale height h is less than or of the same order as the radius r.

Most of these predictions hold in all of the systems with which we test our predictions. The two exceptions are the assumption that convection is subsonic, which is marginal in systems with shallow surface convection zones, and that hlar, which fails to hold in near the centres of coreconvecting stars. However, when these fail to hold they usually do so over a small volume, so we do not expect either of them to prevent us from comparing our predictions with observations or numerical simulations.

OBSERVATIONAL TESTS

In Paper I we studied differential rotation in both slowly and rapidly rotating convection zones in both the hydrodynamic and magnetohydrodynamic (MHD) limits. Our results are given in Table 1. In the limit of slow rotation the scaling is the same for both hydrodynamic and MHD convection so those cases are grouped together.

It is important to emphasize that we have only predicted the scalings of these various quantities but not their actual magnitudes. We expect that in each case there are likely factors of order unity in these relations, such that for a quantity Q scaling as (Ω/|N | 0 ) α we have

Q = λ Q Q 0 Ω |N | 0 β , (6)
for some λ Q of order unity and dimensional Q 0 which gives the characteristic scale in terms of h and |N | 0 . For clarity we have non-dimensionalized all quantities appearing in Table 1, so that in each case Q 0 = 1.

To compute Ω/|N | 0 from observations we used given values where this quantity was provided directly, and in all other cases we first computed stellar models to match Table 1. The scalings of the differential rotation, meridional circulation, baroclinicity, Brünt-Väisälä frequency, convective velocity, and the ratio of magnetic to kinetic energy are given for the three regimes of interest in terms of the non-rotating Brünt-Väisälä frequency

| N | 0 . Case | R∇Ω| Ω | R ∂ R Ω| Ω | R ∂z Ω| Ω |r ∂r Ω| Ω |∂ θ Ω| Ω Slow (Ω | N | 0 ) 1 1 1 1 1 Fast Hydro.(Ω | N | 0 ) Ω | N | 0 -3/5 Ω | N | 0 -3/5 Ω | N | 0 -6/5 Ω | N | 0 -3/5 Ω | N | 0 -3/5 Fast MHD (Ω | N | 0 ) Ω | N | 0 -3/4 Ω | N | 0 -3/4 Ω | N | 0 -3/2 Ω | N | 0 -3/4 Ω | N | 0 -3/4 Case ur h | N | 0 u θ h | N | 0 ξ | N | | N | 0 c h | N | 0 2 A 2 c Slow (Ω | N | 0 ) h r Ω | N | 0 2 Ω | N | 0 2 Ω | N | 2 1 1 1 Fast Hydro.(Ω | N | 0 ) h r Ω | N | 0 -7/5 Ω | N | 0 -7/5 1 Ω | N | 0 2/5 Ω | N | 0 -1/5 N/A Fast MHD (Ω | N | 0 ) h r Ω | N | 0 -1/2 Ω | N | 0 -1/2 1 Ω | N | 0 1/4 Ω | N | 0 -1/2 Ω | N | 0 3/2
the observed characteristics of each star using the Modules for Experiments in Stellar Astrophysics ( MESA Paxton et al. 2011, 2013, 2015, 2018;Paxton et al. 2019) software instrument. We then constructed a volume-weighted average Brünt-Väisälä frequency given as

|N | 0 ≈ |N | avg ≡ exp ∫ dr ln | N | max(h,r) ∫ dr max(h,r) . (7)
A geometric average was chosen because the rotation profile we predict varies as a power-law in |N |. The weighting of dr/max(h, r) was chosen because the shear we calculate is per unit ln r, but near the cores of stars r < h and the convective cells become large relative to r. In that limit h becomes the relevant length-scale, so we weight by h instead of r.

As described in Appendix A, equation (7) is the form we have used whenever a volume average of |N | was needed. In Appendix C we study the sensitivity of our analysis to this averaging and find that, while there is a potentially large systematic offset in |N | 0 depending on the averaging used, up to a factor of 10 2 , the shape of the data are robust to different choices of averaging. An additional factor of order 2 uncertainty arises from uncertainties in matching stellar models to the observations, and particularly in uncertainties in the mixing length parameter. These effects are discussed in more detail in Appendix D.

In what follows we endeavour to test the relations in Table 1 through many different means. We use data from different sources and so must standardize them. Details of how the data from various sources were obtained, standardized and processed are provided in Appendix A, and all scripts and files needed to reproduce this analysis are available on https://zenodo.org/record/3992228.

Latitudinal Shear

We begin by considering observations of latitudinal shear |∂ θ ln Ω|, shown in Fig. 3 alongside our predictions for the MHD and hydrodynamic limits where relevant. The top panel shows observations from a variety of sources including asteroseismology but excluding starspot measurements, while the bottom panel shows only starspot data. In gen-eral these measurements should be interpreted as indicating the average shear over the surface of a star, though different techniques are sensitive to different latitudes. For instance starspot measurements are most sensitive to shear over typical active latitudes. While we have made an attempt to standardize these data, we urge the reader to see Appendix A for details on how data from each source were interpreted.

With the exception of contraints from the Juno mission, the top panel mostly probes the limits of slow and moderate rotation. While there are offsets between measurements made with different techniques, there is general agreement that the shear is of order unity for slow rotators and starts to decline roughly as predicted as the rotation rate increases. The data from Benomar et al. (2018) are notable outliers which lie well above measurements of objects at comparable rotation rates but otherwise the data seem generally consistent both internally and with our predictions.

The Juno data provide our strongest check on the rapidly-rotating limit (Kaspi et al. 2018). The surface data are more consistent with the hydrodynamic regime, while the upper limits from deeper regions are more consistent with the MHD regime. This is in agreement with previous work (Kirk & Stevenson 1987;Liu et al. 2008) and the finding that the transition between shallow and deep differential rotation in Jupiter is associated with a large change in electrical conductivity (Guillot et al. 2018). Nonetheless, even for these rapid rotation rates the difference between these scalings is only a factor of a few and this makes definitive statements about different regimes difficult.

The lower panel of Fig. 3 shows starspot measurements of latitudinal shear. The data appear consistent with our expected scaling relations. In particular the slopes of the data both of Reinhold & Gizon (2015) and Lurie et al. (2017) roughly match our predictions, though with considerable scatter which makes it difficult to say for certain, and the three points from other sources (Frasca et al. 2011;Bonanno et al. 2014;Davenport et al. 2015) scale similarly. There appears to be an overall offset of a factor of 0.1 to 0.3 in the magnitude of the shear between the data and our predictions, so we suggest using our scaling relations with a prefactor of β |R∇Ω | ≈ 0.2. This is consistent with the offset we see with the Sun, though the data of Benomar et and Ammler-von Eiff & Reiners (2012) suggest no such offset. This could be due to differences in which regions of the star are probed by different techniques or it could represent systematic differences of which we are unaware in either the different measurement techniques or the objects observed. For instance different samples may contain different fractions of fully convective stars, stars with deep convection zones, and others with shallow convection zones. While our predictions do not depend on the depth of the convection zone that does not mean that there is no such dependence.

A further test is provided by the empirical scaling relations of Balona & Abedigamba (2016) who find that

∂ θ ln Ω ∝ Ω -n , (8)
where n = 1.1 for K, 0.8 for G and 0.6 for F stars. Their sample is comprised mostly of stars with rotation periods of order 10 d. Except perhaps very near the surface, this is generally more rapid than the convective turnover time in these stars, so they are primarily in the rapid regime in which, for ionized systems, we predict n = -3/4 = -0.75, which falls in the middle of the observed range.

The trend of increasing n, and hence increasing relative shear, as stellar temperature rises could be due to a combination of the shallowing convection zone and an increasing Brünt-Väisälä frequency. As this happens the mean convective motions become faster and the stars approach the slow regime, where we predict n = 0. So it seems likely that the trend in n is a result of the transition from one regime to the other.

Radial Shear

We now turn to the (spherical) radial shear |r∂ r ln Ω|. In this case the observations come from a combination of helioseismology, asteroseismology and gravity multipole measurements of Jupiter (Fig. 4). The asteroseismology measurements are sensitive to the the total shear across the convection zone and so we have converted this into an average shear and plotted them at the average |N | 0 computed by equation (7). The helioseismic and Juno results are plotted with profiles of |N | 0 from models of the Sun and Jupiter described in Appendices A5 and A6.

As before, the fastest rotator for which we have data is Jupiter (Kaspi et al. 2018) and they show a preference in the outer layers for the hydrodynamic regime and in the inner layers for the MHD regime, though in both cases the preference is weak. In the slowly rotating regime the solar data provide our main test and are broadly consistent with our prediction of no significant scaling (Antia & Basu 2001).

The available asteroseismic measurements are mostly for slowly-rotating or moderately-rotating systems which lie close enough to the predicted break in the power-law (Ω ≈ |N |) that it is not clear which regime they probe. Nonetheless they are broadly consistent with our expectations and show slight evidence of a trend of decreasing relative shear with increasing rotation rate. Note that some caution is warranted in the interpretation of these data because the asteroseismic measurements of differential rotation probe a part of the star which often includes a stably stratified region. There is thus some uncertainty in attributing the shear purely to the convection zone, though that is what we have done here.

The solar data as well as the bounds set by Nielsen et al. (2017) and Deheuvels et al. (2015) suggest that the overall scale of our predictions is too large by a factor of 3 or so. This agrees with the offset we saw in the starspot measurements so we think it likely that the scale

β |R ∇Ω | ≈ 0.3.

Baroclinicity

In order to compute the differential rotation in convection zones we also needed to compute the angle ξ between the pressure and entropy gradients. This baroclinic angle is important because it determines the importance of the thermal wind term ∇p× ∇ρ in the vorticity equation. It is also closely related to the temperature difference between the pole and the equator of a star because it is proportional to the temperature gradient along isobars. Because of this, we may compare the observed temperature difference on the Sun to our predictions for ξ.

Unfortunately the pole-equator temperature difference on the Sun has proven difficult to measure due to significant observational uncertainties and the small expected signal (Takeda & UeNo 2017). Rast et al. (2008) summarizes a variety of efforts, most of which either conclude that the difference is zero to within uncertainties or provided a figure of order 1.5 K.

Bearing in mind these uncertainties, in Fig. 5 we have inferred ξ from the measurements made by Altrock & Can- field (1972), who found a temperature difference of 1.5 ± 0.6 K. We have translated done this with two different methods described in Appendix A5. In the first approach we assumed that ξ is determined by the local properties of the convection zone, so we have plotted the resulting data at the local Ω/|N | 0 near the surface of the solar convection zone.

In the second approach we assumed that the relative change in temperature between the pole and the equator persists throughout the convection zone and computed a volumeweighted average ξ and an average Ω/|N | 0 , with the average Brünt-Väisälä frequency given by equation (7).

Bearing in mind the large uncertainties on this measurement, we find good agreement with our predicted Ω 2 scaling for slow rotators. The solar data do not probe the rapidly rotating limit so for that we rely on simulations in Section 4.3. The scaling we find suggests that the baroclinic angle in convection zones is indeed driven by rotation-induced anisotropy in the convective heat flux (Jermyn et al. 2018b).

Magnetic Activity

Our prediction is that magnetic field strength is independent of rotation rate for Ω < |N | and increases with increasing rotation rate for Ω > |N |. This is somewhat at odds with observations of magnetic activity, which show increasing X-ray luminosity with increasing rotation rate for slowly rotating stars (Ω < |N |) and a plateau for rapidly rotating stars (Ω > |N |, Wright et al. 2011). This poses two chal-lenges. First, if our prediction is correct for slowly-rotating stars, how can tracers of magnetic activity such as X-ray luminosity increase without the magnetic field strength increasing? Secondly, if our prediction is correct for rapidlyrotating stars, how can tracers of activity plateau while the magnetic field strength increases?

We think the first challenge is most likely resolved by changes in the geometry of the magnetic field as a function of rotation rate. In numerical simulations of convective dynamos the magnetic field becomes increasingly ordered on large scales as the rotation rate increases (see e.g. Christensen & Aubert 2006). If activity relies on the magnetic field having dipolar or low-order multipolar structure then the activity can increase with rotation rate even if the magnetic field strength is constant, as we predict for slowlyrotating stars. This is supported by observations which suggest that the X-ray luminosity increases with increasing large-scale magnetic field (Vidotto et al. 2014, figure 6).

The second challenge could be explained in a few different ways. First, surface convection zones are often inefficient, meaning that most of the heat flux is carried radiatively rather than by advection. As a result increasing the rotation rate would not increase |N | as we have predicted in the case of efficient convection zones. Because we predict that A ≈ h|N | this suggests that, even when Ω |N |, the surface magnetic field strength does not increase with increasing rotation rate. The entire effect of increasing activity with increasing rotation rate would then be due to the field geometry becoming dipolar. This eventually saturates around Ω ≈ |N |, producing a plateau. Note that this explanation requires that the surface magnetic field be generated near the surface, such that it is sensitive to the properties of near-surface convection.

A further possibility is that the magnetic activity depends primarily on a spot coverage fraction which saturates and so no further increase in field strength can produce more activity (Vilhu 1984). Testing such a proposition is difficult, but if it holds it would resolve the apparent conflict between our predictions and the observations.

NUMERICAL SIMULATIONS

We now turn to numerical simulations. Because these typically report the Rayleigh number, we can directly compute |N | and so we rephrase our scaling relations in terms of the Brünt-Väisälä frequency which is actually realized in the system, rather than the non-rotating one we used in section 3. The reason that |N | is not the same as |N | 0 is that rotation generally makes convection less efficient at transporting heat. This means that the entropy gradient has to be steeper for convection to carry the same heat flux in a rotating system than in a non-rotating one.

The relation between |N | and |N | 0 is derived in detail in Paper I. Briefly, in non-magnetized convection zones, we compute the heat flux as

F ≈ ρc p T µ D|∇s|, (9)
where D ≈2 c /|N | is the turbulent diffusivity. Inserting our predicted scaling for c ∝ |N |3/2 Ω -1/2 then yields

F ≈ ρh 3 |N | 5 Ω -2 . ( 10
)
Setting F to be constant then yields |N | ∝ Ω 2/5 . The analysis is similar in magnetized convection zones, but begins with

F ≈ ρ c 2 A . (11)
Inserting our scaling for c and the scalings A ≈ h|N | yields

F ≈ ρh 3 |N | 4 Ω -1 (12)
and so |N | ∝ Ω 1/4 , as given in Table 1.

Our predictions in terms of |N | are given in Table 2. Once more we emphasize that we have only predicted the scalings of these various quantities but not their actual magnitudes. We expect that in each case there are likely factors of order unity in these relations, such that for a quantity Q scaling as (Ω/|N |) α we have

Q = λ Q Q 0 Ω |N | β (13)
for some λ Q of order unity and dimensional Q 0 which gives the characteristic scale in terms of h and |N |. For clarity we have non-dimensionalized all quantities appearing in Table 1, so that in each case Q 0 = 1. We now test the relations in Table 2. Details of how the data from various sources were obtained, standardized and processed are provided in Appendix A, and all scripts and files needed to reproduce this analysis are available on Zenodo.org.

Hydrodynamic Simulations

Because simulations tend to probe different regions of parameter space they provide a complementary set of tests. In this case simulations tend to fill in the slowly-rotating and moderately-rotating limits more than the observations and this allows us to test whether the scaling indeed changes between these limits.

The results of a variety of different hydrodynamic convection simulations performed for different geometries, with different software instruments, in both the anelastic 2 and fully compressible limits, and with different choices of dimensionless parameters are shown in Fig. 6 There is significant scatter between simulations performed by different groups. Some of this may result from the use of different software instruments, though several of the instruments used do yield identical results on identical Table 2. The scalings of the differential rotation, meridional circulation, baroclinicity, convective velocity, and the ratio of magnetic to kinetic energy are given for the three regimes of interest. Note that the latitudinal and spherical radial differential rotation are each formed of a mixture of the cylindrical vertical and radial differential rotation. Because the cylindrical radial shear is larger than the vertical shear, both spherical components of the differential rotation share the scaling of the former. problems (Jones et al. 2011) It could also be due to imperfections in our standardization of the data. In particular, different symbols in Fig. 6 correspond to different ways in which the differential rotation was reported, and so reflect different kinds of spatial averaging as well as different weightings of the radial and latitudinal components of ∇Ω. We do not see any systematic trends across different reported measures, though, so this also seems unlikely to be the cause of a majority of the scatter.

Case

| R∇Ω| Ω | R ∂ R Ω| Ω | R ∂z Ω| Ω |r ∂r Ω| Ω |∂ θ Ω| Ω Slow (Ω | N | 0 ) 1 1 1 1 1 Fast Hydro.(Ω | N |) Ω | N | -1 Ω | N | -1 Ω | N | -2 Ω | N | -1 Ω | N | -1 Fast MHD (Ω | N |) Ω | N | -1 Ω | N | -1 Ω | N | -2 Ω | N | -1 Ω | N | -1 Case ur h | N | u θ h | N | ξ c h | N | 2 A 2 c Slow (Ω | N |) h r Ω | N | 2 Ω | N | 2 Ω | N | 2 1 1 Fast Hydro.(Ω | N |) h r Ω | N | -3 Ω | N | -3 1 | N | Ω N/A Fast MHD (Ω | N |) h r Ω | N | -1 Ω | N | -1 ≈ Ω | N | -1 1 | N | Ω Ω 2 | N | 2
There is also some scatter within the data from individual sources. Different choices of boundary conditions or geometry can contribute to this scatter (see e.g. Yadav et al. 2016), but it seems likely that a large fraction of it results from the use of different Prandtl, Rayleigh or Reynolds numbers, and we expect such differences to become smaller as the latter two numbers become larger.

To probe these effects across simulations we use the Péclet number

Pe ≡ d α , ( 14
)
where α is the thermal diffusivity, d is a characteristic length scale, is the root-mean square of the velocity in the rotating frame with zero net angular momentum. When the Péclet number is sufficiently large convection is efficient because the heat transport comes to be dominated by advection rather than diffusion. We choose this number despite it being an output from simulations rather than an input because it directly measures whether or not the convection is efficient and this is a critical factor in determining whether or not our predictions are applicable.

The upper panel of Fig. 7 shows the same simulations coloured by Péclet number4 . There is a clear trend towards less shear with lower Pe. In the lower panel we retain only those simulations with Pe > 10 2 . Doing so significantly reduces the scatter and improves the agreement with our predictions.

Some of the remaining scatter is due to density stratification. In particular, the different tracks that can be seen in the data of Gastine & Wicht (2012) in Fig. 6 correspond to different degrees of density stratification. The strength of this effect decreases with increasing Péclet number. This could be because, as the density stratification increases, the velocity required to carry the heat flux becomes large near the top of the domain and small near the bottom. So when the mean Pe remains constant, increasing stratification means that an increasingly large proportion of the domain has a low Pe. As such more work may be required to probe whether this effect holds at the large Pe > 10 6 that apply in the Sun for instance.

If there is an effect of density stratification at large Pe, it could be due to the compressibility torque described by Glatzmaier et al. (2009). Such effects are potentially very physically important if they hold in the astrophysically relevant high-Pe limit particularly because some systems, such as red giants, exhibit extreme degrees of density stratification in their convection zones.

Focusing just on the high-Pe number data shown in the lower panel of Fig. 7, we see an overall trend for the relative shear to plateau towards the slow rotation limit and to fall in the rapidly rotating limit. The existence of the plateau suggests that the scaling of the off-diagonal components of the Reynolds stress is indeed as we have argued in Part I, with the rθ component scaling quadratically and the r φ and θφ components scaling linearly.

In the rapidly rotating limit our scaling appears to approximately bound the results of simulations. That is, some simulations report shears significantly below our predictions while very few report anything significantly above them and the largest shears found in simulations appear to track our scaling reasonably well.

A second test of our predictions comes from the convective velocities. Fig. 8 shows our predictions for the ratio c /h|N | alongside the results of the subset of these simulations which reported enough information for us to calculate this ratio. Note that we have used |N | rather than |N | 0 in the denominator to emphasize the suppression relative to the predictions of non-rotating mixing length theory.

We see that the data all exhibit the same trend of a plateau towards slow rotation and a decline towards rapid rotation. There is again considerable scatter between simulations performed by different groups. This is again primarily due to many of these simulations having a low Péclet number. We can see this in Fig. 9, which shows the full suite of simulations as well as cuts with Pe > 10 2 and Pe > 10 3 . As we restrict the Péclet number to be ever-larger the scatter reduces and the agreement with our predictions improves.

Note that we tend to overestimate the convection speed relative to the simulations. This is just because we have set all dimensionless constants to unity in our theory. In practice we can absorb the apparent constant offset into our predictions by setting the prefactor λ c ≈ 0.2 in our relations.

Finally, while we have focused on simulations in spherical geometries, our results also agree with the cartesian simulations of Barker et al. (2014), who find, as we do, that c ∝ Ω -1/5 at fixed heat flux.

MHD Simulations

While hydrodynamic systems provide a helpful test of our arguments, the most common astrophysical case is the nearideal MHD limit. Fig. 10 shows our predicted differential rotation scaling along with the findings of various MHD simulations of convection. As before we include both anelastic and fully compressible simulations. We employ the same conventions that filled shapes denote spherical domains, open shapes denote shellular domains, squares indicate radial shear, circles indicate latitudinal shear and triangles indicate root-mean square of the shear integrated over the domain.

There is much less scatter between simulations performed by different groups than in the case of hydrodynamic simulations. This is largely because most of our data come from simulations performed with the ASH software instrument and that these data have been reported in a relatively uniform way.

There is also less scatter between the simulations in each individual work. We believe this is partly because these simulations were performed using a smaller set of software instruments and partly because they used a more uniform set of dimensionless parameters5 . We can see the latter in a plot coloured by Péclet number, shown in the upper panel of Fig. 11. There are far fewer simulations at very low Pe than in the hydrodynamic sample and this results in less scatter. When we restrict the sample to just Pe > 10 2 (lower panel) the scatter goes away almost entirely and the remaining simulations lie nicely on our predicted slope, albeit out to slower rotation rates than we expect.

Overall we see a steeper slope than our predicted scaling in the rapidly rotating limit. This conclusion is driven mostly by the results of Yadav et al. (2013a) and Augustson et al. (2016), which suggest that |R∇ ln Ω| ∝ Ω -2 rather than our predicted Ω -1 . This is consistent with the power-law fits provided by Yadav et al. (2013a) for the zonal (azimuthal) and non-zonal (meridional) Rossby number6 . We may translate into predictions for the differential rotation at fixed heat flux 7 , suggesting that |R∇ ln Ω| ∝ Ω -1.328 , Ω -1.449 , Ω -1.210 or Ω -1.2911 . Yadav et al. (2013b) similarly find a scaling of Ω -1.35 . Such steep power laws appears inconsistent with the observations shown in Figs. 3 and4 so we remain uncertain as to whether the slope in nature is truly steeper than our pre-

7 |R ∇ ln Ω | ≈ Ro zonal . Using Ra * Q ∝ Ω -3
we then obtain the scaling for |R∇ ln Ω |. 8 Using their fit for dipolar dynamos with no Pm dependence. 9 Using their fit for dipolar dynamos with a Pm dependence. 10 Using their fit for multipolar dynamos with no Pm dependence. 11 Using their fit for multipolar dynamos with a Pm dependence. diction. However, if this is indeed what happens it points to a difficulty in our analysis. In Section 7 of Paper I we assumed that each of the shear, circulation and baroclinicity are as large as allowed by the conditions of heat and momentum balance. While it must be that in each equation at least one of these saturates its bounds, it need not be the case that they all do. If they do not then our predictions in Tables 1 and 2 are really upper bounds, only one of which must be saturated in any given scenario. So for instance it could be that the cylindrical vertical differential rotation ∂ z Ω saturates our bound and scales like Ω -2 , similar to what we see in Fig. 10. In that case there is no formal requirement that the radial shear ∂ R Ω also saturates its bound, in which case the shear only needs to be as large as

∂ z Ω ∝ Ω -2 .
In the slowly rotating regime the simulations of Varela et al. (2016) hint at a plateau but those of Augustson et al. (2016) instead show a continuing increase of relative shear with decreasing rotation rate. This could indicate that we have misplaced the break point between slow and rapid rotation. On the other hand, all of the simulations which show this trend were performed on a spherical domain and so it could be that, even though the convective turnover is fast relative to rotation on average, deep in the core where convection is slow we might still be in the rapidly rotating limit.

A further test of our predictions comes from the convective velocities. Fig. 12 shows our predictions for the ratio c /h|N | alongside the results of the subset of these simulations for which we were able to calculate this ratio.

We see that the data generally exhibit the same trend of a plateau towards slow rotation and a decline towards rapid rotation. Most of the scatter in these results is directly attributable to not being in the limit of efficient convection. We can see this in the lower panel where we have filtered for Pe > 10 2 and see significantly better agreement with our predicted slope and plateau.

Other authors have found similar, though typically weaker, scalings of MHD convection speed at constant heat flux. We predict c ∝ Ω -1/2 , while others find Ω -0.23 to Ω -0.29 (Christensen & Aubert 2006) 2013a) 13 and Ω -0.32 to Ω -0.47 (Aubert et al. 2017) 14 . This weaker scaling could be a result of the Péclet number dependence we have noted, though of course it could also be that our theory predicts a steeper scaling law than actually occurs.

Finally we examine the ratio of magnetic to kinetic energy, shown as a function of Ω/|N | in Fig. 13. Our predictions are shown as a solid line and the results of simulations are shown as a variety of shapes, following the previous convention that solid shapes indicate spherical domains, open shapes indicate shellular domains, and triangles denote measurements averaged over the simulation domain. Note that our predictions for this energy ratio are not new: they are identical to those made by several authors including Starchenko & Jones (2002) and Augustson et al. (2016).

There is significant scatter between simulations performed by different groups but within the sets of simulations with the most dynamic range in Ω/|N | we see a plateau towards slow rotation and scaling similar to that of our predictions for rapid rotation. 13 From their multipolar fit with magnetic Prandtl number dependence.

14 See their figure 11(a). Their is proportional to Ω -3 . Our predicted scalings are also similar to those that found by several groups in fits to their own numerical simulations. For instance at fixed heat flux we predict that the magnetic field strength scales as B ∝ A ∝ Ω 1/4 and in simulations this scaling has been found to be Ω -0.02 (Christensen & Aubert 2006) 15 , Ω -0.11 (Yadav et al. 2013a) 16 , Ω 0 (Yadav et al. 2013a) 17 , and between Ω -0.02 and Ω 0.

Baroclinicity

As described in appendix A, we have inferred ξ from the reported results of various hydrodynamic and MHD simulations. Our prediction is the same for both kinds of simulations and so we analyze them together. Some authors reported the temperature difference between different latitudes in a model while others reported latitudinal gradients of closely related quantities like the Nusselt number. The results of our inference and standardization are shown in Fig. 14.

We find reasonable agreement with our predicted Ω 2 scaling for slow rotators and, as expected, we see ξ saturating in the limit of rapid rotation. There is considerable scatter among simulations in the saturated regime and even considerable disagreement amongst simulations performed by the same group with the same code. This may reflect a sensitivity of the saturated ξ to the geometry or other dimensionless parameters in the simulation.

These data support our conclusion in Part I that the thermal wind term is one of the dominant terms in both the hydrodynamic and MHD regimes and in the limits of 15 See their equation ( 33). Note that Lo ∝ A /Ω and their Ra * Q ∝ Ω -3 , such that Lo ∝ (Ra * Q ) ( 0.34) means that A ∝ Ω 1-0.34×3 . 16 From their dipolar fit with no magnetic Prandtl number dependence. 17 From their dipolar fit with magnetic Prandtl number dependence. 18 See their figure 11(a). Their is proportional to Ω -3 and their λ ∝ BΩ -1 .

both slow and rapid rotation, in agreement with the arguments of Balbus & Schaan (2012). The scaling we find suggests that the baroclinic angle in convection zones is indeed driven by rotation-induced anisotropy in the convective heat flux (Jermyn et al. 2018b).

FUTURE TESTS

We have tested our predictions under many different circumstances but there are further tests that would be fruitful in the future. First, while we have tested our predictions for the latitudinal and spherical radial shear, we have not made a quantitative test of our prediction that the cylindrical vertical shear ∂ z Ω is suppressed relative to the radial component. This is seen in simulations (see e.g. Mabuchi et al. 2015) but so far as we are aware it has not been quantified. It ought to be straightforward to add ∂ R Ω and ∂ z Ω to the output of future simulations or even to run previous simulations again with these additional outputs.

Secondly, it would be good to test our predictions for the meridional circulation velocities u r and u θ . There have been measurements of the meridional circulation in the Sun (Zhao et al. 2013;Rajaguru & Antia 2015;Schad & Roth 2020) but there are still significant disagreements between the different inversion techniques which make a direct comparison to our theory challenging. On the other hand, such a comparison could be done with simulations by reporting the time-and azimuthal-average of the meridional velocity field. As far as we are aware this has not been reported19 but we see no particular barrier to doing so.

More broadly, further tests from simulations primarily face the challenge of reaching high enough Péclet numbers to be applicable to efficient convection. The primary limitation to this is the availability of computing resources, so we hope that more simulations will be performed at Pe ≈ 10 3 to 10 5 in the near future.

Finally, the most useful observational comparisons are those that probe the most rapidly-rotating systems. These are the systems which are best able to test the slopes we have predicted and to thereby either verify or falsify our theories. In this regard data such as the differential rotation measurements of Davenport et al. (2015) are extremely useful because, while they report observations of a single system, it is well-characterized and extremely rapidly rotating, particularly relative to its slow convective turnover time. This allows it to more strongly constrain the scaling of differential rotation.

PRIOR PREDICTIONS

In Part I we made predictions for the scaling of the differential rotation, meridional circulation, magnetic field and baroclinic angle. We are aware that two of our predicted scalings have been suggested previously. These are that for the magnetic field (Starchenko & Jones 2002) and that of convection speed with rotation rate the hydrodynamic limit (Stevenson 1979). The remaining scalings are new.

However, there have been many previous predictions of the same quantities and it is worth comparing our predictions with these. Our prediction that |R∇ ln Ω| is a constant of order unity for slowly-rotating systems is consistent with predictions based on the Λ effect (e.g. Küker et al. 2019, and references therein). Gilman & Foukal (1979) predicted that for slowly-rotating systems the angular momentum of fluid parcels is conserved in convection and so that Ω ∝ r -2 , or |R∇ ln Ω| = 2, consistent with our predictions. Kissin & Thompson (2015a) make the same prediction of Ω ∝ r -2 for slowly-rotating systems and provide two possible scalings for rapidly-rotating ones, |R∇ ln Ω| = α, 1 < α < 3/2 and Ω ≈ Ω 0 + |N | for some suitable choice of Ω 0 . Both of these are inconsistent with our predictions, particularly the latter, which introduces a new characteristic frequency-scale Ω 0 and predicts a characteristic length-scale for the shear of (∂ r ln |N |) -1 ≈ h.

Our prediction for the scaling of the convection speed ( c ≈ h|N | 2 Ω -1 ) is based on the theory of Stevenson (1979), and has been suggested subsequently along various different lines by Starchenko & Jones (2002), Christensen (2010) and Barker et al. (2014). In the hydrodynamic case, our prediction for the scaling of |N | with Ω matches that of Barker et al. (2014) as well. We previously predicted a somewhat different scaling (Jermyn et al. 2018a) Stevenson (1979), Jermyn et al. (2018a) did not impose the lower bound on the vertical wavenumber. Such a lower bound is physically motivated for stars by the finite scale height, so we favour the scaling of Stevenson (1979) here.

of c ≈ h|N | 3/2 Ω -1/2 because, unlike
Our prediction for the scaling of the the magnetic field with rotation rate

( A ≈ h|N | ≈ h|N | 3/4 0 Ω 1/4
) is based on the relationship between heat flux and magnetic field by Christensen & Aubert (2006) and the bound that if A > ∼ h|N | then the magnetic field begins to inhibit the linear convective instability (Gough & Tayler 1966). This prediction agrees with the theory of Stevenson (1979), andStarchenko &Jones (2002) predict both the same scaling for A and for |N |. However, a variety of other scaling laws have been proposed. These are summarized in Section 4.2. They generally agree that, at fixed |N |, the magnetic field strength varies as a small power of Ω ranging from -0.02 to 0.19, all of which is somewhat smaller than we have predicted.

DISCUSSION

Up to this point we have focused on testing our predictions against observations and simulations. We have found general agreement on the existence and location of plateaus, as well as on the trends away from plateaus. Importantly, in section 4, we found that simulations show a strong dependence on Péclet number and that they only converge to our predictions as Pe becomes large in the limit of efficient convection. Because this is the most common limit in astrophysical systems it matters that simulations attempting to reproduce observed differential rotation operate in this regime.

We now turn to exploring broader astrophysical impli-cations. We have predicted that the relative shear is of order unity for slowly rotating convection zones and falls off quickly once the rotation becomes faster than the convective turnover frequency. In particular, the relative shear per unit log r is at most of order |N |/Ω in the rapidly-rotating limit. A consequence of this is that the shear is approximately bounded by the Brünt-Väisälä frequency, which is of order the convective turnover frequency. This is because if

Ω > |N | then the absolute shear per unit log r is |N | while if Ω < |N | it is Ω, which is less than Brünt-Väisälä frequency |N |.
Hence we expect the maximum absolute shear across a convection zone to be

∆Ω max ≈ ∫ convection zone |N |d ln r. (15)
Note that this upper bound is robust even if, as suggested by Fig. 10, the differential rotation falls off faster with Ω than we have predicted.

Because, in the MHD regime, |N | deviates only weakly from |N | 0 we may estimate the right-hand side of equation (15) using |N | 0 from stellar models. For the red giants studied by Deheuvels et al. (2015) ∆Ω max ≈ 4×10 -6 s -1 . This is on the order of what Cantiello et al. (2014) and Kissin & Thompson (2015b) suggest is needed to match asteroseismic inferences of core-envelope shear in red giants, so an order unity pre-factor in equation ( 15) would suffice to bring our bound into agreement with their calculations, and it seems possible that most or all of the differential rotation in such stars is in their convective envelopes. However, because equation ( 15) is an upper bound, it remains possible that a substantial component lies in the radiative zone (Fuller et al. 2019) and some astereoseismic evidence indicates that there is more shear in the radiative zones of red giants than in their convection zones (Mauro et al. 2018).

A further prediction is that rapidly-rotating convecting stars ought to exhibit significant equator to pole temperature differences. These differences are due to the Coriolis effect acting on convective motions and not the von Zeipel effect. This has been seen in simulations (Raynaud et al. 2018) andobservations (Che et al. 2011), and our scaling relations complement these by providing a theoretical curve with which to interpret the observations. In particular, in the usual notation,

T eff ∝ g β eff , (16)
where g eff is the effective acceleration due to gravity and centrifugal effects and β = 1/4 is the von Zeipel exponent (von Zeipel 1924). Differentiating this with respect to latitude yields

∂ ln T eff ∂θ = β ∂ ln g eff ∂θ ≈ β Ω 2 r cos θ g . ( 17
)
Our findings (equation A108) suggest instead that

∂ ln T eff ∂θ = ξ r |N | 2 g eff . ( 18
)
For slowly-rotating systems, ξ ∝ Ω 2 /|N | 2 and this becomes

∂ ln T eff ∂θ ∝ Ω 2 r g eff , (19)
which is equivalent to the von Zeipel result, possibly with a different β. For rapidly-rotating systems with Ω > |N |, ξ saturates and

∂ ln T eff ∂θ ≈ r |N | 2 g eff . ( 20
)
This saturation could be why observations of late-type stars have inferred β smaller than 1/4: if the systems enter the saturated regime then a fit to the unsaturated scaling would tend to estimate a lower β (Che et al. 2011).

For extremely rapid rotators we also predict that convection is suppressed. In stars this leads to an increase in the Brünt-Väisälä frequency N and thence a steepening of the entropy gradient to carry the nuclear luminosity. This is likely to be difficult to observe directly in measurements such as the stellar radius or effective temperature because the effect is even smaller than centrifugal expansion (see appendix B). On the other hand, it could be observed in asteroseismic measurements. The cores of massive stars likely emit internal gravity waves, with a peak frequency thought to be proportional to the Brünt-Väisälä frequency (Rogers et al. 2013;Couston et al. 2018). Moreover convection in these stars is so slow that nearly all of them are in the rapidly rotating limit. Hence if our predictions hold then there ought to be a systematic shift in the peak frequency of these waves as a function of rotation that scales as Ω 1/4 . Finally, note that in the very rapid limit we expect the Taylor-Proudman constraint to hold. That is we predict that the star forms nested cylinders of constant rotation because |∂ z Ω| |∇Ω|. We could not find data to probe this prediction quantitatively but essentially all of the simulations we examined show this feature in the rapidly rotating limit.

APPENDIX A: DATA PROCESSING

Here we explain the way in which data were processed to produce our figures. All scripts and inlists as well as auxiliary data files used in this data processing are available on https: //doi.org/10.5281/zenodo.3992227.

In many cases the data processing required stellar models. These were produced using revision 11701 of the Modules for Experiments in Stellar Astrophysics MESA. The MESA equation of state (EOS) is a blend of the OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), PTEH (Pols et al. 1995), HELM (Timmes & Swesty 2000) and PC (Potekhin & Chabrier 2010) equations of state. Radiative opacities are primarily from OPAL (Iglesias & Rogers 1993, 1996), with low-temperature data according to Ferguson et al. (2005) and the high-temperature, Compton-scattering dominated regime by Buchler & Yueh (1976). Electron conduction opacities are according to Cassisi et al. (2007). Nuclear reaction rates are from JINA REA-CLIB (Cyburt et al. 2010) plus additional tabulated weak reaction rates (Fuller et al. 1985;Oda et al. 1994;Langanke & Martínez-Pinedo 2000). Electron screening is included via the prescription of Chugunov et al. (2007). Thermal neutrino loss rates are according to [START_REF] Itoh | [END_REF]. Mixing length theory was implemented following [START_REF] Cox | Principles of stellar structure Cyburt R. H[END_REF] with mixing length parameter α = 2. Models were created on the pre-main sequence and evolved from there. All other parameters were set to their defaults.

A1 Latitudinal Shear

We describe the observations.

A1.1 Benomar Observations

Data were taken from tables S2 and S3 of Benomar et al. (2018), who inferred differential rotation through asteroseismology for several solar-type stars and reported it in the form of the latitudinal shear and mean rotation rates. The reported shear Ω 1 is 2/15 times the difference between the pole and the equator. There are π/2 radians between the pole and the equator so we approximate

|R∇Ω| ≈ |∂ θ Ω| ≈ 2 π 15 2 Ω 1 . (A1)
Uncertainties in both the shear and mean rotation rate were reported. These were propagated into the relative shear by the formula

dµ = j ∂ µ ∂ x j 2 dx 2 j , ( A2
)
where µ is a function of interest and x j are variables on which it depends. This amounts to the assumption that errors are uncorrelated and small relative to the scale over which higher-order derivatives are relevant. Stellar models were constructed which matched the ages, metallicities and masses provided in their table S2. The average |N | 0 (Brünt-Väisälä frequency) in the convection zone of each model was calculated following equation (7).

A1.2 Bazot Observations

Data were taken from tables 1 and A3 of Bazot et al. (2019) who inferred differential rotation through asteroseismology for several solar-type stars and reported in the form of the latitudinal shear and mean rotation rates. The reported shear Ω 1 is 2/15 times the difference between the pole and the equator. There are π/2 radians between the pole and the equator so we approximate

|R∇Ω| ≈ |∂ θ Ω| ≈ 2 π 15 2 Ω 1 . (A3)
Uncertainties in both the shear and mean rotation rate were reported for 16 Cyg A. For 16 Cyg B a weighted average of the three modes was used following the recommendation in the caption of table 1. The uncertainty of this average was computed as the weighted root-mean-squared difference between the three modes and the average.

Uncertainties were propagated into the relative shear equation (A2). Stellar models were constructed which matched the ages, metallicities and masses provided in their table A3, and the average |N | 0 in the convection zone of each model was calculated according to equation (7).

A1.3 Ammler-von Eiff Observations

Data were taken from table 2 of Ammler-von Eiff & Reiners (2012). The latitudinal differential rotation was spectroscopically inferred for several A to F stars and reported in the form of the coefficient η/ √ sin i, where η is the pole to equator shear divided by the equatorial angular velocity. For each object our hypothesis is that the latitudinal differential rotation arises owing to the surface convection zone. So we identify

|R∇Ω| ≈ |∂ θ Ω| ≈ 2 π |Ω pole -Ω equator |. (A4)
Stars without a surface convection zone were filtered out of the sample. Many systems were consistent with zero shear to within the error bars. So we only include those with at least a 1σ detection of shear.

The rotation rate is reported as sin i, where is the surface rotation velocity. Inclinations were not reported. So rotation speeds are considered to be lower bounds and shear measurements to be upper bounds.

Eighty stellar models were constructed with metallicity Z = 0.02 and spaced uniformly in mass from 0.2 M to 1.4 M . The models were evolved to an age of (M/M ) -2.5 GYr or 10 GYr, whichever was smaller. For each object the model which best-matched the reported T eff in their table 5 was used to compute Ω = /R and |N | avg .

A1.4 Frasca Observations

Data were taken from Frasca et al. (2011). The differential rotation was inferred from star spot observations of the K star KIC 8429280 and reported in the form of the pole to equator shear, which we treat with equation (A4). The rotation rate is reported as sin i, where is the surface rotation velocity. The inclination is also reported and so we use it to compute . We employ the same grid of main-sequence stellar models from appendix A1.3 to compute |N | avg and R and use these to calculate Ω = /R and Ω/|N | avg .

A1.5 Bonanno Observations

Data were taken from Bonanno et al. (2014) who report differential rotation inferred from star spot observations of the late F star KIC 5955122. The differential rotation was reported in the form of a pole to equator shear, which we treat with equation (A4). We normalise this to the mean rotation rate and take the result to be approximately |R∇Ω|/Ω. A stellar model was constructed which matched the age, metallicity and mass reported for KIC 5955122 and we used this to compute |N | avg and Ω/|N | avg .

A1.6 Donati Observations

Data were taken from Donati et al. (2008). The differential rotation was inferred for a sample of M-dwarfs from spectropolarimetry and reported as the surface equator to pole difference in Ω normalised against that of the equator. We treat this with equation (A4). Donati et al. (2008) report the Rossby number in the form of τ c P rotation , where τ c = 1/|N | 0 (Gilman 1980) was determined from an empirical calibration to non-rotating stellar models (Kiraga & Stepien 2007) and P rotation is the average spin period of the star. So we take Ω

| N | 0 ≈ 2πRo -1 .

A1.7 Reinhold Observations

Data were taken from Reinhold & Gizon (2015), who report differential rotation inferred from star spot observations of Kepler main-sequence and subgiant stars. Stars which were flagged as having highly stable periods were excluded as possible binary or pulsator systems.

The differential rotation was reported in the form of the coefficient η/ √ sin i, which we handle as in appendix A1.4. We employ the same grid of main-sequence stellar models from appendix A1.3 to compute |N | avg and R and use these to calculate Ω = /R and Ω/|N | avg . Objects hotter than the hottest model in the grid or cooler than the coolest model in the grid were excluded from the sample.

A1.8 Lurie Observations

Data were taken from Lurie et al. (2017), who report differential rotation inferred from star spot observations of Kepler eclipsing binary systems. They infer the parameter η by

η = P max -P min P max , (A5)
where the periods P min and P max are the first minimum and maximum in the periodogram of the lightcurve of the system. They report these periods and so we compute η using equation (A5) and translate this into a shear using equation (A4).

Triple-star systems were excluded by means of the catalogues of Rappaport et al. (2013) and Borkovits et al. (2016) on the grounds that the inference of Lurie et al. ( 2017) is based on the assumption of a binary system. We further exclude the false positive systems described by Borkovits et al. (2016) because these have unusual light curves which could interfere with the inference of rotation periods.

Objects were matched by effective temperature to the grid of main-sequence stellar models from appendix A1.3, from which we computed |N | avg and Ω/|N | avg . The effective temperatures of these objects were obtained from the Kepler Eclipsing Binary catalog (Matijevič et al. 2012;Conroy et al. 2014;LaCourse et al. 2015;Abdul-Masih et al. 2016). Objects with no temperature in the catalogue were excluded. Objects hotter than the hottest model in the grid or cooler than the coolest model in the grid were excluded from the sample.

A1.9 Davenport Observations

Data were taken from Davenport et al. (2015), who report differential rotation inferred from star spot observations of the M-dwarf GJ 1243. The spots reveal a shear of 0.0047 rad d -1 . After accounting for spot latitude, they report a pole-equator shear of 0.012 ± 0.002 rad d -1 . The rotation period of the planet is reported as 0.592596 ± 0.00021 d. We convert this into η and then use equation (A4) to obtain a relative shear. Uncertainties were propagated along the way with equation (A2). The Rossby number is also reported, as by Donati et al. (2008) and we handle it in the same manner as in appendix A1.6.

A2 Radial Shear

We describe the observations.

A2.1 Deheuvels Observations

Data were taken from Deheuvels et al. (2015). The differential rotation was seismically inferred for several red giants and reported in the form of the difference between the core and envelope rotation rates divided by the envelope rotation in dimensionful entropy across the simulation domain, L the vertical size of the domain, g is the acceleration owing to gravity, c P is the specific heat at constant pressure, ν is the viscosity and α is the thermal diffusivity.

The ratio ∆S/c P is related to the dimensionless entropy change ∆s by

∆s = ∆S c P γ γ -1 . (A22) Identifying ∆s ≈ L∂ r s, (A23)
we then obtain

r∂ r s ≈ Ra να gL 4 γ γ -1 . (A24)
They further let the Prandtl number Pr be ν/α and the Taylor number Ta be 4Ω 2 L 4 /ν 2 . With these we obtain

r∂ r s ≈ 4 γ γ -1 Ra Pr Ta Ω 2 r g . (A25) So ξ ≈ γ -1 4γ Pr Ta Ra g Ω 2 r ∂ θ s. (A26)
Assuming that isobars are nearly spherical, we see that

∂ θ s ≈ ∂ θ ln T ∂s ∂ ln T P . ( A27
)
With equation ( 2) we find

∂s ∂ ln T P = γ γ -1 , (A28) so ξ ≈ 1 4 Pr Ta Ra g Ω 2 r ∂ θ ln T (A29)
Augustson et al. ( 2012) provide g, R, Pr, Ta and Ra. They also provide a quantity Ω/Ω , which we take to be the mean angular velocity of the simulation normalized by the typical solar angular velocity of 4.3 × 10 -7 cycles s -1 . We then estimate

∂ θ ln T ≈ 3 π ∆T T base CZ , ( A30
)
where we use the temperature T base CZ at the base of the convection zone because the temperature variation is dominated by that region. We are also interested in the convection speed normalized by c /h|N |. To obtain this we write

A3.3 Matt 2011 Simulations

These data come from three-dimensional anelastic hydrodynamic simulations of convection in a rotating spherical shell domain (Matt et al. 2011). The differential rotation is reported as the difference in angular velocity between the equator and a latitude of π/3, normalized to the equatorial angular velocity, so we use equation (A11) to compute |R∇ ln Ω|. The Rossby number actually realised in the flow is also reported, so we follow the procedure we used in appendix A3.1 to calculate Ω/|N |.

The temperature difference between the equator and a latitude of π/3 is also reported. We follow the same procedure as in appendix A3.2 to compute ξ from this, except that we compute

g = GM r 2 , ( A35
)
where G is Newton's constant and M is the mass of the star because Matt et al. (2011) do not report g but do report M and R.

We compute the ratio c /h|N | with the Rossby, Prandtl, Taylor and Rayleigh numbers as in appendix A3.2. Finally, we obtain the Péclet number from equation (A17).

A3.4 Brown 2008 Simulations

These data come from simulations of three-dimensional anelastic hydrodynamic convection in rotating spherical shells, reported in tables 1 and 2 of Brown et al. (2008). The differential rotation is reported as the kinetic energy of differential rotation, defined to be the mass-weighted average of the squared difference between the rotation rate and the mean rotation rate. That is, they report ) where ... denotes the mass-weighted average and Ω is the mean rotation rate. They further report the convective kinetic energy

DRKE ≈ ρR 2 Ω -Ω 2 , ( A36
CKE ≈ ρ 2 c . (A37)
To ensure consistent normalization we estimate c from their reported Rossby number as c ≈ 2ΩRo, (A38) so that the shear is

|R∇Ω| ≈ c R DRKE CKE ≈ 2ΩRo DRKE CKE , ( A39
)
where Ro is the Rossby number actually realised in the flow and of the form of equation (A18). We calculate Ω/|N | from Ro as in appendix A3.1 We compute the ratio c /h|N | with the Rossby, Prandtl, Taylor and Rayleigh numbers as in appendix A3.2. The one difference is that the Rayleigh number is reported as in equation (A14). So we must divide it by (γ -1)/γ before using equation (A34). As before, we obtain the Péclet number from equation (A17). . Different choices of averaging and weighting produce very different overall scales for |N | 0 , by up to a factor of 10 2 , but do a good job at preserving the shape of the data (i.e. slope of shear versus Ω). Because we are primarily comparing our theory with the shape of the data, and admit uncertainty as to the transition between the slowly-rotating and rapidly-rotating limits, the precise choice of averaging therefore makes only a small difference to our comparisons.

APPENDIX D: SENSITIVTY TO STELLAR MODELLING

In addition to the systematic uncertainty in |N | 0 arising from uncertainties in averaging and weighting, there are uncertainties involved in matching stellar models to the observed characteristics of stars. In matching main-sequence models to observations we have relied on the observed effective temperatures, and for giants with asteroseismic observations we have directly used the inferred masses, ages and metallicities.

These matching exercises are subject to several uncertainties. First, different stellar models constructed with different software instruments typically disagree on observable properties at the 1 -5 per-cent level even with identical input physics (Stancliffe et al. 2015). Next, input physics like opacities are uncertain at the 10 -20 per-cent level (Bailey et al. 2015). In addition different choices of age on the mainsequence, or of Helium abundance, can also lead to changes in effective temperature on this level.

Most importantly, though, 1D stellar models rely on mixing length theory, which is a highly simplified prescription for convection. One known source of uncertainty in this theory is the mixing length parameter α MLT , which is seen to vary by as much as 20 per-cent between stars which differ by 0.2M (Joyce & Chaboyer 2018). Naively extrapolating this over the range of stellar masses we have considered, this perhaps leads to a factor of 2 uncertainty in α MLT . In limit of efficient convection, |N | 0 ∝ α -2/3 MLT (see e.g. Cantiello & Braithwaite 2019, and note that |N | 0 ∝ √ ∇ -∇ ad ∝ √ ∇ -∇ ), so the dominant source of uncertainty in |N | 0 would seem to be that in the mixing length, and more broadly in using such a simplified theory for convection in the first place.

This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .

 1 Figure 1. The rotation, coordinate system and differential rotation are shown schematically: (top) the mean angular velocity Ω; (upper-right) the cylindrical radius R, vertical direction along the rotation axis z, spherical radius r and polar angle θ; (upper-left) an example of cylindrical radial differential rotation (∂ R Ω); (lower) an example of cylindrical vertical differential rotation (∂ z Ω).

Figure 2 .

 2 Figure 2. Key concepts in our theory are shown schematically: (upper-right) the meridional circulation components u r and u θ ; (upper-left) turbulent eddies move at the convection speed c and with the Brünt-Väisälä frequency N ; (lower-left) surfaces of constant pressure and entropy meet at an angle of approximately the baroclinicity ξ, resulting in surface temperature variations; (lower-right) a magnetic field with Alfvén velocity A .

Figure 3 .

 3 Figure 3. The relative latitudinal shear |∂ θ ln Ω | is shown as a function of Ω/ | N | 0 for observed convecting stars and Jupiter alongside our prediction, which is constant for Ω < | N | 0 and scales as (Ω/ | N | 0 ) -3/4 (MHD) and (Ω/ | N | 0 ) -3/5 (hydrodynamic) for Ω > | N | 0 . Shapes indicate the origin of the data: squares are for asteroseismology, triangles are for spectroscopy, stars are for spectropolarimetry and circles are for star spot measurements. For clarity the star spot measurements are shown separately in the lower panel. Helioseismic data and constraints from Juno are shown as solid lines. Juno measurements of the shear at a depth of less than about 3000km are shown separately from Juno upper limits on that at greater depths.

Figure 4 .

 4 Figure 4. The relative radial shear |r ∂ r ln Ω | ≈ |R ∇ ln Ω | is shown as a function of Ω/ | N | 0 for observed convecting stars and Jupiter alongside our prediction, which is constant for Ω < | N | 0 and scales as (Ω/ | N | 0 ) -3/4 (MHD) and (Ω/ | N | 0 ) -3/5 (hydrodynamic) for Ω > | N | 0 . Squares denote asteroseismic results. Helioseismic data and constraints from Juno are shown as solid lines. Juno measurements of the shear at a depth of less than about 3000km are shown separately from Juno upper limits on that at greater depths.

Figure 5 .

 5 Figure5. The angle ξ between the pressure gradient and the entropy gradient is shown from our predictions. It was also inferred for the Sun from the equator-pole temperature difference and the resulting upper bound is shown with a circle for the surface and a triangle for the inferred average(Altrock & Canfield 1972). Our prediction scales as (Ω/ | N | 0 ) 2 for slow rotation (Ω < | N | 0 ) and is constant for rapid rotation (Ω > | N | 0 ).

  . Squares indicate radial shear, circles indicate latitudinal shear and triangles indicate the root-mean square of the shear integrated over the domain. Filled shapes indicate simulations performed on deep spherical shell domains 3 , which are meant to mimic full spheres, while open shapes indicate those performed in thinner shellular domains.

  al. (2012) Augustson et al. (2013) Augustson et al. (2016) Aurnou et al. (2007) Brown et al. (2008) Browning et al. (2004) Brun et al. (2009) Brun et al. (2017) Featherstone et al. (2015) Gastine et al. (2012) Gastine et al. (2013) Gilman (1979-79) Guerrero et al. (2013) Käpylä et al. (2011) Kaspi et al. (2009) Mabuchi et al. (2015) Matt et al. (2011)Rogers (2015)

Figure 6 .

 6 Figure 6. The relative differential rotation |R∇ ln Ω | is shown alongside our prediction as a function of Ω/ | N | for a variety of hydrodynamic convection simulations. Squares indicate radial shear, circles indicate latitudinal shear and triangles indicate the root-mean square of the shear integrated over the domain. Filled shapes indicate simulations performed on deep spherical shell domains while open shapes indicate those performed in shellular domains. Our prediction is constant for slow rotation (Ω < | N |) and scales as (Ω/ | N |) -1 for rapid rotation (Ω > | N |).

Figure 7 .

 7 Figure 7. The relative differential rotation |R∇ ln Ω | is shown alongside our prediction as a function of Ω/ | N | for a variety of hydrodynamic convection simulations, coloured by Péclet number Pe. The upper panel shows data from all simulations for which we could compute Pe while the lower shows just those simulations with Pe > 10 2 . The slopes of the power-law predictions are indicated in the relevant regimes.

Figure 8 .

 8 Figure 8. The normalized convection speed c /h | N | is shown alongside our predictions as a function of Ω/ | N | for a variety of hydrodynamic convection simulations. Filled shapes indicate simulations performed on spherical domains while open shapes indicate those performed in shellular domaints. Triangles denote outputs which were averaged over the simulation domain. Our prediction is constant for slow rotation (Ω < | N |) and scales as (Ω/ | N |) -1 for rapid rotation (Ω > | N |).

Figure 9 .

 9 Figure 9. The normalized convection speed c /h | N | is shown alongside our predictions as a function of Ω/ | N | for a variety of hydrodynamic convection simulations, coloured by Péclet number Pe. The upper panel shows data from all simulations for which we could compute Pe, the lower shows just those simulations with Pe > 10 2 , and the lower panel shows just those with Pe > 3 . Our prediction is constant for slow rotation (Ω < | N |) and scales as (Ω/ | N |) -1 for rapid rotation (Ω > | N |).

Figure 10 .

 10 Figure 10. The relative differential rotation |R∇ ln Ω | is shown alongside our predictions as a function of Ω/ | N | for a variety of MHD convection simulations. Circles indicate latitudinal shear and triangles indicate root-mean-squared shear integrated over the domain. Filled shapes indicate simulations performed on spherical domains while open shapes indicate those performed in shellular domaints. Our prediction is constant for slow rotation (Ω < | N |) and scales as (Ω/ | N |) -1 for rapid rotation (Ω > | N |).

Figure 11 .

 11 Figure 11. The relative differential rotation |R ∇ ln Ω | is shown alongside our predictions as a function of Ω/ | N | for a variety of MHD convection simulations, coloured by Péclet number Pe. The upper panel shows data from all simulations for which we could compute Pe while the lower shows just those simulations with Pe > 10 2 . Our prediction is constant for slow rotation (Ω < | N |) and scales as (Ω/ | N |) -1 for rapid rotation (Ω > | N |).

Figure 12 .

 12 Figure 12. The normalized convection speed c /h | N | is shown alongside our predictions as a function of Ω/ | N | for a variety of MHD convection simulations, coloured by Péclet number Pe. Our prediction is constant for slow rotation (Ω < | N |) and scales as (Ω/ | N |) -1 for rapid rotation (Ω > | N |).

Figure 13 .

 13 Figure 13. The ratio 2 A / 2 c of magnetic energy to convective kinetic energy is shown alongside our predictions as a function of Ω/ | N | for a variety of MHD convection simulations. Filled shapes indicate simulations performed on spherical domains while open shapes indicate those performed in shellular domaints. All of these data were averaged over the simulation domain. Our prediction is constant for slow rotation (Ω < | N |) and scales as Ω/ | N | ≈ (Ω/ | N |) 2 for rapid rotation (Ω > | N |).

Figure 14 .

 14 Figure 14. The angle ξ between the pressure gradient and the entropy gradient is shown for a variety of simulations. Our prediction scales as (Ω/ | N | 0 ) 2 for slow rotation (Ω < | N | 0 ) and is constant for rapid rotation (Ω > | N | 0 ).

  obtain the Péclet number by equation (A17).

Figure C1 .

 C1 Figure C1. The relative latitudinal shear |∂ θ ln Ω | is shown as a function of Ω/ | N | 0 for observed convecting stars inferred from star spot measurements by Lurie et al. (2017) along with our predictions. The same data are shown with different colours representing different averaging schemes for | N | 0 which weight by different factors in radius and average different functions of | N | 0 .

  12 , Ω -0.41(Yadav et al.

See Appendix B of Paper I. MNRAS 000, 1-27 (2020)

MNRAS 000, 1-27 (2020)

The anelastic approximation is the limit as the Mach number and Froude number both go to zero, meaning that characteristic velocities are much less than both the sound speed and the freefall speed across the domain(Masmoudi 2007). This allows the pressure and gravity terms in the Navier-Stokes equation to be written in terms of an entropy perturbation away from a background state, and imposes the constraint that ∇ • ( ρv) = 0, where ρ is the background density

profile.3 These are domains which just exclude a small region in the centre of the sphere for numerical or algorithmic reasons.MNRAS 000, 1-27 (2020)

For details of how this was calculated for each simulation see Appendix A.

The greater uniformity of dimensionless parameters in these simulations is likely because the simulation must be fairly turbulent (advection-dominated) to see a dynamo in the first place, which effectively sets a floor on the Péclet number.

 6 See alsoDavidson (2013) for more discussion of these scaling relations.

See their equations (30) and (31).

Some authors (e.g.Augustson et al.

2012) report the root-mean square of the meridional circulation. Despite the similar name this is not the quantity we predict. By squaring before taking a time-average they are measuring the strength of axisymmetric convective motions rather than the slow mean meridional flow.MNRAS 000, 1-27 (2020)

ACKNOWLEDGEMENTS

The Flatiron Institute is supported by the Simons Foundation. ASJ thanks the Gordon and Betty Moore Foundation (Grant GBMF7392) and the National Science Foundation (Grant No. NSF PHY-1748958) for supporting this work. ASJ also acknowledges financial support from a UK Marshall Scholarship as well as from the IOA, ENS and CEBS to work at ENS Paris and CEBS in Mumbai. PL acknowledges travel support from the french PNPS (Programme National de Physique Stellaire) and from CEBS. CAT thanks Churchill College for his fellowship. SMC is grateful to the IOA for support and hospitality and thanks the Cambridge-Hamied exchange program for financial support. ASJ, SMC and PL thank Bhooshan Paradkar for productive conversations related to this work. The authors thank Jim Fuller, Douglas Gough, Steven Balbus, Chris Thompson for comments on this manuscript. We are also thankful to the anonymous referee, whose suggestions led us to include more data as well as the sensitivity analysis in Appendices C and D. ASJ is grateful to Lars Bildsten for helpful suggestions regarding the data analysis, to Frank Timmes and Matteo Cantiello for suggestions regarding the presentation, and to Daniel Lecoanet, Evan Anders and Matthew Browning for comments regarding the interpretation of simulations.

DATA AVAILABILITY

Details of how the data from various sources were obtained, standardized and processed are provided in Appendix A, and all scripts and files needed to reproduce this analysis are available on https://doi.org/10.5281/zenodo.3992227.

rate. We take this to equal

|R∇ ln Ω|d ln R, (A6) on the assumption that most of the differential rotation is developed in the convection zone. We thereby infer

Stellar models were constructed with MESA which matched the masses and radii in table 2 of Deheuvels et al. (2015). These were used to compute |N | avg and ∫ convection zone d ln R.

A2.2 Klion Observations

Data were taken from the text of Klion & Quataert (2017) for the red giant Kepler-56. The differential rotation was seismically inferred and reported in the form of a power-law

with β ≈ 1, so we take |R∇Ω|/Ω ≈ 1. The rotation period is provided in the text as 74 ± 3d. A solar-metallicity stellar model was constructed with MESA which matched the mass and radius reported for Kepler 56. This was used to compute |N | avg .

A2.3 Nielsen Observations

Data were taken from table 1 of Nielsen et al. (2017). The differential rotation was seismically inferred for several solarlike stars and reported as envelope and interior angular frequencies. We then approximate the shear as

where ∆ ln R was computed as

as we have done in equation (7). Uncertainties were propagated via equation (A2). Objects were matched by effective temperature to the grid of main-sequence stellar models from appendix A1.3, from which we computed |N | avg and Ω/|N | avg .

A3 Hydrodynamic Simulations

We describe how we extract the relevant data from various hydrodynamic simulations.

A3.1 Browning 2004 Simulations

These data come from three-dimensional anelastic hydrodynamic simulations of convection in a rotating spherical domain (Browning et al. 2004). The differential rotation is reported as the difference in angular velocity between the equator and a latitude of π/3, normalized to the equatorial angular velocity, so we write

Ω equator .

(A11) Browning et al. (2004) report the Reynolds number

where ˜ is the root-mean square of the convective velocity, L is the radial extent of the convective region and ν is the microscopic viscosity. They also report the Prandtl number

and the Rayleigh number

where α is the thermal diffusivity, ρ is the density, s is the entropy, ∆s is the entropy change across the domain and g is the acceleration owing to gravity. We assume that the partial derivative is taken at constant pressure,

With this information we obtain the convective velocity by

and the Péclet number as

The Rossby number actually realised in the flow is also reported. This they define as

where ω φ,convective ≈ c /h is the root-mean square of the vorticity of the convective flow. So we find

where the term in parentheses is given by equation (A16).

A3.2 Augustson 2012 Simulations

These data come from three-dimensional anelastic hydrodynamic simulations of convection in a rotating spherical shell domain (Augustson et al. 2012). The differential rotation is reported as the difference in angular velocity between the equator and a latitude of π/3, normalized to the equatorial angular velocity, so we process this as equation (A11). The Rossby number actually realised in the flow is also reported, as equation (A18), so we handle this as before.

The temperature difference ∆T between the equator and a latitude of π/3, averaged over depth, is reported. To obtain the baroclinic angle ξ from ∆T we write ) This is valid in the limit of a radial pressure gradient and a mostly-radial entropy gradient, so that ξ is the small projection of ∇s/|∇s| along e θ . Augustson et al. (2012) define

where Ra is the Rayleigh number, ∆S is the radial change

A3.5 Brun 2009 Simulations

These data come from three-dimensional anelastic hydrodynamic simulations of convection in a rotating spherical domain reported by Brun & Palacios (2009). The differential rotation is reported as the mass-weighted average of the squared difference between the rotation rate and the mean rotation rate. We translate this into a shear as we did in appendix A3.4. The Rossby number actually realised in the flow is reported in the form of equation (A18) and we compute Ω/|N | from this as in appendix A3.1. We compute the ratio c /h|N | in the same way as we did for the data from Augustson et al. (2012), using the Rossby, Prandtl, Taylor and Rayleigh numbers. Two differences are that Brun & Palacios (2009) used a Prandtl number of 1 and that they report the Rayleigh number as in equation (A14). So we must divide Ra by (γ -1)/γ before using equation (A34). They also do not report γ so we assume that it is 5/3. As before, we obtain the Péclet number from equation (A17).

A3.6 Brun 2017 Simulations

These data come from three-dimensional anelastic hydrodynamic simulations of convection in a rotating spherical domain (Brun et al. 2017). The differential rotation is reported as the mass-weighted average of the squared difference between the rotation rate and the mean rotation rate. We translate this into a shear as we did in appendix A3.4. The Rossby number actually realised in the flow is reported in the form of equation (A18) and we compute Ω/|N | from this as in appendix A3.1.

The temperature difference between the equator and a latitude of π/3 is also reported. We follow the same procedure as in appendix A3.2 to compute ξ from this, except that we compute g as in appendix A3.3.

We compute the ratio c /h|N | in the same way as we did in appendix A3.5. As before, we obtain the Péclet number from equation (A17).

A3.7 Featherstone 2015

These data come from three-dimensional anelastic hydrodynamic simulations of convection in rotating spherical shells, reported by Featherstone & Miesch (2015). The differential rotation is reported as the mass-weighted average of the squared difference between the rotation rate and the mean rotation rate. We translate this into a shear as we did in appendix A3.4. The Rossby number actually realised in the flow is reported in the form of equation (A18), and we compute Ω/|N | from this as in appendix A3.1. Featherstone & Miesch (2015) report a modified Rayleigh number

where these quantities are all as defined previously. Using equations ( 2), ( 1) and (A22) we find

from which we compute Ω/|N |.

To compute c /h|N | we use the Rossby number actually realised in the flow. With this we find

We compute the Reynolds number as Re = 2ΩLRo/ν from the given Ω, L and viscosity ν. We then compute the Péclet number from equation (A17).

A3.8 Augustson 2013 & 2016 Simulations

These data come from anelastic three-dimensional anelastic hydrodynamical simulations of rotating convection in a spherical domain (Augustson et al. 2016) and spherical shells (Augustson et al. 2013). The differential rotation is reported as the mass-weighted average of the squared difference between the rotation rate and the mean rotation rate. We translate this into a shear as we did in appendix A3.4. The Rossby number actually realised in the flow is reported in the form of equation (A18) and we compute Ω/|N | from this as in appendix A3.1.

The temperature difference between the equator and a latitude of π/3 is also reported by Augustson et al. (2013). We follow the same procedure as in appendix A3.2 to compute ξ from this. We likewise compute the ratio c /h|N | in the same way as we did in appendix A3.2. We obtain the Péclet number from equation (A17).

A3.9 Guerrero Simulations

These data come from simulations of three-dimensional anelastic convection in rotating spherical shells, reported in table 1 of Guerrero et al. (2013). The reported differential rotation rate is the surface equator-pole difference in Ω normalised to that of the equator. Because this spans π/2 radians we divide by π/2 to give it the same units as ∂ θ Ω/Ω. We take this as to be approximately |R∇Ω|/Ω.

They report R * as in equation (A40). We compute Ω/|N | from this as in appendix A3.7.

They also report the Rossby number in the form

where ṽ is the root-mean square of the fluctuating velocity, L is the vertical size of the domain and Ω is the rotation rate of the frame with zero total angular momentum. We use this to find

We were unable to determine the Péclet number in these simulations because the EULAG software instrument they employ does not use explicit diffusivities and so makes it difficult to determine precise values for dimensionless parameters such as Pe.

A3.10 Mabuchi Simulations (HD)

These data come from three-dimensional fully compressible hydrodynamic simulations of rotating convection in spherical shells, reported by Mabuchi et al. (2015). Data were taken from their table 1. The reported differential rotation is in the form

where Ω is evaluated at the surface of the star on the equator.

The convective Rossby number is reported as

where

d is the vertical size of the domain and

With these definitions, the convective Rossby number simplifies to

The Rossby number actually realised in the flow is also reported as

where d is the radial extent of the convection zone. From this we obtain

We obtain the Péclet number from equation (A17). We calculate the Reynolds number using the height of the unstable layer (0.3), the thermal diffusivity (3.9×10 -4 ) and the reported root-mean-square velocities, all in code units.

A3.11 Käpylä Simulations

These data come from simulations of three-dimensional fully compressible hydrodynamical convection in rotating spherical shells (Käpylä et al. 2011). Data were taken from their table 1. The reported differential rotation rate is the surface difference of Ω over 2π/3 radians normalised to that of the equator. We thus divide by 2π/3 to estimate ∂ θ Ω/Ω and take this as to be approximately |R∇Ω|/Ω.

The Rayleigh number is reported as in equation (A48). The Coriolis number was reported as Ro -1 , with Ro defined by equation (A50). The Prandtl number is reported as well. Finally, the Reynolds number is reported as

where d is the vertical size of the domain and ν is the viscosity.

With these pieces we can write the Taylor number, defined by equation (A47), as

and find

from which we compute Ω/|N |. We further obtain

We calculate the Péclet number from equation (A17).

A3.12 Gastine 2012 Simulations

These data come from three-dimensional hydrodynamic simulations of anelastic convection in rotating spherical shells, reported by Gastine & Wicht (2012).

The polytropic index is reported as m = 2, so that γ = 3/2. The same Ekman number was used in each of these simulations, and reported as

where Ω frame is the angular velocity of the simulation reference frame, d is the vertical size of the domain and ν is the microscopic viscosity.

The Prandtl number was also held constant and reported as

where α is the microscopic thermal diffusivity. Finally, the aspect ratio η was held constant and reported as

This is related to the vertical size d of the domain by

The Rayleigh number is reported as

where c p is the specific heat at constant pressure, g is the acceleration due to gravity and ∆S is the difference in specific dimensionful entropy across the domain. The ratio ∆S/c p may be written as

where s is the dimensionless entropy defined in equation ( 2). So

Approximating the entropy gradient as radial, we write

and, with these pieces, we may write the ratio

The ratio of the axisymmetric azimuthal kinetic energy in the frame rotating at Ω frame to the total kinetic energy Convective Differential Rotation II 23 was reported. In the language of Brown et al. (2008) this ratio is

from which we can compute DRKE/CKE. The Rossby number was also reported as

where we have inferred the definition of Re from the text surrounding their equation ( 26) and 2 1/2 is the time-and volume-averaged root-mean square of the velocity. Putting these two numbers together we find

We computed the ratio of the convection speed to h|N | as

where the final square root factor corrects the Rossby number for the fraction of energy actually in convection, as opposed to differential rotation, and

is the number of pressure scale heights in the domain, which we cut off at unity because that is the relevant mixing length in the limit of no density stratification. Here N ρ is the reported number of density scale heights in the domain. We obtain the Péclet number using the relation

where Ek is the Eckman number, and use this Reynolds number in equation (A17).

A3.13 Gastine 2013 Simulations

These data are from three-dimensional hydrodynamic simulations of anelastic convection in rotating spherical shells, reported by Gastine et al. (2013). The Prandtl number, aspect ratio, domain size, poltropic indices are the same as those in appendix A3.12. The Rayleigh number and Eckman number are reported in the same manner but vary from simulation to simulation. We compute the Rossby number from these in the same way as in appendix A3.12. Gastine et al. (2013) report the time-averaged Rossby number at the surface on the equator (Ro e ). We interpret this directly as a measure of differential rotation, so that

where Ω frame is defined as in appendix A3.12.

We compute the ratio of the convection speed to h|N | with equation (A64) to be

where n is defined by equation (A72) and Re is the Reynolds number of fluctuations in the flow.

We obtain the Péclet number from equation (A17). Here we use the fluctuating Reynolds number, because the other Reynolds numbers reported by Gastine & Wicht (2012) are broken down into specific spatial components.

A3.14 Aurnou Simulations

These data are from simulations of three-dimensional hydrodynamic Boussinesq convection in rotating spherical shells (Aurnou et al. 2007). They appear in the second summary figure of Gastine et al. (2014) and were extracted with automated graphic data extraction software to provide α as in equation (A45). Points with α < 0.02 were removed from the sample because the plot was scaled linearly and below this limit the figure resolution was insufficient to ensure accuracy. To measure the rotation we note that α is provided as a function of the convective Rossby number, which Gastine et al. (2014) define

Inspection of equations (2-4) of (Aurnou et al. 2007) reveals that this quantity is just |N |/Ω, from which we compute Ω/|N |.

We obtain the Péclet number from equations (A17) and (A73). To find the Ekman number we note that the quantity Ra * reported by Gastine et al. (2014) is not just the convective Rossby number, but rather

From this we can solve for Ek and thence for Pe.

A3.15 Kaspi Simulations

These data are from three-dimensional anelastic general circulation simulations over a thick spherical shell (Kaspi et al. 2009). They appear in the second summary figure of Gastine et al. ( 2014) and were extracted with automated graphic data extraction software providing the quantity α given in equation (A45). These were extracted and processed as in appendix A3.14.

A3.16 Gilman Simulations

These data come from simulations of three-dimensional Boussinesq convection in rotating spherical shells reported in Gilman (1977) and Gilman (1979). The results appear in the second summary figure in Gastine et al. (2014). These were extracted and processed as in appendix A3.14.

A3.17 Rogers Simulations

These data are from three-dimensional hydrodynamic simulations of anelastic convection in a spherical domain reported by Rogers (2015). The results were taken from their table.

The table lists c for the two extremal convective forcing functions and so we linearly interpolate these as c = [2.9 + (q -1.5) * (4.5 -2.9)] km s -1 , (A78)

where q equals their Q/ c and is the convective forcing. We take the scale height to be the size of their convection zone, which is 0.3R . Using the Ω values in their table we obtain c /hΩ. We then employ the scaling in table 2 to write

when Ω > |N | and

when Ω < |N |. Note that this definition assumes that our scaling for the convective velocity holds, and so we cannot use these data to test that scaling.

The differential rotation is reported as the ratio of the core-to-envelope angular velocity. We estimate the shear from this with equation (A9).

To obtain the Reynolds number we use their reported viscosity ν = 4 × 10 13 cm 2 s -1 and let

We then determine the Péclet number through equation (A17) using the Prandtl number

where the thermal conductivity α = 5 × 10 11 cm 2 s -1 .

A3.18 Raynaud Simulations

These data are from simulations of three-dimensional anelastic hydrodynamic convection in rotating spherical shells, reported in table B.1 of Raynaud et al. (2018). The table contains the number of density scale heights in the simulation domain (N ρ ), the ratio of the outer radius of the domain to the inner radius ( χ) and the Nusselt number (Nu) at both the equator and pole. The Nusselt number is proportional to the radial entropy gradient. The latitudinal derivative of this is then

If we estimate that the radial derivative in the numerator to produce a factor of order the pressure scale height then

By equation (A20) we then obtain

The polytropic index n is not specified, so we assume n = 3/2. The pressure scale height is then related to the density scale height by a factor of γ = 5/3, so

where h ρ is the density scale height. The average density scale height in the simulation domain is just the size of the domain divided by the number of scale heights it contains, so

Inserting this into equation (A86) we find

We then estimate

and with equation (A89) obtain ξ.

The table contains the Ekman number

The Rayleigh number is reported as

where we have used α = κc p to convert from their notation to ours. Using d/r = 1 -χ we see that this Rayleigh number differs from that of equation (A62) by a factor of γ(γ-1) -1 (1χ) -2 , so we may use equation (A65) to find

A4 MHD Simulations

We describe how we extract the relevant data from various MHD simulations.

A4.1 Soderlund Simulations

These data are from simulations of three-dimensional Boussinesq MHD convection in rotating spherical shells of varying thickness (Soderlund et al. 2013). These appear in the second summary figure in Gastine et al. (2014) and were extracted with automated graphic data extraction software, providing the quantity α given in equation (A45). Points with α < 0.02 were removed from the sample because the plot was scaled linearly and below this limit the figure resolution was insufficient to ensure accuracy.

To measure the rotation we note that α is provided as a function of the convective Rossby number, defined as |N |/Ω. We obtain the Péclet number from equations (A17) and (A73).

A4.2 Mabuchi Simulations (MHD)

These data are from three-dimensional fully compressible MHD simulations of rotating convection in spherical shells, reported by Mabuchi et al. (2015). Data were taken from their table 1. The reported differential rotation is in the form of α as defined by equation (A45). The data analysis was performed in the same way as for the hydrodynamic case in appendix A3.10. Mabuchi et al. (2015) also report the magnetic energy and the kinetic energy in the meridional plane and we take the ratio of these to be 2 A / 2 c . We obtain the Péclet number from equation (A17). We calculate the Reynolds number using the height of the unstable layer (0.3), the thermal diffusivity (3.9 × 10 -4 ), and the reported root-mean-square velocities, all in code units.

A4.3 Brun 2005 Simulations

These data are from three-dimensional anelastic MHD simulations of convection in a rotating spherical domain, reported by Brun et al. (2005). The differential rotation is reported as the mass-weighted average of the squared difference between the rotation rate and the mean rotation rate. We translate this into a shear as we did for the data of Brown et al. (2008), taking into account the scaling of c in MHD systems given in Table 1.

The Rossby number actually realised in the flow is also reported. This they define as

For this calculation we take h ≈ R because the average depth is comparable to the radius. We then calculate Ω/|N | from this as in appendix A3.1. The magnetic energy and the convective kinetic energy are also reported. We use this to calculate 2 A / 2 c . We compute the ratio c /h|N | in the same way as we did in appendix A3.1. The one difference is that Brun et al. (2005) do not report the adiabatic exponent γ, so we assume it is 5/3. Brown et al. (2008) report the Péclet number so we do not need to compute it.

A4.4 Augustson 2013 & 2016 Simulations

These data come from anelastic three-dimensional MHD simulations of rotating convection in a spherical domain (Augustson et al. 2013) and spherical shells (Augustson et al. 2016). The differential rotation is reported as the mass-weighted average of the squared difference between the rotation rate and the mean rotation rate. We translate this into a shear as we did in appendix A3.4, taking into account the scaling of c in MHD systems given in Table 1. The data analysis was performed in the same way as for the hydrodynamic case. The magnetic energy and the convective kinetic energy are also reported. We use this to calculate 2 A / 2 c . We obtain the Péclet number from equation (A17).

A4.5 Varela 2016 Simulations

These data are from anelastic three-dimensional MHD simulations of rotating convection in a spherical shell domain, reported in tables 1 and 2 of Varela et al. (2016). The differential rotation is reported as the difference in angular velocity between the equator and a latitude of π/3, normalized to the equatorial angular velocity, so we follow equation (A11) to compute |R∇ ln Ω|.

The Rossby number actually realised in the flow is reported in the form of equation (A18). We calculate Ω/|N | from this as in appendix A3.1. The magnetic energy and the convective kinetic energy are also reported and used to calculate 2 A / 2 c . We compute the ratio c /h|N | in the same way as we did in appendix A3.2, using the Rossby, Prandtl, Taylor and Rayleigh numbers. The one difference is that the Rayleigh number is reported as in equation ( A14), so we must divide it by (γ -1)/γ before using equation (A34). We obtain the Péclet number from equation (A17).

A4.6 Yadav 2013 Simulations

These data are from Boussinesq three-dimensional MHD simulations of rotation convection in a spherical shell domain, reported in table A.2 of Yadav et al. (2013a). The Rayleigh number is reported as Finally, the Lorentz number is reported as ) where ... denotes a volume average and d is the domain thickness. From this we obtain

A5 Solar Data

A helioseismic inversion of the solar profile was obtained from Antia (2016). It corresponds to that appearing in Antia et al. (2008). This includes Ω as a function of position throughout the convective envelope as well as in the upper portions of the radiative envelope. This profile was then supersampled on to a grid running from 0.5R to R in the radial direction, with 100 uniformly spaced points, and from 0 to π in the latitude, with 70 uniformly spaced points. By applying a differentiating Gaussian filter with width equal to five grid points in each of the radial and latitudinal directions we computed R∇Ω/Ω from this profile. We then averaged each component of this over latitude, weighted by sin θ, and took the square root of the result to produce a measure of the mean radial and latitudinal differential rotation at each radial slice.

The radial profile of |N | was obtained from a MESA model of a 1 M star with metallicity Z = 0.02 at an age of 4.6 Gyr. This was used to compute Ω/|N | everywhere in the solar convection zone. We then averaged Ω 2 /|N | 2 over latitudes, weighting by sin θ, and computed the square root to determine the mean Ω/|N | in each radial slice. Altrock & Canfield (1972) provide an upper bound on the temperature difference between the solar pole and equator. Following equation (A20) we write ) We then note that, from equation (1),

and inserting equation (A27) we find

Because the temperature measurement is at the surface we compute ξ using |N | 2 , g and r at the surface, estimating

The factor of 2/π accounts for the number of radians involved in estimating the latitudinal derivative from a finite difference. To plot this we take the rotation rate of the Sun to be its mean rate and compute Ω/|N | for the surface |N |.

A6 Jupiter Data

The profile of differential rotation in the surface layers of Jupiter was taken to be that given by Kaspi et al. (2018). This yields the variation of the characteristic zonal flow rate with depth but does not provide the shear itself. To compute this we first note that roughly 30 per-cent of Jupiter's surface is covered in latitudinal bands with velocities of order 100 km s -1 , while the remaining latitudes contain bands with velocities of order 25km s -1 (Kaspi et al. 2018). The former are also approximately twice as wide per band as the latter. There are approximately 15 bands in total, and so we estimate the typical shear at the surface to be

where R J ≈ 7×10 9 cm is the radius of Jupiter. We then assume that the band velocities scale with depth following the profile

(1 -α)e (r-α)/H + α

found by Kaspi et al. (2018), where α, H, a and ∆H are given therein. The shear then scales as

We further assume that the mean rotation period remains 9.92hr (Kaspi et al. 2018).

For regions deeper than 3000 km we bound the band velocity above by the inferred 6 m s -1 (Guillot et al. 2018) and use the band structure above to compute the shear. Furthermore we treat the number of bands as varying linearly in radial coordinate between 1 near the centre of the planet and 15 at the surface. This is an approximation of the cylindrical nature that Kaspi et al. (2018) assumed for Jupiter's differential rotation. Following this assumption we consider this estimate to be for both the spherical-radial and latitudinal shears.

All that remains is to compute |N | 0 . To do this we ran a MESA model of a Jupiter-mass planet irradiated by 12500 erg cm -2 s -1 at a pressure of 1 bar. This produced a model within 3 per cent of Jupiter's radius and an effective temperature of 136 K, close to the Jovian 1 bar temperature of 165 K [START_REF] Ingersoll | Dynamics of Jupiter's atmosphere[END_REF]). To ensure efficient convection we exclude the photosphere and the top scale height of the convection zone from our analysis.