Adam S Jermyn

Shashikumar M Chitre

Pierre Lesaffre

Christopher A Tout

Convective Differential Rotation in Stars and Planets I: Theory

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

We derive the scaling of differential rotation in both slowly-and rapidly-rotating convection zones using order of magnitude methods. Our calculations apply across stars and fluid planets and all rotation rates, as well as to both magnetized and purely hydrodynamic systems. We find shear |R∇Ω| of order the angular frequency Ω for slowly-rotating systems with Ω |N |, where N is the Brünt-Väisälä frequency, and find that it declines as a power-law in Ω for rapidly-rotating systems with Ω |N |. We further calculate the meridional circulation rate and baroclinicity and examine the magnetic field strength in the rapidly rotating limit. Our results are in general agreement with simulations and observations and we perform a detailed comparison with those in a companion paper.

INTRODUCTION

Differential rotation is one of the key complications in the study of stars. It involves breaking symmetries, structure formation and both heat and momentum transport. Importantly the origin of magnetic fields (Miesch & Toomre 2009), the transport of angular momentum (Cantiello et al. 2014) and the transport of material (Chaboyer & Zahn 1992) are all critically influenced by the scale and geometry of differential rotation. Moreover, differential rotation plays a significant role in setting the spins of stellar cores, which are of significant interest for understanding the spins of white dwarfs, neutron stars and recently black holes (Fuller & Ma 2019).

Over the past several decades helioseismology has permitted studies of the rotation profile of the solar convection zone (Christensen-Dalsgaard & Schou 1988). With the passage of time the data have become more precise and detailed, providing information on the time-variability of the rotation profile (Antia & Basu 2001;Thompson et al. 2003) as well as that of its gradients (Antia et al. 2008). Similarly, related quantities such as the meridional circulation (Rajaguru & Antia 2015) have now been characterized. The overall picture that has emerged for the Sun is one in which the solar differential rotation is of order |∇Ω| ≈ Ω/R, where Ω is E-mail: adamjermyn@gmail.com the angular frequency and R is the distance from the spin axis. This differential rotation reflects velocities which are large relative to the meridional circulation yet, depending on depth, may be either large or small relative to the convective velocity. These observations present a challenge: what phenomenon sets the scale of differential rotation in the solar convection zone?

Complementing the depth of solar observations, asteroseismic observations have begun to produce information on the rotation profiles of other convecting stars (Beck et al. 2012). Limits on the rotation profiles for several red giants and Sun-like stars are now known, and remarkably tell a similar picture of shear comparable to the rotation rate (Schunker et al. 2016;Klion & Quataert 2017;Benomar et al. 2018). This is in line with results from studies on a wider range of red giants. These do not attempt to localize the difference rotation but find |∇Ω| typically of order Ω/R (Deheuvels et al. 2015). While the precision of these data and analyses continues to improve an outstanding question is whether the bulk of the differential rotation in red giants lies in their convective envelopes or in their radiative interiors (Cantiello et al. 2014). This highlights the importance of determining the magnitude of differential rotation in convection zones.

Unfortunately, despite these observational successes, we still lack a fully explanatory theory of differential rotation. Early theories suggested solid body rotation (Stewartson 1966). As turbulence and the microscopic viscosity together serve to dissipate energy, fluid bodies without any forcing are expected to come to rigid-body rotational equilibrium. The characteristic time-scale for this is the rotation period because this is the only such scale for rotationally-driven kinetic turbulence. This is extremely rapid in the context of stars and planets and so should preclude differential rotation. That this is not observed is evidence of processes which inject energy into differential rotation. For instance, in convecting bodies, turbulence may be strongly anisotropic. Such anisotropy drives differential rotation and, if the turbulence is powered by heating rather than by the shear itself, it can maintain such a state indefinitely (Unno 1957;Kippenhahn 1963).

In addition to the expectation of solid-body rotation there was also the expectation of cylindrical rotation. This was due to the Taylor-Proudman theorem, which states that the rotation profile of an efficiently convecting (e.g. isentropic) region ought at least to have translational symmetry along the rotation axis (Hough 1897). This arises from a balance between the Coriolis and centrifugal effects, and so is very different in origin from the solid body expectation. Despite this clean result, observations indicate that the Sun does not obey any such constraint (Di [START_REF] Mauro | Structure and Dynamics of the Interior of the Sun and Sun-like Stars[END_REF][START_REF] Gough | The inversion problem[END_REF]. This could be due to a variety of effects which the theorem neglects, including viscosity and turbulent stresses and MHD phenomena. Additionally, convection zones are not perfectly isentropic and this leads to the so-called thermal wind correction to the Taylor-Proudman state. One aim in this work is to determine which effects serve to break the Taylor-Proudman state and under what circumstances.

More recently, and in part owing to dramatic improvements in observational capabilities, differential rotation has garnered substantial theoretical attention. Some authors argue that thermal wind balance and entropy gradients dominate the solar rotation profile (Miesch et al. 2006;Balbus & Schaan 2012), while others have called this into question (Brun et al. 2010;Brun & Toomre 2002). Observations suggest that the thermal wind term is substantial, though there remain uncertainties as to its precise contribution (Caccin et al. 1976;Rast et al. 2008;Teplitskaya et al. 2015). Other models suggest that turbulent anisotropy is the most relevant factor in setting the differential rotation [START_REF] Ruediger | Differential rotation and stellar convection. Sun and the solar stars Sackur O[END_REF]Kueker et al. 1993;[START_REF] Kitchatinov | IAU Symposium[END_REF], and more complex models with various parameterizations have also been proposed (Tuominen & Ruediger 1989;Kissin & Thompson 2015;Brun & Rempel 2009). One of our goals in this work is to understand which of the proposed effects matter the most and under what circumstances.

Numerical investigations of these issues have proven more successful at reproducing details of the solar rotation profile (Thompson et al. 2003;Miesch & Toomre 2009), though typically not for configurations that reflect realistic convective velocities or luminosities. These differences may reflect uncertainties in sub-grid physics, or could be due to the fact that physically realistic resolution and diffusivites remain out of reach (Miesch 2005;[START_REF] Käpylä | [END_REF]). Nevertheless, the results which have been found in this way are intriguing. For instance, red giants are found in both simulation and asteroseismic inference to exhibit significant differential rotation, including cases where the angular velocity changes sign (Brun & Palacios 2009).

In this work we aim to understand the magnitude of differential rotation in the convection zones of stars and gaseous planets, and to understand how these solutions connect to the Keplerian limit of an accretion disk. That is, we aim to determine the approximate magnitude and scaling of |∇Ω| in these systems.

It is important to emphasise that our arguments are purely from an order-of-magnitude scaling perspective. In particular, we generally assume that dimensionless geometric factors are of order unity rather than being very large or small. We believe that this is likely in most cases because the alternative is significant coincidence in the geometries of various fields which are determined by a variety of fundamentally dissimilar physical processes. Thus, for instance, we are agnostic on whether baroclinic pumping or turbulent stresses play a greater role in the limit of slow rotation (c.f. Miesch et al. 2006) because we find that they exhibit identical scaling and are related by a dimensionless factor of order unity.

We further caution that there remain significant uncertainties in the precise scaling of turbulent velocities and stresses. We have done our best to estimate these scalings from a combination of symmetry arguments and mixing length approaches. The agreement we find with observations and simulations in the companion manuscript suggests that these tools have been useful, but they are inherently simplifications and we feel compelled to point out their limitations.

We begin in Section 2 with a discussion of our assumptions. In Section 3 we examine the vorticity equation and derive the form which we use in all subsequent analysis. We then consider in turn magnetic fields (Section 4), the condition of thermal equilibrium (Section 5) and the thermal wind contribution (Section 6), deriving helpful expressions relating the magnitudes of different effects. In Section 7 we introduce an asymptotic scaling approach which helps to organise the remainder of the calculation.

The remaining Sections focus on our key results. In Sections 8, 9 and 10 we examine the differential rotation in both rapidly and slowly rotating convection zones in both the MHD and non-magnetic limits. In Appendices F and G we show how these results relate to the inverse cascade in density-stratified systems and how the limit of rapid rotation transitions into a Keplerian state. These are separated from the main text because their results are applicable in somewhat more limited circumstances.

Taken together our results provide a unified theory of the magnitude of convective differential rotation which covers the full range from accretion discs to gas planets to the most massive stars. We provide a detailed comparison of our model with both observations and simulations in a companion manuscript, so we conclude by discussing the limitations of our analysis (Section 12), summarizing our results (Section 13, Table 4) and commenting on their astrophysical implications (Section 14).

ASSUMPTIONS

For simplicity, we make 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.

Small violations of these assumptions do not undermine our conclusions. For instance so long as the entropy is logarithmic in pressure and density, and the sound speed is of order P/ρ, corrections owing to radiation pressure and degeneracy are not a problem.

Similarly, the assumption of axisymmetry is meant not in each instant but rather in a time-averaged sense. This is a much weaker condition and importantly does not run afoul of Cowling's anti-dynamo theorem (Cowling 1933;Parker 1955). So for example turbulence and dynamo cycles may produce temporary deviations from axisymmetry but the long-term average behaviour must be axisymmetric.

The most important assumption lies in our treatment of geometric factors. We do not treat the effects of spherical geometry in detail: we approximate latitude-dependent effects with their averages over θ, and take angular derivatives to produce factors of r -1 . Moreover, aside from considering the ratio of the pressure scale height h to the radius r we do not consider the effects of the aspect ratio or depth of the convection zone, and do not differentiate between spherical or shellular geometries. On the other hand we are highly concerned with the consequences of spherical symmetry and the ways in which it breaks. For instance a non-rotating selfgravitating system with no external or fossil magnetic field is spherically symmetric and so, even though the convective stresses do not vanish, the angular momentum they transport does. Hence in the slowly-rotating regime there is a small, rotation-dependent geometric factor associated with rotation breaking this symmetry which relates the stress that transports angular momentum to that which does not. We pay significant attention to such terms. More broadly, we effectively assume that, apart from any symmetries, the boundary conditions and geometry of the system are relatively generic. A consequence of this is that we generally avoid assuming that terms in the solutions are tuned to be irrelevant or that they contain geometric factors which scale with rotation rate.

VORTICITY

First, we derive the laws governing the evolution of angular momentum in axisymmetric systems, paying particular attention to the distinction between meridional and longitudinal components. We then argue that these systems are likely to be in a quasi-steady state. We begin by defining vorticity of the fluid as

r θ z R Ω ∂ z Ω ∂ R Ω u θ u r
ω ≡ ∇ × , (1)
where is the velocity. The vorticity is the angular velocity of the fluid about a point, and so is closely related to the differential rotation. In particular in the limit where rotation dominates = RΩe φ and

ω = Ωe z + R∇ × (Ωe φ ), ( 2
)
where Ω is the local angular velocity about the z axis and R is the cylindrical radial coordinate (Fig. 1). We do not impose this limit but it is worth keeping in mind because it gives an intuitive connection between vorticity and the force of rotation. By taking the curl of the Navier-Stokes equation we have that, in the absence of external forcing (e.g. tides), the vorticity evolves according to

∂ω ∂t = ω•∇ -ω∇• -•∇ω+ 1 ρ 2 ∇ρ×∇P+∇× 1 ρ ∇ • T +∇× F B ρ , (3)
where F B is the magnetic force, ρ is the density, p is the pressure and T is the turbulent fluid stress such that T i j is the flux of e i momentum along the direction e j , such that

(∇ • T) i = j ∂T i j ∂ x j . (4)
We are interested in axisymmetric systems, where there is a natural distinction between the meridional and longitudinal components of the flow (Fig. 1). To make this explicit we write the meridional flow as

u(R, z) ≡ (R, z) -Ω(R, z)Re φ .
(5)

Inserting this into equation (3), we then obtain (Appendix A)

∂ω ∂t = ω • ∇u -ω∇ • u -u • ∇ω + e φ (ω • ∇)(ΩR) -e φ Ωω R + 1 ρ 2 ∇ρ × ∇P + ∇ × 1 ρ ∇ • T + ∇ × F B ρ . (6)
The first line of this equation describes kinematic effects associated with the rotation and circulation while the second describes the effects of the thermal wind, turbulent stresses and magnetic stresses. Astrophysical systems are likely to be near angular momentum equilibrium because the time-scale over which shear turbulence transports momentum is quite short. To see this note that the diffusivity of shear turbulence over the whole system is of order

1 ν ≈ R 2 |∇ |, (7)
where R is the radius of the star. This is because the shear time-scale is set by |∇ | and its associated length-scale is set by the distance over which it persists. The diffusive timescale associated with a shear is therefore

τ diff ≈ R 2 ν ≈ |∇ | -1 , (8) which means that d|∇ | dt ≈ - |∇ | τ diff ≈ |∇ | 2 , (9)
so at times longer than the initial

τ diff |∇ | ∝ (|∇ |t) -1 . ( 10
)
What this implies is that transients ring down on a timescale comparable to their size. So after astrophysical timescales, transient shears are likely very small, and the system is well-approximated by instantaneous momentum equilibrium. This applies even in the presence of secular evolution owing to nuclear processes and wind losses and suggests that those effects are generally not enough to violate angular momentum equilibrium. There are two cases in which this argument fails. The first is in accretion discs, where the diffusivity is suppressed relative to equation (8), resulting in less relative angular momentum transport and hence longer equilibration timescales. However, in this case transients must ring-down on time-scales of order Ω -1 because non-Keplerian motion results in the fluid centrifugally adjusting its orbit on this time-scale. So, once more we find that transient effects in the differential rotation decay quickly.

The second case is in systems involving waves and turbulence. Equation ( 6) typically supports linear oscillatory motions, such as Alfvén waves and sound waves, as well as instabilities such as those associated with a dynamo [START_REF] Käpylä | [END_REF] and convection (Böhm-Vitense 1958). We take these motions to have been averaged over time, such that mean effects appear in the stress tensor T and so we neglect them wherever they appear explicitly in the vorticity equation. We also specifically neglect gravity waves because these are evanescent in convection zones (Fuller et al. 2014).

We now set the time derivative in equation ( 6) to zero and find

0 = ω • ∇u -ω∇ • u -u • ∇ω + e φ ω • ∇(ΩR) -e φ Ωω R + 1 ρ 2 ∇ρ × ∇P + ∇ × 1 ρ ∇ • T + ∇ × F B ρ . (11)
Furthermore, in steady state conservation of mass requires

∇ • (ρu) = 0 (12)
Inserting this into the second term of equation ( 11) and separating this vector-valued equation into its meridional and e φ components we find

0 = ω m • ∇u + ω m u • ∇ ln ρ + u • ∇ω m + ∇ × 1 ρ ∇ • T m + ∇ × F B ρ m (13)
and

0 = R -1 ω φ u R + ω φ u • ∇ ln ρ + u • ∇ω φ + ω • ∇(ΩR) -Ωω R + ρ -2 e φ • ∇P × ∇ρ + e φ • ∇ × 1 ρ ∇ • T + e φ • ∇ × F B ρ . ( 14
)
where we denote the projection of a vector into the meridional plane by the subscript m, i.e.

u ≡ m . ( 15
)

MAGNETIC FIELDS

The appearance of the magnetic force in equation ( 11) means that to close this equation we must address the source of magnetism. There are typically three sources of magnetic fields in astrophysical contexts, fossil fields, turbulent dynamos and externally imposed fields. A fossil field is present at the time of formation of the body. From the perspective of a scaling analysis there is little difference between a fossil field and an externally imposed field: both may break symmetries, including axisymmetry, and in both cases there is no guarantee of a relationship between the magnitude of the field and other properties of the body. Hence, they may be analyzed together.

Such an analysis awaits future study. For now, we consider the many cases in which the dominant magnetic field is generated by turbulent dynamo processes. In such cases the field obeys the same symmetry, on average, as the turbulence which drives it 2 , so its contribution vanishes in highly symmetric situations just as does the contribution of the turbulence. What remains then is due to the fluctuations, which we absorb into the turbulent stress, such that

0 = ω m • ∇u + ω m u • ∇ ln ρ + u • ∇ω m + ∇ × 1 ρ ∇ • T m ( 16
)
and

0 = R -1 ω φ u R + ω φ u • ∇ ln ρ + u • ∇ω φ + ω • ∇(ΩR) -Ωω R + ρ -2 e φ • ∇P × ∇ρ + e φ • ∇ × 1 ρ ∇ • T , (17)
where T now includes the turbulent magnetic stress

T B ≡ 1 4π B ⊗ B - 1 2 IB 2 . ( 18
)
So what we must determine is how this compares with the characteristic scale of the fluid stress. There are two regimes of interest, namely that of slow rotation and that of rapid rotation, which we discuss in sections 8 and 10.

THERMAL EQUILIBRIUM

Equations ( 16) and ( 17) involve the meridional circulation and the baroclinicity, so we must close the system of equations with a study of heat transport. We do so through the equation of thermal equilibrium, which reads

µ -1 ρc p T u • ∇s + ∇ • F = 0, ( 19
)
where µ is the mean molecular weight, ρ is the density, c p is the specific heat at constant pressure, T is the temperature,

s ≡ 1 γ -1 (ln P -γ ln ρ) (20)
is a dimensionless entropy for an ideal gas (see appendix B) and

F = µ -1 ρc p TQ • ∇s (21)
is the convective flux with diffusivity tensor Q. For an ideal gas P ∝ µ -1 ρc p T so equation ( 19) may be written as

Pu • ∇s + ∇ • (PQ • ∇s) = 0. (22)
In the non-rotating limit this possesses barotropic solutions with ∇s and ∇P both radial. Outside of this limit that is generally not true, both because ∇P is distorted by centrifugal forces [START_REF] Eddington | [END_REF]) and because Q becomes anisotropic [START_REF] Kitchatinov | IAU Symposium[END_REF].

Coordinate System

To proceed we define e p to be the unit vector along the pressure gradient and e q ≡ e φ × e p to be a unit vector perpendicular to e p in the meridional plane (Fig. 2, upper-right).

breaks one of the symmetries of the system in a time-averaged sense. Hence it is not sufficient that the solar magnetic field at any moment points in a particular direction because when averaged over the solar cycle the field vanishes. As far as we are aware, no such spontaneous symmetry breaking has been observed in stars or seen in MHD simulations, so we assume that this does not occur. (lower-right) a magnetic field with Alfvén velocity A ; (upperright) surfaces of constant pressure and density meet at an angle of approximately λ, the unit vector e p points along the pressure gradient and the unit vector e q is perpendicular to it.

Meridional

We denote components of vectors by the subscripts p and q to mean the components along these unit vectors. It is also useful to define the baroclinicity

ξ ≡ e φ • (∇ ln P × ∇s) |∇ ln P||∇s| , (23)
which measures the extent to which the pressure and entropy gradients are misaligned (Fig. 2, lower-left). When ∇P and ∇s are nearly aligned, ξ measures the small angle between them. When they are further misaligned it approaches ±1.

For convenience, we also define

ξ ≡ ± 1 -ξ 2 , ( 24
)
where the branch of the square root is taken such that

∇s = |∇s| e p ξ + e q ξ . (25)
With this definition u = e p u p + e q u q .

(26)

We can relate ξ to the pressure and density gradients by defining

λ ≡ |∇ ln P × ∇ ln ρ| |∇ ln P||∇ ln ρ| , (27)
which is directly proportional to the thermal wind term in equation ( 17) (Fig. 2, upper-right). When ∇P and ∇ρ are nearly aligned, λ measures the small angle between them.

When they are orthogonal λ approaches unity. In analogy with ξ we also define

λ ≡ ± 1 -λ 2 , ( 28
)
where the sign is chosen to match that of ξ.

Meridional Circulation

We next consider the meridional circulation. This operates under two constraints: conservation of mass (equation 12) and thermal equilibrium (equation 22). In Appendix C2 we expand the first of these in the coordinate system we have just defined and obtain (equation C28)

|u p | ≈ λ + h r |u q |, ( 29
)
where h is the pressure scale height defined in equation (C12) as

h ≡ |∇ ln P| -1 = P ρg ( 30
)
and g is the acceleration owing to gravity. Equation ( 29) says that the vertical flow u p is suppressed relative to the horizontal flow u q by a factor involving the density stratification. In Appendix C3 we combine this with the condition of thermal equilibrium and find (equation C33)

ξ + h r u + |∇ • (PQ • ∇s)| P|∇s| ≈ 0, (31)
Equation ( 31) relates the magnitude of the meridional circulation to the conditions of thermal equilibrium, which in turn are related to the baroclinicity by the diffusivity tensor. This does not intrinsically mean that thermal equilibrium drives the circulation, just that thermal equilibrium demands that this equation be satisfied. Thus, for instance, if the vorticity equation is more sensitive to the meridional circulation than the momentum equation it could be that the momentum equation drives a circulation, in which case the causation runs from left-to-right and the circulation determines the baroclinicity (Miesch & Toomre 2009;Brun & Rempel 2009). Likewise, in the reverse case baroclinicity drives and determines the circulation [START_REF] Eddington | [END_REF]Osaki 1982;Maeder & Zahn 1998). In fact we show that both possibilities likely occur occur depending on the context. A related result from recent simulations is that, in the slowrotation limit, these two effects actually scale with rotation in a similar manner (Miesch & Toomre 2009). The momentum imbalance generally provides a stronger impetus. This may be why different approaches, focusing on one or the other, have made similar predictions for the circulation rate.

THERMAL WIND

We now turn to the thermal wind balance, which enters via the baroclinic term ρ -2 (∇P × ∇ρ) φ appearing in equations (17). Noting that ∇P and ∇ρ lie in the meridional plane we find

| ρ -2 (∇P × ∇ρ) φ | = P ρ |∇ ln P × ∇ ln ρ| (32) = γN 2 ξ ξ , (33)
See Appendix C4 for more detail on the algebra between equations ( 32) and ( 33). Here N is the Brünt-Väisälä frequency defined in equation (C2) as (Fig. 2, upper-left)

N 2 ≡ - γ -1 γ g • ∇s. ( 34
)
Because γ and ξ are of order unity we drop these and obtain 16) and azimuthal vorticity balance (17), and reflect their scaling in the relevant limits. The exponents α i,k,l are independent of the rotation rate. We solve these equations for ξ, u, R∂ R Ω and R∂ z Ω as functions of Ω. Note that because the meridional component of the vorticity equation has two components we have four equations and four unknowns, so the system is determined.

| ρ -2 (∇P × ∇ρ) φ | ≈ N 2 ξ, (35
The solution {φ l (Ω)} may be expanded as a Puiseux series (Aroca et al. 2008). In the limit of asymptotically large or small Ω these series are well approximated by power-laws. That is,

φ l ∝ Ω |N | β l . ( 37
)
Our aim is to determine the exponents β l . To do this we note that asymptotic solutions to equation ( 36) must have at least two of the terms in the equation be of comparable magnitude. Were this not the case there would be a single term which is asymptotically larger than all the others, which would preclude the sum of all terms vanishing. We therefore seek the β l which cause each equation to have two or more terms which are larger than the rest and scale in the same manner as one another.

It is important to note that there need not be a single unique solution. This may arise either because there is an unphysical branch or because either of multiple terms suffices to balance another. In the former case we employ physical arguments to eliminate the extraneous solutions. In the latter case we assume that solutions are not fine-tuned to select just one of the various terms which may balance, so that all such terms exhibit the same asymptotic behavior. So for instance in the slowly-rotating limit we shall find that either baroclinicity or a meridional circulation suffices to balance the equation of thermal equilibrium, and that these terms enter in the same manner in the vorticity equation. Without a reason to believe that the system fine-tunes to have the baroclinicity vanish and the meridional circulation carry the full burden, or vice-versa, we assume that they share the responsibility. Note that this means we often cannot say anything about the relative signs of the {φ l (Ω)} we obtain, because there could be more than two terms involved in the overall heat and momentum balance.

In each limit we compute the terms and present an heuristic argument for the scalings β l . To confirm our results we also perform a search over all possible pairs of terms in each equation which might balance and compute the required asymptotic scalings3 .

We discard solutions for which the chosen pair is not dominant. Because each φ vanishes as Ω → 0 by symmetry in the slowly rotating limit we discard any solution with any β l ≤ 0. Similarly, in the rapidly-rotating limit we discard any solution for which u/h|N | or |R∇Ω|/|N | diverge as Ω → ∞. This is because shears are dissipative, with energy loss per unit mass

E ≈ h 2 ω 3 , ( 38
)
where ω is the shear time-scale. The only source of energy in these systems is the entropy gradient, which may act either through baroclinic pumping or convective forcing. This has characteristic power scale

E ≈ h 2 |N | 3 , (39)
so any shear with ω > |N | dissipates faster than it is forced and cannot be sustained over long time-scales.

SLOW ROTATION

In the special case of a non-rotating body with no fossil field, every term in the vorticity equation vanishes. This follows from symmetry because the only preferred direction is radial and this implies that the system is spherically symmetric. All vector fields of interest must therefore be radial4 . Every term in the vorticity equation results from the curl of a vector field and the curl of a radial field with spherical symmetry vanishes. Consequently every term in equation ( 11) vanishes. This provides a useful starting point for perturbation theory in the slow-rotation limit. Note that here slow is with respect to the Brünt-Väisälä frequency |N |, so we assume that Ω |N |.

Magnetic Field (Ω |N |)

When the rotation is slow relative to the convective turnover time, the field approaches equipartition with the turbulent flow (Antia et al. 2000;Roberts & Glatzmaier 2000;Sreenivasan & Jones 2006) and the Alfvén speed (Fig. 2, lowerright)

A ≡ B √ 4πρ (40)
is comparable to the convection speed. This is what is typically found in simulations (Hotta et al. 2015;[START_REF] Augustson | Astrophysical Dynamics: From Stars to Galaxies[END_REF]) even up to rotation rates comparable to the turnover time [START_REF] Augustson | [END_REF]. So in this regime the magnetic contribution to the turbulent stress scales with the kinetic term and we can focus on the latter.

Stress (Ω |N |)

We turn now to the contributions from the convective stress.

When the rotation is slow the turbulence is primarily convective, with characteristic length scale h and characteristic time scale |N | -1 . As a result (Böhm-Vitense 1958)

T ≈ ρh 2 |N | 2 , ( 41
)
where T is the typical magnitude of entries in T.

In the non-rotating limit the stress is constrained by spherical symmetry to be of the form

T ≈ T rr 0 0 0 T θθ 0 0 0 T φφ , (42)
where we have written this tensor in spherical coordinates and the diagonal terms are of the same order of magnitude Gough (1978). In this limit we know that

∇ × 1 ρ ∇ • T Ω = 0. (43)
In the slowly-rotating limit then we can compute the contribution of the stress to the vorticity balance by only consider how rotation perturbs T. Using symmetry arguments about the way different perturbations behave under reflections we estimate these perturabtions to be (Appendix D2)

∇ × 1 ρ ∇ • T r ≈ λ + h r |N | (Ω + |R∇Ω|) 1 + ξ + u h|N | , (44)
∇ × 1 ρ ∇ • T θ ≈ |N | (Ω + |R∇Ω|) 1 + ξ + u h|N | (45)
and

∇ × 1 ρ ∇ • T φ ≈ |N | 2 ξ + u h + Ω 2 |N | 2 + Ω|R∇Ω| |N | 2 , (46)
Here we have assumed that the stress is analytic about Ω = 0 and that the leading order terms are always the lowest-order ones which are allowed by the symmetries of the problem.

Advective Terms (Ω |N |)

The advective terms in the vorticity equation are estimated in Appendix E following the approximations introduced in Appendix C.

Baroclinicity (Ω |N |)

To proceed further we must determine the baroclinicity. We do this by analyzing equation ( 22), accounting for the fact that the diffusivity tensor is no longer isotropic in the presence of rotation. We again employ symmetry arguments as well as the approximations of section 5. The result is equation (C52)

u + h |N | (|R∇Ω| + Ω) 2 + |N | 2 ξ ≈ 0, (47)
Term Magnitude Meridional ( 16)

ω m • ∇u u h (Ω + |R∇Ω |) ω m u • ∇ ln ρ u h (Ω + |R∇Ω |) λ + h r u • ∇ω m u h (Ω + |R∇Ω |) ∇ × 1 ρ ∇ • T m | N | (Ω + |R∇Ω |) 1 + ξ + u h | N | Azimuthal (17) R -1 ω φ u R u 2 h R ω φ u • ∇ ln ρ u 2 h 2 λ + h r u • ∇ω φ u 2 h 2 ω • ∇(ΩR) Ω |R ∇Ω| Ωω R Ω |R ∇Ω| ρ -2 (∇P × ∇ρ) φ ξ | N | 2 ∇ × 1 ρ ∇ • T φ | N | 2 ξ + u h | N | + Ω 2 | N | 2 + Ω| R∇Ω| | N | 2
Table 1. The magnitudes of the terms in equations ( 16) and ( 17) are summarised here. Factors of order unity have been dropped for simplicity. See appendix E for derivations of each of these.

which serves to relate the meridional circulation, differential rotation, and baroclinicity through the condition of thermal equilibrium.

Results (Ω |N |)

We now have enough information to estimate the magnitude of each term in the heat equation ( 19) and the meridional and azimuthal vorticity equations (16 and 17). Putting it all together and dropping sub-dominant terms these equations may be written in terms of the magnitudes of their terms as

0 = u h |N | + ξ |N | 2 + (Ω + |R∇Ω|) 2 (48) 0 = u h (Ω + |R∇Ω|) + |N | (Ω + |R∇Ω|) 1 + ξ + u h|N | (49)
and

0 = u 2 h 2 + Ω|R∇Ω| + |N | 2 ξ + Ω 2 . ( 50
)
Consider starting with Ω = 0 and gradually spinning the system up. The first equation indicates that the perturbation to the heat transport is of order Ω 2 . This may be balanced either by a meridional circulation scaling similarly, baroclinicity scaling similarly or rotational shear scaling as Ω. Suppose for the moment that all three of these are realized. Then 1

Ω/|N | ξ, |N | u/h, so the meridional component reduces to 0 = |N | (Ω + |R∇Ω|) , (51)
which is indeed solved asymptotically by |R∇Ω| ≈ Ω. Finally in the azimuthal vorticity equation with u/h|N | ≈ Ω 2 /|N | 2 we see that the meridional circulation drops out, leaving

0 = Ω|R∇Ω| + |N | 2 ξ + Ω 2 , ( 52
) which is consistent with ξ ≈ Ω 2 /|N | 2 and |R∇Ω| ≈ Ω.
An exhaustive search of asymptotic solutions as described in section 7 reveals three other consistent choices of scalings. Each of these has |R∇Ω| ≈ Ω, and they all have one of ξ ≈ Ω 2 /|N | 2 or u/h|N | ≈ Ω 2 /|N | 2 . The other scales as a higher power of Ω/|N |. That is, they all require the same shear but allow for tuning such that equation ( 19) is satisfied with just one of the meridional circulation or baroclinicity. Without imposing further assumptions we therefore have

|∇Ω| ≈ Ω R , (53)
u θ < ∼ u ≈ h|N | Ω 2 |N | 2 , ( 54
) ξ < ∼ Ω 2 |N | 2 . ( 55
)
With equation (C28) we then obtain

u r ≈ h r u θ , (56)
and with equation (C10) we find

λ < ∼ Ω 2 h g . ( 57
)
As discussed in Section 7 we suggest that both the meridional circulation and baroclinicity scale similarly so that actually

|∇Ω| ≈ Ω R , (58)
u r ≈ h r u θ , (59)
u θ ≈ u ≈ h|N | Ω 2 |N | 2 , ( 60
) ξ ≈ Ω 2 |N | 2 (61) and λ ≈ Ω 2 h g . ( 62
)

HYDRODYNAMIC RAPID ROTATION

In this section we examine the case of rapid rotation with respect to the Brünt-Väisälä frequency, such that Ω |N |. We neglect magnetic fields, which we shall consider in the next section. However it is important to be careful because the relationship between the Brünt-Väisälä frequency and various convective quantities is altered in the limit of rapid rotation so, to be clear, we take |N | to be the actually realized Brünt-Väisälä frequency and |N | 0 to be what the Brünt-Väisälä frequency would be were the rotation slow and all else held constant. In this notation, the rapid rotation limit is that in which Ω |N | 0 . We are not interested in arbitrarily large rotation. In particular the system must remain primarily pressure supported and so we also require that Ω g/r. The combination of these two limits is only sensible because, unlike in radiative zones, in convection zones |N | may be significantly smaller than gr -1 .

Convection Speed

The Coriolis effect stabilizes motion perpendicular to the rotation axis. This means that in a rapidly rotating system fewer modes are unstable to convection, and those which remain unstable likely saturate at a smaller amplitude.

Many attempts have been made to estimate the strength of this effect through both closure models and numerical simulations, producing a variety of results including c ∝ Ω -1/2 (Boubnov & Golitsyn 1990;Showman et al. 2011) and c ∝ Ω -1 (Stevenson 1979;Barker et al. 2014). In section 11 we shall adopt the latter of these scaling relations, motivated by the recent convincing suite of numerical simulations by Currie et al. (2020) 5 . However, at this stage we do not need to pick a scaling and deferring this decision keeps our analysis more general so we simply write

c ≈ h|N |k, (63)
where k(Ω/|N |) is a continuous function which is order unity for Ω < |N | and which decreases asymptotically at least as fast as (Ω/|N |) -1/2 and no faster than (Ω/|N |) -1 .

Stress (Ω |N |)

In Appendix D3 we estimate the magnitude of the shear using our estimates of the convective velocity and Brünt-Väisälä frequency. In the limit of rapid rotation spherical symmetry is strongly broken, so it is no longer a guide as to how different components of the stress scale. We assume that, because there is no symmetry protection, every component of the stress feels the baroclinicity and differential rotation and shear at linear order. We thus obtain (equations D34, D35 and D36)

∇ × 1 ρ ∇ • T r ≈ ξ |N | + u h + |R∇Ω| (|N | + |R∇Ω|) + |N | 2 k 2 × h r + λ λ , ( 64
)
where k is defined by equation ( 63). Likewise, equations (D26), (D29) and (D32) give

∇ × 1 ρ ∇ • T θ ≈ ξ |N | + u h + |R∇Ω| (|N | + |R∇Ω|) + |N | 2 k 2 , (65)
and finally equations (D27), (D30) and (D33) produce

∇ × 1 ρ ∇ • T φ ≈ ξ |N | + u h + |R∇Ω| (|N | + |R∇Ω|) + |N | 2 k 2 .
(66)

Advective Terms (Ω |N |)

We now evaluate the terms which depend on and its derivatives but which are not a part of the stress. None of our analysis to determine the terms involving u in the meridional vorticity equation depended on the rotation being slow, so those may be found in appendix E. Only two terms couple

Term Magnitude

Meridional ( 16)

ω m • ∇ m u u h (Ω + |R∇Ω |) ω m u • ∇ ln ρ u h (Ω + |R∇Ω |) h r + λ λ u • ∇ω m u h (Ω + |R∇Ω |) ∇ × 1 ρ ∇ • T m k 2 | N | | N | + ξ | N | + u h + |R∇Ω | Azimuthal (17) R -1 ω φ u R u 2 h R ω φ u • ∇ ln ρ u 2 h 2 h r + λ λ u • ∇ω φ u 2 h 2 ω • ∇(ΩR) Ω |R ∇Ω | Ωω R Ω |R ∇Ω | ρ -2 (∇P × ∇ρ) φ ξ | N | 2 ∇ × 1 ρ ∇ • T φ k 2 | N | | N | + ξ | N | + u h + |R∇Ω |
Table 2. The magnitudes of the terms in equations ( 16) and ( 17) are summarized here. Factors of order unity have been dropped for simplicity.

the rotation to the e z component of the differential rotation, so

|ω • ∇(ΩR)| ≈ |RΩ∂ z Ω| (67)
and

|Ωω R | ≈ RΩ|∂ z Ω|. ( 68
)
This is the only place where we encounter an intrinsic directional preference in the coupling of differential rotation to the vorticity equation.

Baroclinicity (Ω |N |)

As well as its contribution to the stress tensor, the baroclinicity ξ enters into the vorticity balance both in relation to the scale of the meridional circulation and through the thermal wind term. So we must estimate ξ. Once more we cannot rely on symmetry arguments so we make this estimate assuming that all symmetries are maximally broken.

In Appendix C6 we thus find that (equation C60)

(ξ + k) u h|N | + k 1 + ξ + |R∇Ω| Ω ≈ 0, (69)
which is a form of the heat equation ( 22) in this limit.

Results (Ω |N |)

Our results thus far are summarized in Table 2. To proceed we must balance equations ( 19), ( 16) and ( 17) for u, ξ and R∇Ω. We do this using the methods of section 7, beginning with a heuristic argument and ending with an exhaustive search of the possibilities. The heat equation ( 19) is satisfied by having either |R∇Ω| ≈ Ω, u/h|N | ≈ k/ξ or ξ ≈ 1. We argued in section 7 that the first of these is not allowed in the limit of rapid rotation, so one of the other two must hold.

In the meridional vorticity equation u enters one order higher in Ω/|N | than all other terms. This suggests that it must asymptotically decrease like k 2 |N |/Ω in order for that equation to be balanced. This is smaller than required to satisfy the heat equation ( 19) with the meridional circulation dominant, so we find ξ ≈ 1.

In the azimuthal vorticity balance the term Ω|R∇Ω| is the Taylor-Proudman term (Hough 1897). As we show in appendix E this piece is actually Ω|R∂ z Ω|. Other occurrences of the differential rotation in this equation are sensitive to both components of the differential rotation. In this equation then, as in the other two, the shear is comparable to the dominant terms we have identified if

|R∇Ω| ≈ |N | and |R∂ z Ω| ≈ |N | 2 /Ω.
An exhaustive search of the consistent scalings with the extremal options of k ∝ Ω -1/2 and k ∝ Ω -1 reveals several other solutions. Each of these has ξ of order unity, not scaling with Ω/|N |, and has |R∂ z Ω| ≈ |N | 2 /Ω. The differences permit various tradeoffs of tuning u and the other component of the shear. Assuming as before that no tuning occurs, so that each of these is as large as it can be

ξ ≈ 1, ( 70
)
|R∇Ω| ≈ |N |, (71)
|R∂ z Ω| ≈ |N | 2 Ω (72) and u ≈ h|N | |N | Ω k 2 . ( 73
)
Note that the scaling of |R∇Ω| is consistent with the findings of Showman et al. (2011) in what they term the asymptotic regime, and corresponds to our rapidly rotating regime when all microscopic diffusivities vanish.

MAGNETIC RAPID ROTATION

We now repeat the analysis of the previous section including the effects of magnetic fields.

Magnetic Fields (Ω |N |)

If the magnetic diffusivity is vanishingly small, then the growth of the dynamo is limited only by the fact that above equipartition the field begins to quench convection (Moreno-Insertis & Spruit 1989). This simple argument predicts

A ≈ h|N |. (74)
However, both analytical arguments and numerical simulations show considerable dispersion in the scaling of magnetic field strength with buoyancy and rotation. Stevenson (1979) predicted that with fixed heat flux the magnetic field scales as Ω Given this uncertainty, we parameterize the scaling of the magnetic field instead by

A ≈ h|N |q, (75)
where q Ω | N | is a continuous function which is approximately constant when Ω |N |. This form is consistent with observations that magnetic activity is principally a function of Rossby number (|N |/Ω) [START_REF] Lehtinen | Nature Astronomy[END_REF]. In all suggested scalings of which we are aware q > ∼ k when Ω |N |, so the magnetic field becomes super-equipartion in this limit and the magnetic field energy exceeds the convective energy by a factor of (q/k) -2 .

Modifications

The arguments of sections 9.2 apply to this scenario with just one modification, namely that the scale of the turbulence is set by the Alfvén speed, which means that the stress is proportional to q 2 |N | 2 rather than k 2 |N | 2 . With this equations ( 64), ( 65) and (66) become

∇ × 1 ρ ∇ • T r ≈ q 2 |N | 2 1 + ξ + |R∇Ω| |N | + u h|N | (76) × h r + λ λ , ∇ × 1 ρ ∇ • T θ ≈ q 2 |N | 2 1 + ξ + |R∇Ω| |N | + u h|N | (77)
and

∇ × 1 ρ ∇ • T φ ≈ q 2 |N | 2 1 + ξ + |R∇Ω| |N | + u h|N | . ( 78
) 10.3 Results (Ω |N |)
The results of sections 9.3 and section 9.4 still apply. The same heuristic argument applies as in section 9, except that now the stress is a factor of (q/k) 2 larger in the meridional vorticity equation, so the meridional circulation u is enhanced by a factor of (q/k) 2 . Note that as before the meridional circulation is limited most strongly by the meridional vorticity equation and not the heat equation, so the change in the heat equation from the hydrodynamic case does not enter into our calculations. The remainder of the argument is unchanged, so

ξ ≈ 1, ( 79
) |R∇Ω| ≈ |N |, (80)
|R∂ z Ω| ≈ |N | 2 Ω ( 81
)
and

u ≈ h|N |q 2 |N | Ω . ( 82
)
Once more we performed an exhaustive search for other consistent scaling relations, this time with all combinations of k ∝ Ω -1/2 , k ∝ Ω -1 , q ∝ Ω -1/2 and q ∝ Ω 0 . All of the other consistent relations have ξ ≈ 1. The result above is recovered by assuming, as before, that no tuning occurs such that each of these achieves their greatest allowed value. Note that with the total shear scaling as |R∇Ω| ≈ |N | and the stress scaling as Ω|R∇Ω|, the overall stress is similar to that of the saturated magnetorotational instability (MRI) (Wheeler et al. 2015).

HEAT FLUX

We must account for the fact that the Coriolis effect arrests convective motions, so that at fixed |N | increasing Ω decreases the heat flux. Equivalently, at fixed heat flux |N | must increase with increasing Ω. If the rotation does not significantly affect the heat source in a star, or the outer boundary condition in a planet, |N | must change to keep F constant.

To estimate this effect when Ω |N | we use

k ≈ |N | Ω . (83)
This scaling was obtained by Stevenson (1979) through analytic closure arguments and has been observed in the extensive suites of simulations by Barker et al. (2014) and Currie et al. (2020).

Hydrodynamic

The modulus of equation ( 21) is

F = µ -1 ρc p T |Q • ∇s|. (84)
Assuming the inner product produces a factor of order unity we obtain

F ≈ µ -1 ρc p T |Q||∇s|. (85)
We may approximate the diffusivity Q by that of a random walk with velocity equal to the convective velocity and timescale given by the Brünt-Väisälä frequency so that

F ≈ µ -1 ρc p T 2 c |N | |∇s|. ( 86
)
Inserting equation (C2) and using ξ ≈ 1 we find

F ≈ ρ 2 c p T µp 2 c h|N |. ( 87
)
With the ideal gas law this simplifies to

F ≈ ρ 2 c h|N |. ( 88
)
Inserting equation ( 63) we find

F ≈ ρh 3 |N | 3 k 2 . ( 89
)
With equation ( 83) this becomes

F ≈ ρh 3 |N | 5 Ω -2 . ( 90
)
So for fixed heat flux

|N | ≈ |N | 0 Ω |N | 0 2/5 (91) and c ≈ h|N | 0 Ω |N | 0 -1/5 , (92)
consistent with what is seen in simulations (Gastine et al. 2016).

Magnetohydrodynamic

The magnetic case is more complicated. A consistent finding in both numerical simulations (Stelzer & Jackson 2013) and analytic arguments (Christensen & Aubert 2006) is that

F ≈ ρ c 2 A , (93)
which can be interpreted as the heat flux coming about primarily from advection of magnetic energy rather than entropy. In addition to simulations, equation ( 93) also agrees with observations of Jupiter (Christensen & Aubert 2006).

There are more significant differences when compared with the observed surface magnetic field of Saturn but there is reason to believe that this is not reflective of its interior (Stevenson 1982).

Next note that from equation ( 75) the magnetic acceleration is of order h|N | 2 q 2 . By contrast the Coriolis effect acting on the convection speed given by equation ( 63) produces an acceleration of order h|N |Ωk. So long as

q 2 < ∼ k Ω |N | (94)
the magnetic force is no larger than the Coriolis force, so the convection speed is limited by rotation and we expect equation ( 63) to remain valid. Equation ( 94) holds for most of the scalings we are aware of in the literature, so we shall assume that k is unchanged. We now let q = 1, and shall show that this generates the same scalings as those found by Stevenson (1979) and Starchenko & Jones (2002). Inserting q = 1 into equation (75) we expect the ratio of magnetic to kinetic energy to scale like

2 A 2 c ≈ Ω 2 |N | 2 . ( 95
)
This is consistent with simulations by Augustson et al. (2016b) and Augustson et al. (2019), which show that the ratio of magnetic to kinetic energy scales as the inverse Rossby number

2 A 2 c ≈ Ro -1 ≈ Ωh c ≈ Ω 2 |N | 2 , ( 96
)
where we have inserted equation ( 83) in the last step.

We proceed with equation ( 83). Combining this with equations ( 93), ( 63) and ( 93) we find

F ≈ ρ c 2 A ≈ ρh 3 |N | 3 k ≈ ρh 3 |N | 4 Ω -1 , (97)
so at fixed heat flux

|N | ≈ |N | 0 Ω |N | 0 1/4 , ( 98
) c ≈ h|N | 0 Ω |N | 0 -1/2 , ( 99
)
and

A ≈ h|N | 0 Ω |N | 0 1/4 . ( 100
)
This scaling for A was predicted by Stevenson (1979) and more recently by Starchenko & Jones (2002), who further produced the scalings Brünt-Väisälä frequency and convection speed in equations ( 98) and ( 99). There is significant scatter in the corresponding scalings produced by simulations, but generally they suggest a weaker scaling for c than what we find here. For instance the convection speed at constant heat flux is seen to scale as Ω -0.23 to Ω -0.29 (Christensen & Aubert 2006)10 , Ω -0.41 (Yadav et al. 2013) 11 and Ω -0.32 to Ω -0.47 (Aubert et al. 2017)12 . Similarly the scalings of the magnetic field which we discuss in section 10.1 are generally a bit weaker than that in equation ( 100), so it is possible that q ought to mildly decrease with Ω/|N |.

To bracket the possibilities note that a scaling of the form q ∝ Ω β for β < ∼ -1/2 produces an equal amount of tension in the opposite direction, and is ruled out by the scaling of the ratio of magnetic to kinetic energy Augustson et al. (2019). So it may be that the truth lies between these. Given the precision of current observations of differential rotation the difference between these choices of q is small, producing relative shear scaling like either Ω -3/4 or Ω -3/5 . If future numerical simulations pin down the scaling of magnetic field strength more narrowly we can always revisit this scaling, but for now we take q = 1 as the simpler option.

Note that with this choice we may combine equations ( 93) and ( 21) to find

Q|∇s| ≈ ρ P c 2 A . (101)
Rearranging equation (C2), inserting equation (C12) and dropping factors of order unity yields

|∇s| ≈ ρ P h|N | 2 . ( 102
)
Inserting this into the previous relation we obtain

Q ≈ c 2 A h|N | 2 (103)
and recalling that A ≈ h|N | we find

Q ≈ c h, (104)
which is the same diffusivity we used in the hydrodynamic case.

LIMITATIONS

To reiterate from section 2, we have made the following 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.

We now consider each of these assumptions and explain how they limit our results.

We have already discussed geometric factors extensively at various points. The main limitation they introduce is that we cannot easily incorporate further information about boundary conditions or the scale of the convection zone. For instance, we cannot readily predict what ought to happen in a convecting shell surrounding a differentially-rotating sphere. Near the sphere boundary effects dominate and our theory is inapplicable.

The next assumption is that perturbing forces such as tides or external sources of heating may be neglected. In single systems this is valid but in binary or planetary systems it may not be. In order for tides to be relevant the angular momentum transport they induce must be at least of order the steady state flux which is transported by the various terms which balance in the vorticity equation. Likewise, in order for heating to be relevant, it must be at least of the order that arises from the rotational perturbations to the equation of thermal equilibrium. In both cases our assumption is unlikely to be violated for stars but could fail in, for instance, hot Jupiter systems where the heat flux owing either to tides or to insolation may exceed that emerging from the centre (Jermyn et al. 2017). Similarly in a highly eccentric low-mass binary system the instantaneous tidal torque could be significant relative to convective angular momentum flux. Such scenarios are rare but likely exist.

Our third assumption enters the analysis only insofar as it allows us to use an ideal gas-type equation of state. In particular, in several places, we have used the fact that the pressure depends on both temperature and density. Removing the dependence on temperature changes the structure of these arguments significantly and so we have simply ignored such cases. This means that our analysis cannot be applied robustly to compact objects, rocky or otherwise solid bodies, or to degenerate planetary cores. In such systems though the microscopic viscosity may be quite large and convection may not be fully developed and so we would need to exclude them anyway. We further cannot apply our results to radiation-dominated regions of massive stars.

The fourth assumption is principally one of convenience: by neglecting microscopic diffusivities, we achieve significant simplifications of the equations of thermal and vorticity equilibrium. Indeed, in systems for which the third assumption holds the true momentum diffusivity is expected to be extremely small [START_REF] Spitzer | Physics of Fully Ionized Gases Sreenivasan[END_REF]). On the other hand the thermal diffusivity may not be small and near a radiative-convective boundary this assumption definitely fails. Nevertheless, for systems in which a convecting region is large enough to matter for the rotation of the system, we expect it to also be large enough that such boundaries do not dominate its dynamics. In effect this is an extension of the assumption that geometry, and hence boundary effects, are not too important.

Related to this, there is one place in which the microscopic thermal diffusivity cannot be neglected, namely near the limit of breakup rotation. As we have mentioned, the scaling laws we have derived do not hold all the way to the breakup velocity. We have shown that this is because a system which does not reconfigure to follow Keplerian rotation cannot continue to convect as the rotation approaches breakup. However the way in which convection is disrupted is by reducing the effective gravity such that a radiative temperature gradient may be convectively stable. This requires a finite radiative gradient and thence that the diffusivity does not vanish. It may be arbitrarily small because reducing the diffusivity just shifts the rotation rate at which convection ceases closer to the breakup rate but for any non-zero diffusivity there is a rotation rate at which convection fails. This may seem like a purely technical point but it is important to note because it precludes smoothly connecting the rapidly-rotating convecting solution to that of a Keplerian disc. Notably this points to one of the key open problems in understanding heat transport in W UMa 13 systems (Li et al. 2004), namely that there must be regions in which convection fails because the effective gravity vanishes.

Along similar lines, in Jupiter the microscopic conductivity changes dramatically at the depths at which ionization occurs. This could result in a transition from hydrodynamic to magnetohydrodynamic scaling. Near the transition region the behavior could be more complex than in either limit.

The fifth assumption, that of axisymmetry, is a strong one. It is responsible for a myriad of simplifications in our equations and, in particular controls, the orders of various perturbations which are protected by this symmetry. As a result any phenomena which break this symmetry may introduce new modes of heat and momentum transport which violate our calculations by an amount which is proportional to the symmetry breaking. This is a concern for systems which exhibit tides or non-axisymmetric external sources of heating. In many cases these effects are either very low in amplitude, as in a long-period binary, or very high in frequency, as in a short-period binary. In the former case they may be neglected owing to their amplitude, while in the latter they may produce no leading order effect because they are not well-matched frequency-wise to the turbulence. This was noted by Goldreich & Keeley (1977) in the context of tides, where at high frequencies relative to |N | convection only couples weakly to the tidal potential. Nevertheless, 13 i.e. low-mass contact binary there are cases in which axisymmetry strongly fails, such as in W UMa systems (Li et al. 2004) and so this assumption is worth considering carefully when applying our results.

We have already discussed the assumption that convection is subsonic and so merely note that this would only result in incorrect scaling relations if the Mach number were to exceed unity by a factor which depended strongly on Ω/|N |. That is, if the Mach number exceeds unity by a factor of a few which is set by thermodynamic considerations our analysis is unchanged but if the Mach number can increase without bound as Ω → 0 or ∞ we have a problem because then we cannot bound the convection speed by thermodynamic considerations. In fact we have shown that this is not the case because the convection speed is largely independent of Ω as Ω → 0 and decreases for fixed |N | as Ω → ∞. This assumption is therefore not one which we expect to be violated in any significant way.

Finally we must consider chemistry. We have assumed everywhere that the system is chemically homogeneous. This is actually quite likely because convection rapidly mixes chemical composition and so we do not expect to find substantial violation of this assumption unless, for instance, material is being injected into a convection zone at a rate comparable to the convective mass flux. This is a rather exotic scenario though and, with a few notable exceptions14 , does not reflect a system which is undergoing evolution on secular or nuclear time-scales so we suffer no great loss by excluding it.

SUMMARY

In this work we studied differential rotation in both slowly and rapidly rotating convection zones in both the hydrodynamic and magnetohydrodynamic limits. We obtained scaling laws for the differential rotation, the baroclinicity ξ and the meridional circulation velocity. These are summarized in Table 3. In section 11 we then incorporated scaling relations for the convection speed and magnetic field to obtain the scaling of the Brünt-Väisälä frequency with rotation rate at fixed heat flux. This allows us to put our scaling relations in terms of the non-rotating Brünt-Väisälä frequency |N | 0 which is obtained from standard stellar evolution calculations. These modified scaling relations are summarized in Table 4.

Our findings for slowly rotating systems are consistent with the solar rotation profile (Rajaguru & Antia 2015) as well as the rotation profiles of other slowly-rotating systems (Brun & Palacios 2009;Käpylä et al. 2011). Physically this results from a balance between turbulent viscosity and the Λ-effect in both equations, in agreement with arguments by [START_REF] Ruediger | Differential rotation and stellar convection. Sun and the solar stars Sackur O[END_REF]). In the azimuthal vorticity equation a comparable amount is also contributed by the thermal wind term. This is in good agreement with the work of Balbus et al. (2012), who find that thermal wind balance produces a good match to the solar rotation profile in the bulk of the solar convection zone.

The scaling we obtain for the differential rotation differs Table 3. 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. 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.

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 |) k 2 h r Ω | N | -1 ≈ h r Ω | N | -3 Ω | N | -1 k 2 ≈ Ω | N | -3 1 k ≈ | N | Ω N/A Fast MHD (Ω | N |) q 2 h r Ω | N | -1 ≈ h r Ω | N | -1 Ω | N | -1 q 2 ≈ Ω | N | -1 1 k ≈ | N | Ω q 2 k 2 ≈ Ω 2 | N | 2 MNRAS 000, 1-26 (2020)
Table 4. 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 and expanding q ≈ 1 and k ≈ | N |/Ω. 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.

Case Rüdiger et al. (1998) by a factor of Ω/|N |. This appears to be because they neglect both the thermal wind term and the effect of anisotropy on T r θ . That neglect causes their meridional vorticity balance to favour weaker relative shears for slowly-rotating systems.

| 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 from that of
The baroclinicity λ is in agreement with (albeit uncertain) measurements of the solar pole-equator temperature difference (Teplitskaya et al. 2015). Unfortunately while there have been measurements of the meridional circulation in the Sun (Zhao et al. 2013;Rajaguru & Antia 2015;Schad & Roth 2020), there are still significant disagreements between the different inversion techniques which make a direct comparison to our theory challenging.

In the rapidly rotating limit our scaling for the differential rotation is in good agreement with three-dimensional MHD simulations of rapidly rotating solar-type stars (Brun et al. 2017) and stars with convecting cores (Augustson et al. 2016a). In particular, we obtain differential rotation which increases sub-linearly with Ω (Mabuchi et al. 2015). We also find that the convection becomes increasingly magnetically dominated towards lower Rossby number. This is typically seen in these simulations (Brun et al. 2017) and arises because the Coriolis effect arrests convective motions but does not impede the growth of the magnetic field.

We obtain a similar result in the hydrodynamic regime, where we again find qualitative agreement between our predicted shear and what is found in three-dimensional hydrodynamic simulations of rapidly rotating solar-type stars (Brown et al. 2008;Matt et al. 2011;Käpylä et al. 2011), though our results disagree with at least some twodimensional simulations (Sun & Schubert 1995).

In Fig. 3 we summarize our results for the differential rotation visually alongside data for the Sun and Jupiter. These generally agree with the trends we find. Interestingly, Jupiter seems to follow our hydrodynamic scaling law in the upper regions (left, red). Further down (right, purple) Juno provides only upper bounds. The transition between these regimes occurs at a depth quite similar to that at which ionization occurs. The upper bounds are consistent with both scaling laws but [START_REF] Guillot | [END_REF] suggest that the tran-sition between the region with measurements and that with upper bounds is relatively sharp. This could mean that at the depth where the atmosphere ionizes it also transitions from hydrodynamic to MHD scaling, resulting in a steep drop in shear.

In the slow-rotation limit we did not find any order of magnitude preference for different shear directions. In that limit the angular momentum balance is chiefly between the turbulent viscosity, the thermal wind term and the Λ-effect, with all other terms scaling more slowly with angular velocity.

By contrast in the rapid-rotation limit the shear is preferentially in the cylindrical radial direction, and that preference is enforced by a factor of Ω/|N |. Specifically, in the azimuthal vorticity equation inertial or advective terms come to balance the thermal wind term, but the former preferentially couple to shear in the e z direction and hence the system tends to preferentially shear orthogonal to this. By contrast in the meridional vorticity equation both directions of the shear appear and balance against turbulent stresses, which have no such direction dependence. As a result the shear perpendicular to the rotation axis is less constrained than that along it, tending to a Taylor-Proudman state (Hough 1897) of rotation on cylinders.

CONCLUSIONS

We have made predictions for the scaling of the differential rotation, meridional circulation, magnetic field and baroclinic angle. We are aware that three of our predicted scalings have been suggested previously. These are that for the magnetic field energy (Christensen & Aubert 2006), that of the magnitude of differential rotation in the rapid hydrodynamic limit (Showman et al. 2011) and that of convection speed with rotation rate the hydrodynamic limit (Stevenson 1979). The remaining scalings are new.

Our predictions suggest that a great many details of a convecting system are irrelevant to the question of the magnitude of its differential rotation, baroclinicity and meridional circulation. All that matters to leading order is the bulk ratio Ω/|N | 0 , which is readily computed from observable parameters such as the rotation rate, surface temperature and mean density. This allows a wide range of systems to be approximately characterised with minimal data and should help to support future studies of momentum and magnetic flux transport in convection zones.

Our finding that the shear is at most of order unity in slowly-rotating systems and is suppressed in rapid rotators places a strong bound on the amount of shear even a deep convection zone may support. In particular, even for very deep regions, it is difficult to arrange for the angular velocity to vary across the zone by more than of order |N |. This is because shear of that order requires that one boundary of the zone or the other is rotating rapidly, and so the achievable shear drops off quickly in that limit.

This bound is particularly important in red giants, for which the Brünt-Väisälä frequency corresponds to periods of order 100 d or more. If the observed core-envelope shear in these systems is primarily in the convective envelope and not the interior radiative layer then the relative shear needs to be of order unity (Kissin & Thompson 2015), and this requires that the absolute shear be less than of order |N |. Because the time-scale |N | -1 is so long and the observed shear is so large, this favors scenarios which place significant shear in the radiative zones of red giants (Fuller et al. 2019).

Treating µ as a constant, we may absorb all factors other than ρ and P into the additive constant. We then set this to zero by appropriate choice of units for P and ρ, leaving just

s = 3 2 ln P ρ 5/3 , (B4)
which for a monatomic ideal gas is equal to

s = 1 γ -1 ln P ρ γ , ( B5
)
which holds more generally for any ideal gas of constant γ.

APPENDIX C: BAROCLINICITY AND CIRCULATION C1 ξ and λ

We begin with equations ( 23), ( 27) and ( 20), which may be combined to yield

λ = |∇ ln P × γ-1 γ ∇s| |∇ ln P||∇ ln ρ| = (γ -1)|∇s| γ|∇ ln ρ| ξ. (C1)
We next define the Brünt-Väisälä frequency

N 2 ≡ - γ -1 γ g • ∇s. (C2)
Expanding equation (C2) we see that

N 2 = - γ -1 γ g • ∇s = - γ -1 γ ρ ∇P • ∇s. (C3)
In the case of a convectively unstable entropy gradient this may be written as

N 2 = - (γ -1)P γ ρ |∇ ln P||∇s| ξ. (C4)
Rewriting equation (C1) in terms of N using equation ( C4) yields

λ = -(γ -1)ρN 2 P|∇ ln ρ||∇ ln P| ξ ξ . (C5)
When ∇P and ∇ρ are aligned, |∇ ln ρ| ≈ γ -1 |∇ ln P| because convection enforces a near-adiabatic relation. When they are not the density gradient projected along the pressure gradient remains adiabatic because convective motions are primarily along the pressure gradient 15 , so we instead obtain

|∇ ln ρ| = 1 γ |∇ ln P| e p + λ λ e q (C6) = 1 γ λ |∇ ln P| (C7)
With equation (C5) we obtain

λ λ = - γ ρN 2 p|∇ ln P| 2 ξ ξ . (C8)
With equation (C11) we find

λ λ = - γh 2 N 2 ρ P ξ ξ . ( C9
)
15 This is because the pressure gradient is the source of the convective restoring force. See e.g. Fig. 16 of Jermyn et al. (2018).

Inserting equation (C12) then yields

λ λ = γ hN 2 g ξ ξ , (C10)
where

h ≡ |∇ ln P| -1 (C11) = P ρg (C12)
is the pressure scale height 16 .

A further simplification we shall often use is to let ξ < ξ and λ < λ, which imply respectively that ξ ≈ 1 and λ ≈ 1. This follows because were it not the case ∇P and ∇ρ would need to be nearly perpendicular, which can only be the case near breakup rotational velocities. We address that case separately in section G.

C2 Mass Conservation

We expand the equation of mass conservation (12) as

0 = ∇ • (ρu) (C13) = u • ∇ ln ρ + ∇ • u. (C14)
We may evaluate the first term as

u • ∇ ln ρ = (u p e p + u q e q ) • ∇ ln ρ (C15) = h(u p ∇ ln p + u q e φ × ∇ ln p) • ∇ ln ρ (C16) = |∇ ln ρ|(λu q + u p λ) (C17) = (hγ) -1 λ λ u q + u p (C18)
The second may be expanded as ∇ • u = ∇ • e p u p + e q u q (C19) = e p • ∇u p + e q • ∇u q + u p ∇ • e p + u q ∇ • e q .

(C20)

Putting these together we find

0 = (hγ) -1 λ λ u q + u p + e p • ∇u p + e q • ∇u q + u p ∇ • e p + u q ∇ • e q . (C21)
Because e p is parallel to the gravitational field its variation is set by the overall scale and symmetry of the system, rather than local thermodynamic properties. The unit vector e q is then determined as a linear function of e p and so varies over the same scales. Thus, we expect that ∇ • e p ≈ ∇ • e q ≈ 1/r, where r is the spherical radial scale and enters owing to the large-scale structure of the system. Note that in the simple case where e p = e r we have

∇ • e p = 2 r (C22) ∇ • e q = cot θ r . ( C23
)
On the other hand u p and u q generally vary more rapidly than this because equation (C30) shows that the velocity is related to the local condition of thermal equilibrium. As a result we expect that the variation of u p and u q is due to the variation of local thermodynamic properties as well as global effects having to do with the scale and symmetry of the system. Along the pressure gradient then

|e p • ∇ ln u p | ≈ |e p • ∇ ln u q | ≈ |∇ ln P| + r -1 = h -1 + r -1 , (C24)
where once more r arises due to the large-scale structure of the system. In general the pressure scale height is somewhat less than the radius, so

|e p • ∇ ln u p | ≈ h -1 . (C25)
Perpendicular to the pressure gradient only the density varies, so

|e q • ∇ ln u p | ≈ |e q • ∇ ln u q | ≈ |e q • ∇ ln ρ| + r -1 ≈ λ h + 1 r , (C26)
where as before we used λ ≈ 1. Inserting equations (C24), ( C25) and (C26) into equation (C21) we find

1 h |u p | + λ + h r |u q | ≈ 0, (C27)
where we have neglected the signs of terms and just written down their magnitudes. Note that we have used γ ≈ λ ≈ 1, h r and λ 1 but have not assumed anything about the relative magnitudes of λ and h/r. Because we assume minimal geometric tuning it must be that motion along e p balances that along e q in this equation, so

|u p | ≈ λ + h r |u q |. (C28)

C3 Thermal Equilibrium

We expand u as u = e p u p + e q u q .

(C29)

Inserting this and equation (25) into equation ( 22) we obtain

u p ξ + u q ξ = - ∇ • (PQ • ∇s) P|∇s| . ( C30
)
Taking ξ ≈ 1 and inserting equations (C10) and ( 29) into equation (C30) we find

|u q |ξ 1 + h ξr + hN 2 g + |∇ • (PQ • ∇s)| P|∇s| ≈ 0. (C31)
We have assumed that the convection is efficient so |N | 2 g/h and we may drop that term to obtain

|u q |ξ 1 + h ξr + |∇ • (PQ • ∇s)| P|∇s| ≈ 0. (C32)
Examining equation ( 29) and using h < r and ξ < 1 we see that up to factors of order unity u ≈ u q , so

ξ + h r u + |∇ • (PQ • ∇s)| P|∇s| ≈ 0. (C33)

C4 Thermal Wind

To obtain equation (33) from equation ( 32) we first insert equation ( 27) to find

| ρ -2 (∇P × ∇ρ) φ | = P ρ |∇ ln P||∇ ln ρ|λ. (C34) Equation (C7) then gives | ρ -2 (∇P × ∇ρ) φ | = P ρ |∇ ln P| 2 λ γ λ . (C35)
Recalling equations ( C11) and (C12) we find

| ρ -2 (∇P × ∇ρ) φ | = g h λ γ λ . (C36)
Inserting equation (C10) we obtain equation ( 33).

C5 Slowly Rotating Perturbations (Ω |N |)

Next we analyze the effect of slow rotation on equation ( 22).

To do so we note that the diffusivity tensor is of the form (Rüdiger et al. 2005;[START_REF] Lehtinen | Nature Astronomy[END_REF])

Q ≈ h 2 |N | -1 |N | 2 Ω 2 Ω 2 |N | 2 , ( C37
)
where we have neglected multiplicative factors of order unity, the first column and row reflect e p and the second of each reflects e q . There is a second-order contribution of the form Ω|R∇Ω|/|N | 2 , analogous to that in T, but we shall argue that |R∇Ω| is no greater than Ω and so absorb that contribution into the Ω 2 terms. Similarly there is a contribution from the baroclinicity proportional to ξ, but we shall show that this is at most of order Ω 2 /|N | 2 and so may likewise absorb it into the Ω 2 terms. With equation (C37) we find

∇ • (PQ • ∇s) = ∇ • h 2 |N | P|∇s| e p (|N | 2 ξ + Ω 2 ξ) + e q (Ω 2 ξ + |N | 2 ξ) . (C38)
When Ω = 0 the system is spherically symmetric and this vanishes, so

∇ • (PQ • ∇s) = ∇ • h 2 |N | P|∇s|e p |N | 2 = 0. (C39)
Subtracting this from equation (C38) we find

∇ • (PQ • ∇s) = ∇• h 2 |N | P|∇s| e p (|N | 2 ( ξ -1) + Ω 2 ξ) +e q (Ω 2 ξ + |N | 2 ξ) . (C40)
We expand ξ ≈ 1 and ξ -1 ≈ -ξ 2 /2 as above to find

∇ • (PQ • ∇s) = ∇• h 2 |N | P|∇s| e p - 1 2 ξ 2 |N | 2 + Ω 2 ξ +e q (Ω 2 + |N | 2 ξ) . (C41)
Each term on the right-hand side of equation ( C41) is of the form

h 2 |N | P|∇s|∇ • (e i f (ξ, |N |, Ω, h, |N |, P, |∇s|)) = h 2 |N | P|∇s| (∇ • e i + e i • ∇ f (ξ, |N |, Ω, h, |N |, P, |∇s|)) , (C42)
where f is a product of powers of its arguments and i is one of p or q. As argued in the text following equation (C21), the divergence of our unit vectors is of order r -1 . The other term in equation (C42) may be expanded as

e i • ∇ f = e i • ∇ ln p ∂ f ∂ ln p ln ρ + ∇ ln ρ ∂ f ∂ ln ρ ln p (C43) = e i • - e p h ∂ f ∂ ln p ln ρ + ∇ ln ρ ∂ f ∂ ln ρ ln p . ( C44
)
Next we evaluate

∇ ln ρ ≈ λe p h + e q λ h + e q r . (C45)
The first two terms just come from our definition of λ and λ in equations ( 27) and ( 28) respectively, along with the fact that |∇ ln ρ| is of order h -1 . The final term arises because the quantities being differentiated are perturbations driven by rotation and so are sensitive to the spherical geometry. In principle there is also a term of order e p /r for the same reason, but we omit that because it is smaller than the main contribution along the pressure gradient. With λ ≈ 1 we then find

∇ ln ρ ≈ e p h + e q λ h + e q r . (C46)
Inserting this into equation (C47) we find

e i • ∇ f = 1 h e i • e p + e q λ + h r . (C47)
With this equation (C41) may be written as

∇ • (PQ • ∇s) = h |N | P|∇s| - 1 2 ξ 2 |N | 2 + Ω 2 ξ + λ + h r × Ω 2 + |N | 2 ξ . (C48)
Inserting this into equation ( 31), dropping factors of order unity and neglecting signs we obtain

ξ + h r u + h |N | ξ + λ + h r Ω 2 + |N | 2 ξ ≈ 0, (C49)
Because λ ξ in convection zones where |∇s| |∇ ln P| this reduces to

ξ + h r u + h |N | ξ + h r Ω 2 + |N | 2 ξ ≈ 0, (C50)
Dividing through by ξ + h/r we find

u + h |N | Ω 2 + |N | 2 ξ ≈ 0, (C51)
which is our condition of thermal equilibrium. Note that we could have also considered perturbations to Q owing to the differential rotation. Like those owing to the rotation these also enter at second order, and break the same symmetries, so we may modify our result to

u + h |N | (|R∇Ω| + Ω) 2 + |N | 2 ξ ≈ 0. (C52)

C6 Rapidly Rotating Perturbations (Ω |N |)

Next we analyze the effect of rapid rotation on equation ( 22).

In addition to its contribution to the stress tensor, the baroclinicity ξ enters into the vorticity balance both in relation to the scale of the meridional circulation and by means of the thermal wind term. We therefore wish to understand the various terms in the heat equation which depend on ξ.

The meridional circulation is given by equation ( 31) as

ξ + h r u + |∇ • (PQ • ∇s)| P|∇s| ≈ 0. ( C53
)
Our aim is to evaluate the second term as a function of u, ξ, Ω and |R∇Ω|.

Unlike the limit of slow rotation the divergence of the convective flux no longer vanishes by symmetry considerations, so when ξ = u = |R∇Ω| = 0 we expect that

|∇ • (PQ • ∇s)| P|∇s| ≈ Q h . ( C54
)
Without any symmetry to preclude this we expect the baroclinicity, meridional circulation and shear to perturb it at first order, so when these are non-zero

|∇ • (PQ • ∇s)| P|∇s| ≈ Q h 1 + ξ + |R∇Ω| |N | + u h|N | , (C55)
where, as usual, we have normalized the perturbations by the convective time-scale |N |.

We use mixing length theory and find the diffusivity to be

Q ≈ h c . (C56)
Using equation ( 63) we find

Q ≈ h 2 |N |k (C57) So equation (C55) becomes |∇ • (PQ • ∇s)| P|∇s| ≈ k h|N | 1 + ξ + |R∇Ω| |N | + u h|N | . ( C58
)
Inserting this into equation (C53) we obtain

ξ + h r + k u h|N | + 1 + ξ + |R∇Ω| Ω k ≈ 0. (C59)
We shall later determine that ξ is of order unity in this limit, so the factor of h/r may be dropped to yield

(ξ + k) u h|N | + 1 + ξ + |R∇Ω| Ω k ≈ 0. ( C60
)
APPENDIX D: STRESS TERMS

D1 Derivatives

In spherical coordinates we write

T = T rr T r θ T r φ T θr T θθ T θφ T φr T φθ T φφ . (D1)
Using Mathematica we expand the stress terms appearing in equation ( 3) and find

∇ × 1 ρ ∇ • T r = 1 r 2 ρ 2 ρ∂ θ T φr + T φr (cot θ ρ -∂ θ ρ) + cot θ ρ∂ θ T φθ -T φθ (cot θ∂ θ ρ + ρ) -r∂ θ ρ∂ r T rφ + 2ρ∂ θ T rφ + r ρ∂ r ∂ θ T r φ + r cot θ ρ∂ r T r φ -2∂ θ ρT r φ + 2 cot θ ρT r φ -∂ θ ρ∂ θ T θφ + ρ∂ 2 θ T θφ + 2 cot θ ρ∂ θ T θφ -cot θ∂ θ ρT θφ -ρT θφ , (D2) ∇ × 1 ρ ∇ • T θ = 1 r ρ 2 ∂ r ρ T φr + cot θ(T φθ + T θφ ) +r∂ r T rφ + 2T rφ + ∂ θ T θφ -ρr ∂ r T φr + cot θ ∂ r T φθ + ∂ r T θφ +3∂ r T rφ + r∂ 2 r T r φ + ∂ r ∂ θ T θφ (D3)
and

∇ × 1 ρ ∇ • T φ = 1 r 2 ρ 2 ρ -2∂ θ T rr -∂ 2 θ T θr -cot θ∂ θ T θr + ∂ θ T θθ + ∂ θ T φφ +r -∂ r ∂ θ T rr + ∂ r T θr + 3∂ r T r θ + r∂ 2 r T r θ +∂ r ∂ θ T θθ + cot θ ∂ r T θθ -∂ r T φφ + ∂ θ ρ r∂ r T rr + 2T rr + ∂ θ T θr -T θθ -T φφ +T θr -r∂ r ρ + cot θ∂ θ ρ + csc 2 (θ)ρ -r∂ r ρ (r∂ r T r θ + 2T r θ + ∂ θ T θθ + cot θ(T θθ -T φφ ) . (D4)
Inspection of these equations reveals that only the offdiagonal components of T contribute to the r and θ components of the vorticity equation, while all but the r φ, θr, φr, and rθ components of T contribute to the φ component of the vorticity equation. The terms highlighted in green are the ones we find to be dominant in the slowly-rotating limit.

D2 Slow Scaling (Ω |N |)

When there is rotation, shear or baroclinicity, spherical symmetry is broken. The stress tensor is then perturbed away from the symmetric form

T ≈ T rr 0 0 0 T θθ 0 0 0 T φφ . ( D5
)
It is this perturbation which gives the contribution of the stress to equations ( 13) and ( 14). However these effects break the symmetry in different ways and so enter as perturbations at different orders. In this section we shall evaluate the perturbations to ∇ × (ρ -1 ∇ • T) owing to rotation in the limit where Ω |N |. We begin with perturbations owing to rotation directly, then consider those from shear in the |R∇Ω| |N | limit, those from the meridional flow in the |∇u| |N | limit, and finally those from baroclinicity in the ξ 1 limit.

D2.1 Rotation

To leading order, the perturbation owing to rotation is of the form (Canuto et al. 1994;[START_REF] Kitchatinov | IAU Symposium[END_REF]Jermyn et al. 2018)

δT ≈ T Ω 2 | N | 2 Ω 2 | N | 2 Ω | N | Ω 2 | N | 2 Ω 2 | N | 2 Ω | N | Ω | N | Ω | N | Ω 2 | N | 2 , ( D6
)
where we have neglected dimensionless factors of order unity which multiply the various factors of Ω/|N |. The first order perturbations arise in the r φ, φr, φθ and θφ components because the Coriolis effect couples motion along other directions to motion along e φ , while the perturbations to the remaining components are second order because a second application of the Coriolis effect is required to couple motion along two directions neither of which is e φ . Finally, the perturbations on the diagonal are second order, both because of the centrifugal effect and because it takes two applications of the Coriolis effect to couple motion in a given direction to itself.

In appendix D1 we expanded the curl of the divergence of the turbulent stress. Inspection of equations (D2), ( D3) and (D4) reveals that only the off-diagonal components of T contribute to the r and θ components of the vorticity equation, while all but the r φ, θr, φr, and rθ components of T contribute to the φ component of the vorticity equation. As before we take radial gradients to produce factors of h -1 and latitudinal gradients to produce factors of r -1 +h -1 λ/ λ ≈ r -1 + λh -1 . Because λ 1 and h r, terms with fewer latitudinal derivatives are dominant.

In equation (D2) there are several terms involving one radial derivative which are perturbed at first order, so these are the most important. In equation (D3) there is a term involving two derivatives which is perturbed at first order, so that is the dominant term. In equation (D4) there is one term which is perturbed at second (leading) order and which has only radial derivatives, so that is the important term. These terms, which actually dominate in all symmetry-constrained circumstances we consider, are highlighted in green in equations (D2), (D3) and (D4). With that, we obtain

∇ × 1 ρ ∇ • T r ≈ 1 ρh 2 λ + h r T Ω |N | , (D7) ∇ × 1 ρ ∇ • T θ ≈ 1 ρh 2 T Ω |N | (D8)
and

∇ × 1 ρ ∇ • T φ ≈ 1 ρh 2 T Ω 2 |N | 2 . (D9)

D2.2 Differential Rotation

Differential rotation likewise breaks spherical symmetry. The r φ, φr, θφ and φθ components of T break this symmetry at first order because they directly couple to the shear, while the remaining terms break the symmetry at second order both by coupling to the shear twice and by coupling once to the shear and once to the rotation itself. To see this note that the differential rotation acts in the plane of e φ and the shear direction, so if motions along e r , e θ and e φ are not correlated initially then motions between e r and e θ cannot be coupled at first order by the differential rotation. Thus either the Coriolis effect is needed to couple these components or else a higher order perturbation is needed. For the same reason the differential rotation cannot perturb the diagonal components of T to first order, and two applications are needed to turn motion along, say, e r into motion along another axis and back into motion along e r . Another way to understand this is to note that the mapping φ →φ also maps Ω → -Ω and ∇Ω → -∇Ω, which may be undone by then letting φ → -φ. This spatial transformation negates the components of T which involve the direction e φ an odd number of times but not those involving it an even number of times 17 , and so the latter must be even functions of Ω. It follows that they must be at least quadratic in Ω and hence only terms of the form Ω|R∇Ω| and |R∇Ω| 2 are allowed at leading order. This is in agreement with the model of Kitchatinov et al. (1994). Based on these symmetry arguments we write the leading order term in the differential rotation

∇ × 1 ρ ∇ • T r ≈ 1 ρh 2 λ + h r T |R∇Ω| |N | , (D10) ∇ × 1 ρ ∇ • T θ ≈ 1 ρh 2 T |R∇Ω| |N | (D11)
and

∇ × 1 ρ ∇ • T φ ≈ 1 ρh 2 T Ω|R∇Ω| |N | 2 . (D12)

D2.3 Meridional Circulation

A similar argument produces the couplings to the meridional circulation. Mapping φ → -φ leaves both the meridional circulation and its shear unchanged, yet it negates all components of the stress along e φ . All terms in equations ( D2) and (D3) therefore vanish unless φ → -φ is a broken symmetry, which requires that one of Ω or |R∇Ω| be non-zero. It follows that the meridional circulation perturbs the meridional vorticity equation via the stress following

∇ × 1 ρ ∇ • T r ≈ 1 ρh 2 λ + h r T |∇u| |N | Ω |N | + |R∇Ω| |N | (D13)
and

∇ × 1 ρ ∇ • T θ ≈ 1 ρh 2 T |∇u| |N | Ω |N | + |R∇Ω| |N | (D14)
On the other hand, the same mapping of φ → -φ does not negate every term in equation (D9), for example ∂ 2 r T θr is left unchanged. Those terms may therefore be non-zero. Under the mapping e θ → -e θ these are negated, just as the relevant component of the meridional flow is negated, so they can couple to the meridional circulation at first order. It follows that the meridional circulation can appear in this 17 One might ask why this transformation is not also undone by letting z → -z. The reason is that while is a vector, Ω is a pseudovector generated by a cross-product with e φ , so its component along the e z axis is left invariant upon reflection about that axis.

stress contribution at first order without any intermediating effects, so

∇ × 1 ρ ∇ • T φ ≈ 1 ρh 2 T |∇u| |N | . (D15)
That is, u only couples to the stress via the perturbed terms in the meridional vorticity equation whereas it couples directly in the azimuthal component. Note that the couplings in equations ( D7) and (D8) are the leading order contributions to the Λ-effect, while those in equations (D10) to (D15) produce an effective turbulent viscosity with magnitude |N | -1 T [START_REF] Ruediger | Differential rotation and stellar convection. Sun and the solar stars Sackur O[END_REF][START_REF] Kitchatinov | IAU Symposium[END_REF]. These expansions are consistent with standard closure models such as those used by Gough (2012) and [START_REF] Lehtinen | Nature Astronomy[END_REF], as well as with simulations of slowly rotating convection (Käpylä et al. 2011).

D2.4 Baroclinicity

There is one further effect which may contribute to the stress at leading order, namely baroclinicity [START_REF] Ruediger | Differential rotation and stellar convection. Sun and the solar stars Sackur O[END_REF]Jermyn et al. 2018). This effect lies inside the meridional plane and hence is invariant with respect to the mapping φ → -φ. It follows that it only contributes at leading order to components of T which incorporate the direction e φ an even number of times. This is the same case we dealt with for the coupling to the meridional circulation. As a result

∇ × 1 ρ ∇ • T φ ≈ 1 ρh 2 Tξ. (D16)
In order to couple ξ into the remaining components of the vorticity equation we need other effects to break this symmetry. Hence

∇ × 1 ρ ∇ • T r ≈ 1 ρh 2 λ + h r Tξ Ω |N | + |R∇Ω| |N | (D17) and ∇ × 1 ρ ∇ • T θ ≈ 1 ρh 2 Tξ Ω |N | + |R∇Ω| |N | . (D18)

D2.5 Overall Perturbation

Putting it all together with equation ( 41) we find

∇ × 1 ρ ∇ • T r ≈ λ + h r |N | (Ω + |R∇Ω|) 1 + ξ + u h|N | , (D19) ∇ × 1 ρ ∇ • T θ ≈ |N | (Ω + |R∇Ω|) 1 + ξ + u h|N | (D20) and ∇ × 1 ρ ∇ • T φ ≈ |N | 2 ξ + u h + Ω 2 |N | 2 + Ω|R∇Ω| |N | 2 , ( D21
)
where we have used the fact from section 5 that |∇u| ≈ u/h and have added additional terms at higher order to permit a more compact representation. These are not necessarily present, though there is no symmetry which prohibits them.

D3 Rapid Scaling (Ω |N |)

When |R∇Ω| c /h the turbulent stress is dominated by convective motions. In the opposite limit it is dominated by the shear. Taking |R∇Ω| to be the frequency scale of the shear forcing and h to be its characteristic length, we write the shear turbulent velocity as so the shear never dominates the stress. We therefore have the usual convective result (Gough 1978)

T ≈ ρ 2 c . (D24)
In the rapidly rotating regime the stress is not symmetry protected, because the rotation is rapid and the offdiagonal components of T are of the same order as the diagonal [START_REF] Kitchatinov | IAU Symposium[END_REF]Jermyn et al. 2018). Hence, in the absence of shear and baroclinicity, we write

∇ × 1 ρ ∇ • T r ≈ h r + λ |N | 2 k 2 , (D25) ∇ × 1 ρ ∇ • T θ ≈ |N | 2 k 2 (D26) and ∇ × 1 ρ ∇ • T φ ≈ |N | 2 k 2 , ( D27
)
where we have followed the prescription in section 8 to evaluate the derivatives. Similarly, the convective turbulence couples viscously to the differential rotation and any meridional shear. Because there is no symmetry protection, we assume that this produces a coupling at first order. The relevant diffusivity is just the stress divided by its characteristic time-scale, which is |N | for convectively-dominated turbulence. So,

∇ × 1 ρ ∇ • T r ≈ h r + λ |N |k 2 |R∇Ω|, (D28) ∇ × 1 ρ ∇ • T θ ≈ |N |k 2 |R∇Ω| (D29) and ∇ × 1 ρ ∇ • T φ ≈ |N |k 2 |R∇Ω|. (D30)
There is likewise a first-order viscosity-like coupling to the meridional circulation of the form

∇ × 1 ρ ∇ • T r ≈ h r + λ |N |k 2 |∇u|, (D31) ∇ × 1 ρ ∇ • T θ ≈ |N |k 2 |∇u| (D32) and ∇ × 1 ρ ∇ • T φ ≈ |N |k 2 |∇u|, ( D33
)
where |∇u| is the magnitude of the tensor formed of derivatives of the velocity components. Finally, we must consider the contribution of baroclinicity to the stress. Its presence breaks no symmetries in this limit, and over its possible range from -1 to 1 it completely changes the character of the convection, so we approximate its effect as being linear and of order unity. Combining this with equations (D25), (D28) and (D31) and letting |∇u| ≈ u/h we find

∇ × 1 ρ ∇ • T r ≈ k 2 |N | |N | + ξ |N | + u h + |R∇Ω| h r + λ , (D34)
Likewise, equations (D26), (D29) and (D32) give

∇ × 1 ρ ∇ • T θ ≈ k 2 |N | |N | + ξ |N | + u h + |R∇Ω| , (D35)
and finally equations (D27), ( D30) and (D33) produce

∇ × 1 ρ ∇ • T φ ≈ k 2 |N | |N | + ξ |N | + u h + |R∇Ω| . ( D36
)

APPENDIX E: ADVECTIVE TERMS

In this appendix we evaluate the terms in the vorticity equations ( 16) and ( 17) which capture the kinematic effects associated with the rotation and circulation

E1 Meridional Equation (16)

We begin with the term ω m • ∇ m u. We may evaluate ω m by expanding equation (1) with equation ( 5) and projecting into the meridional plane to obtain

ω m = e r 2 cos θΩ + sin θ ∂Ω ∂θ -e θ sin θ 2Ω + r ∂Ω ∂r , (E1)
which does not depend on u because the vorticity is a curl and u is itself meridional. Next, we write the meridional circulation u in the basis formed by the pressure gradient and the perpendicular unit vector in the meridional plane we obtain u = u p e p + u q e q . (E2)

When the system is slowly rotating, e p ≈ e r and e q ≈ e θ , with corrections to both of order λ. Even in the limit of rapid rotation λ ≈ h|N | 2 /g 1, so we may neglect such corrections and write

ω m • ∇u ≈ e r [∂ r u r (2Ω cos θ + ∂ θ Ω sin θ) + sin θ (u θ -∂ θ u r ) (2Ω + r∂ r Ω)] + e θ [∂ r u θ (2Ω cos θ + ∂ θ Ω sin θ) -sin θ (u r + ∂ θ u θ ) (2Ω + r∂ r Ω)] . (E3)
Making the approximation that radial derivatives of u produce factors of h -1 while latitudinal ones produce factors of r -1 , taking h r and ignoring factors of order unity we find

|ω m • ∇u| ≈ 1 h [Ωu r + Ωu θ + u r |R∇Ω| + u θ |R∇Ω|] , (E4)
where we have replaced r∂ r Ω by |R∇Ω|. Equation ( 29) then tells us that u ≈ u θ u r so

|ω m • ∇u| ≈ u h [Ω + |R∇Ω|] . (E5)
The relative corrections to this are at least of order h/r, λ and Ω/|N |. Because the non-rotating system is spherically symmetric λ must be at least of order Ω/|N | too. Hence this expansion is accurate to leading order in both h/r and Ω/|N |.

We next turn to the term ω m u • ∇ ln ρ, the magnitude of which is

|ω m u • ∇ ln ρ| = |ω m ||u • ∇ ln ρ|. (E6)
The first term we may find using equation (E1) to be

|ω m | ≈ |R∇Ω| + Ω. (E7)
The remaining term we have already computed in equation (C18) and found to be

|u • ∇ ln ρ| ≈ 1 hγ u p + u q λ . (E8)
Inserting equation ( 29) and dropping factors of order unity we find

|u • ∇ ln ρ| ≈ u h h r + λ , (E9)
hence

|ω m u • ∇ ln ρ| ≈ u h (|R∇Ω| + Ω) h r + λ . (E10)
Finally, we examine the term u • ∇ω m . Again using e p ≈ e r and e q ≈ e θ we find

u • ∇ω m ≈ e r [u r (2 cos θ∂ r Ω + sin θ∂ r ∂ θ Ω) + u θ r 3 cos θ∂ θ Ω + sin θ ∂ 2 θ Ω + r∂ r Ω -e θ u r sin θ 3∂ r Ω + r∂ 2 r Ω + u θ r (r cos θ∂ r Ω + sin θ (∂ θ Ω + r∂ r ∂ θ Ω)) . (E11)
Making the same approximations as in equation (E5) we find

|u • ∇ω m | ≈ u h [Ω + |R∇Ω|] . (E12) E2 Azimuthal Equation (17)
The first term is R -1 ω φ u R . Expanding the azimuthal vorticity we find

ω φ = ∂ r u θ + r -1 u θ -r -1 ∂ θ u r . (E13)
Noting that u θ ≈ u q ≈ u we may approximate the first term by u/h. This is larger than the remaining terms, so

ω φ ≈ u h . (E14) Hence R -1 ω φ u R ≈ u R u hR . ( E15
)
Averaged over latitudes u R projects comparably on to both u q and u p , so u R ≈ u and

|R -1 ω φ u R | ≈ u 2 hR . (E16)
The next term is ω φ u • ∇ ln ρ. The first factor we have already evaluated in equation (E14) while the latter we have computed in equation (C18), so

|ω φ u • ∇ ln ρ| ≈ u 2 h 2 |N | 2 h r + λ . (E17)
We have only incurred errors of order h/r and λ/ λ ≈ λ in computing this term because we have approximated u θ by u so the expansion is accurate to leading order in both factors. The next term is u • ∇ω φ . Expanding the vorticity as in equation (E13) we find

u • ∇ω φ ≈ e φ r 2 u θ ∂ θ u θ -∂ 2 θ u r + r∂ r ∂ θ u θ +u r ∂ θ u r -u θ + r ∂ r u θ -∂ r ∂ θ u r + r∂ 2 r u θ . (E18)
This contains a term involving two radial derivatives of u θ ≈ u q ≈ u so that term dominates the expression abd we gave

|u • ∇ω φ | ≈ u 2 h 2 . (E19)
The next term is ω • ∇(ΩR). Using equation (E1) we find

ω • ∇(ΩR) = ω m • ∇(Ωr sin θ) = Ω sin θ (∂ θ Ω sin θ -r cos θ∂ r Ω) . (E20)
Neglecting the geometric factors this is just

|ω • ∇(ΩR)| ≈ Ω|R∇Ω|. (E21)
The final advective term is Ωω R e φ . Writing the vorticity in cylindrical coordinates we see that

ω = e R (-R∂ z Ω) + e z (2Ω + R∂ R Ω) + e φ (∂ z u R -∂ r u z ). (E22) Hence, |Ωω R e φ | = RΩ|∂ z Ω| ≈ Ω|R∇Ω|.
(E23)

APPENDIX F: INVERSE CASCADE

In our analysis we have assumed that solutions to the governing equations contain no geometric factors which differ significantly from being of order unity. In effect, we have assumed that any structures which form are in some fashion generic and do not depend specifically on Ω/|N |, though they may depend on the regime in which the system lies. While this assumption is usually sound, there is a known exception in the case of rapidly-rotating two-dimensional turbulence (Rhines 1973;Sukoriansky et al. 2006). This phenomenon is known as the inverse cascade (also the Rhines or enstrophy cascade) and results from the Coriolis effect preferentially scattering waves into large-scale modes. In both cases the result is feedback between small-scale convective motions and the overall geometry of the solution.

The reason that two-dimensional turbulence is relevant for our purposes is that stars and planets exhibit significant density stratification. This makes turbulence effectively twodimensional by restricting motion along the density gradient.

The Rhines cascade has been found both analytically and numerically (Danilov & Gurarie 2002) to result in the formation of alternating bands of differential rotation, also known as jets. In particular, the number of jets is seen to scale as (see equation 21 of Verhoeven & Stellmach 2014) n ≈ ΩR w , (F1)

where we have used d = R as the relevant vertical lengthscale and w is the characteristic velocity scale of turbulence in the system, given by c when the primary means of energy transport is convective. A similar effect has been observed in MHD systems and leads to a similar scaling but with w = A [START_REF] Diamond | β-Plane MHD turbulence and dissipation in the solar tachocline[END_REF], though it leads to much weaker jets (Naulin et al. 2005). The appearance of jets and the scaling in equation (F1) has been seen in a broad array of simulations (Gastine et al. 2013;Verhoeven & Stellmach 2014) and agrees well with observations of the four gas giant planets in the solar system (Ingersoll & Pollard 1982;Galperin et al. 2001).

When n > 1 the inverse cascade enhances latitudinal derivatives of quantities perturbed by the jets by a factor of n. So for instance, if the stress T is perturbed a fractional amount by the presence of jets then

1 r ∂T ∂θ ≈ T r + n T r . (F2)
Similarly, it enhances radial derivatives such that

∂T ∂r ≈ T h + n T r . ( F3
)
Taking to be no greater than 1, we see that radial derivatives are only affected when n > r/h ≈ R/h. In the slowlyrotating limit our analysis was dominated by radial derivatives, generally by a factor of R/h, so at most this effect may serve to make the latitudinal derivatives comparably important. We therefore expect no change to the scaling of any terms when Ω < ∼ |N |, because then w ≈ h|N | and n < R/h. In the rapidly-rotating limit the situation is more complicated. The dominant terms in the meridional vorticity and heat equations again contain radial derivatives, so any enhancement just rescales those equations and leaves their solution unchanged.

By contrast, in the azimuthal vorticity equation there are several terms which are not enhanced by the formation of jets. These are the Taylor-Proudman term (RΩ∂ z Ω) and the thermal wind term (∇p×∇ρ). The former is not enhanced because it involves only vertical (e z ) derivatives, while the latter is not enhanced because in this limit it is already maximized with ξ |N | 2 ≈ |N | 2 . The advective terms cannot be enhanced enough to dominate the azimuthal vorticity balance because they scale as u 2 < ∼ h 2 |N | 4 /Ω 2 , which decreases too rapidly with Ω to matter even after multiplying by an enhancing factor of (nh/R) 2 . Thus the only pieces which may be enhanced enough to matter are the stress terms.

We are not certain as to how the fluid stress responds to the introduction of a new, smaller length-scale. It is possible that the scale of the stress is unchanged but its derivatives are enhanced, leading to a greater contribution and more differential rotation. It is also possible that the stress is reduced in magnitude because its characteristic scale is shortened. We favour the latter view because, far from the critical Reynolds or Rayleigh numbers, the only sensitive length-scales are those imposed by the geometry and background gradients of the system.

Taking this to be the case, we find that the stress terms with components in the cylindrical radial direction are diminished by the same factor of nh/R by which derivatives are enhanced. Therefore the only terms whose contributions to the vorticity equation are enhanced are those which contain more derivatives along e R than stress indices along e R .

There are just two such terms which enter into the azimuthal vorticity equation as

∂ 2 R T Rz (F4)
and

R -1 ∂ R ∂ φ T φφ . (F5)
In each case there is one more radial derivative than radial index (subscript), so these contribute on the whole a factor of nh/R more.

In the hydrodynamic limit w ≈ c ≈ h|N |k and the stress terms contribute of order |N | 2 k 2 to the azimuthal vorticity equation. So after enhancement the stress contributes a net of

nh R |N | 2 k 2 ≈ h R |N | 3/2 Ω 1/2 k 3/2 . (F6)
The dominant terms are of order |N | 2 , so the ratio of these enhanced terms to the dominant ones is

h R |N | 3/2 Ω 1/2 k 3/2 |N | 2 = h R |N | -1/2 Ω 1/2 k 3/2 (F7)
Because k falls at least as fast as Ω -1/2 and h < R this ratio is less than unity, so there is no effect on our results in section 9.

In the MHD limit the jets are driven by the Reynolds (non-magnetic) stress, not by the total turbulent stress (Diamond et al. 2007, see fig. 1.3 of). This means that ≈ 2 c 2 a .

(F8)

Using equation ( 63) we see that

≈ k 2 q 2 (F9)
The MHD stress contributes of order q 2 |N | 2 to the azimuthal vorticity equation, so the terms which are enhanced are of order q 2 |N | 2 . The dominant terms are of order |N | 2 , so the ratio of the enhanced terms to the dominant ones is q 2 ≈ k 2 times the enhancement factor. With w ) so the ratio of the enhanced terms to the dominant ones is of order

≈ A ≈ qh|N |, this is nh R ≈ hΩ qR|N | , ( F10
k 2 hΩ qR|N | < ∼ h 1/2 |N | 3/2 R 1/2 Ω 3/2 1. (F11)
With Ω > |N | and h < R we see that this is less than unity, so there is no effect on our results in section 10.

APPENDIX G: BREAKUP ROTATION

Stable systems cannot rotate faster than the Keplerian velocity

Ω K ≡ g R (G1)
without invoking pressure profiles that increase outwards.

As a system approaches this velocity it nears the state of an accretion disk, in which

R∂ z Ω = 0 (G2)
and

R∂ R Ω = - 3 2 Ω. (G3)
The derivation of this state is straightforward, so we do not dwell on it. However, we are interested in how this limit is approached.

To see why this is not simply an extension of the rapidly rotating limit, note that |R∇Ω| increases sub-linearly in that limit (equations 71 and 80), such that

|R∇ ln Ω| → 0 (G4)
as Ω becomes large. This is incompatible with the Keplerian limit of -3/2 and so something different must happen in between the two limits.

The first important point is that the window of rotation rates we considered in section 10 is not infinite. In particular, we considered

|N | Ω g R = Ω K , (G5)
so the window has width

Ω max Ω min = g R|N | 2 . (G6)
For Ω > Ω max certain terms which could be ignored because h r become significant however this does not explain the difference between our results and the Keplerian limit because such terms do not scale sufficiently quickly with Ω as to produce |R∇Ω| ∝ Ω.

The second point to note is that, as the rotation rate increases, so does |N |. To incorporate both the hydrodynamic and the magnetized limits discussed in the previous section we write

|N | = |N | 0 Ω |N | 0 α (G7)
for some α > 0. Using equation (C2) we may relate this to the entropy gradient and find that

|∇s| = |N | 2 0 g ξ Ω |N | 0 2α , ( G8
)
where g is the effective gravitational field accounting for the centrifugal acceleration and the factor of ξ captures the misalignment between this and the entropy gradient. In the simple case of an ideal gas with aligned pressure and density gradients this may be related to the temperature gradient using equation ( 20

= γ h [∇ -∇ a ] , (G12)
where ∇ a and ∇ are the adiabatic and actual logarithmic temperature gradients respectively. But ∇ is bounded above by the radiative temperature gradient ∇ R , so

|∇s| ≤ γ h [∇ R -∇ a ] . (G13)
Combining equations ( G8) and (G13) we find

|N | 2 0 Ω |N | 0 2α ≤ γg ξ h [∇ R -∇ a ] . (G14)
Using equation (C12) and recalling that for an ideal gas c 2 s = γp/ρ we see that

|N | 2 0 Ω |N | 0 2α ≤ ξ γg c s 2 [∇ R -∇ a ] . (G15)
To proceed note that g in this equation is the effective 18 gravity, such that

g = g 0 -e R Ω 2 R (G16)
where g 0 is the actual acceleration of the gravitational interaction. Hence, on the equator, where g 0 is parallel to e R , g = g 0 1 -

Ω 2 R g 0 . (G17) So ξ γg 0 c s |N | 0 2 1 - Ω 2 R g 0 2 Ω |N | 0 -2α [∇ R -∇ a ] ≥ 1. (G18)
As Ω increases this breaks down because the left-hand side vanishes, both because Ω -2α and as a result of the centrifugal term. In other words at some point the system must become radiative. Of course this is only a problem in regions near the equator but it may have significant consequences for the transport of angular momentum and hence may suffice to explain the discrepancy between the Keplerian and sub-Keplerian regimes.

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 Ω); (lower-left) an example of a meridional circulation u.

Figure 2 .

 2 Figure 2. Key concepts in our theory are shown schematically: (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 ; (upperright) surfaces of constant pressure and density meet at an angle of approximately λ, the unit vector e p points along the pressure gradient and the unit vector e q is perpendicular to it.

Figure 3 .

 3 Figure 3. Our prediction for the differential rotation given by Table 4 is shown as a function of Ω/ | N | 0 and normalised by Ω. This is constant for Ω < | N | 0 and scales like (Ω/ | N | 0 ) -3/4 (MHD) and (Ω/ | N | 0 ) -3/5 (hydrodynamic) for Ω > | N | 0 . Also shown are data for the Sun(Antia et al. 2008) and Jupiter(Kaspi et al. 2018;[START_REF] Guillot | [END_REF]). The volume-weighted root-mean square shear in the Sun is shown at the average Ω/ | N | 0 . Juno measurements of the surface, at a depth of less than 3000 km, shear are shown separately from Juno upper limits on the shear deeper down. Details of the data analysis may be found in a companion paper where we focus on observational tests[START_REF] Jermyn | Convective Differential Rotation in Stars and Planets II: Observational and Numerical Tests[END_REF]).

  analogously to the scaling of convection forced at |N | in a rapidly rotating system (equation 63). The extra factor of k again comes from the Coriolis effect stabilizing motion perpendicular to the rotation axis. In section 7 we argued that |R∇Ω| < ∼ |N |. Because of this,

  ), so|∇s| = |∇ ln p -γ∇ ln ρ| (G9) = |(1 -γ)∇ ln p + γ∇ ln T |

  Aubert 2006) 6 , Ω -0.11(Yadav et al. 2013) 7 , Ω 0 (Yadav et al. 2013) 8 , and between Ω -0.02 and Ω 0.19 (Aubert et al. 2017) 9 .

1/4

from analytic growth-rate arguments. Several recent arguments suggest similar scaling laws

(Starchenko & Jones 2002;Aubert et al. 2017)

, though others obtain B ∝ Ω 0 (Davidson 2013). In numerical simulations, the field strength has been found to scale as Ω -0.02

(Christensen &

MNRAS 000, 1-26 (2020)

This estimate is for neutrally or unstably stratified media and does not generalize to the case of stable stratification.

The only exception to this is if the dynamo spontaneously MNRAS 000, 1-26 (2020)

The software used to perform this search is available at https: //doi.org/10.5281/zenodo.3967763.

The only exception to this occurs if there is turbulence which exhibits spontaneous symmetry breaking, as mentioned previously.

Jermyn et al. (2018) obtained c ∝ Ω -1/2 . The difference between our calculations there and those ofStevenson (1979) is that we did not impose the lower bound on the vertical wavenumber which they do. Such a lower bound is physically motivated for stars by the finite scale height, so we favour their scaling here.

See their equation (33).

From their dipolar fit with no magnetic Prandtl number dependence.

 8 From their dipolar fit with magnetic Prandtl number depen-

dence. 9 See their figure11(a). Their is proportional to Ω -3 and their λ ∝ BΩ -1 . MNRAS 000, 1-26 (2020)

See their equations (30) and (31).

From their multipolar fit with magnetic Prandtl number dependence.

 12 See their figure11(a). Their is proportional to Ω -3 .

For example consider hot bottom burning in Asymptotic Giant Branch stars, planet injestion, He-flashes, etc.MNRAS 000, 1-26 (2020)

Equation (C12) comes from hydrostatic equilibrium. In systems with circulation currents, stresses, and differential rotation there are deviations from this, but these are only large when the rotation is nearly Keplerian, and so we may ignore such corrections except in that regime.MNRAS 000, 1-26 (2020)

This paper has been typeset from a T E X/L A T E X file prepared by the author.18 i.e. centrifugally-corrected MNRAS 000, 1-26 (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 the IOA, 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, Chris Thompson, Douglas Gough and Steven Balbus for comments on this manuscript. ASJ is grateful to Frank Timmes and Matteo Cantiello for suggestions regarding the presentation of this work.

DATA AVAILABILITY

The software used in this analysis is available at https: //doi.org/10.5281/zenodo.3967763.

APPENDIX A: VORTICITY

To pass from equations (3) to equation ( 6) note that with equation ( 5) we have

and

Inserting these into equation (3) recovers equation ( 6).

APPENDIX B: ENTROPY

The dimensionful entropy of an ideal gas per particle is given by the Sacker-Tetrode equation as

where u is the internal energy of the gas per particle, ρ is its mass density, h is Planck's constant, k B is Boltzmann's constant and µ is the mean molecular weight of the particles (Sackur 1913;Tetrode 1912). We non-dimensionalize this by letting

The internal energy of an ideal gas is proportional to its temperature, and hence to P/ρ, so )