Rotation Reduces Convective Mixing in Jupiter and Other Gas Giants

Recent measurements of Jupiter’s gravitational moments by the Juno spacecraft and seismology of Saturn’s rings suggest that the primordial composition gradients in the deep interior of these planets have persisted since their formation. One possible explanation is the presence of a double-diffusive staircase below the planet’s outer convection zone, which inhibits mixing across the deeper layers. However, hydrodynamic simulations have shown that these staircases are not long-lasting and can be disrupted by overshooting convection. In this Letter, we suggests that planetary rotation could be another factor for the longevity of primordial composition gradients. Using rotational mixing-length theory and 3D hydrodynamic simulations, we demonstrate that rotation significantly reduces both the convective velocity and the mixing of primordial composition gradients. In particular, for Jovian conditions at t ∼ 108 yr after formation, rotation reduces the convective velocity by a factor of 6, and in turn, the kinetic energy flux available for mixing gets reduced by a factor of 63 ∼ 200. This leads to an entrainment timescale that is more than 2 orders of magnitude longer than without rotation. We encourage future hydrodynamic models of Jupiter and other gas giants to include rapid rotation because the decrease in the mixing efficiency could explain why Jupiter and Saturn are not fully mixed.


Introduction
The history of the solar system is written in the interiors of the giant planets.In particular, since Jupiter accounts for 75% of the Solar System's planetary mass, understanding its interior (or more precisely, determining the distribution and total mass of heavy elements within the planet) can provide important clues on the first stages of the formation of the solar system (see the recent review by Miguel & Vazan 2023).Traditionally, Jupiter has been considered a fully convective planet, with a well-defined interface separating a dense core made of heavy elements and the hydrogen-helium envelope (Pollack et al. 1996).However, new data from the Juno mission indicates that Jupiter does not have a well-defined central core of heavy elements, nor is it homogeneously mixed (e.g., Bolton et al. 2017b,a;Wahl et al. 2017;Debras & Chabrier 2019;Howard et al. 2023).Interestingly, ring seismology suggests that Saturn is likewise not homogeneously mixed (Mankovich & Fuller 2021).These new observational constraints on the internal structure of gas giants require updates to formation and evolution models (Helled et al. 2022).Composition gradients decrease the efficiency of the heat transport throughout the planet's interior and thus affect its cooling, and mass-radius relationship at a given age (e.g., Chabrier & Baraffe 2007;Leconte & Chabrier 2012, 2013).
Primordial composition gradients are expected to naturally result from planet formation (e.g., Helled & Stevenson 2017;Müller et al. 2020;Stevenson et al. 2022) or collisional events during the planet's evolution (Liu et al. 2019), but the persistence of composition gradients against strong convectivemixing over evolutionary timescales is poorly understood.One mechanism for maintaining a composition gradient appears in one-dimensional evolution models, where convective staircases (i.e., multiple convective layers separated by thin, stably-stratified, diffusive interfaces) form underneath the outer convective envelope, preventing mixing deeper in the interior (e.g., Vazan et al. 2015Vazan et al. , 2018;;Stevenson et al. 2022).However, 3D simulations of layered convection in incompressible flows have shown that convective staircases fully mix on short timescales (e.g., Wood et al. 2013;Garaud 2018).Recently, Fuentes et al. (2022) explored the formation of a convective staircase beneath a growing convection zone, and found that compositional mixing across the interface between the convection zone and the stable region prevents the formation of the staircase.Similar results were found by Anders et al. (2022) in the context of stellar convection, concluding that compositional gradients are ineffective barriers against convective mixing over evolutionary timescales.
Rotation is another factor that could prevent convective mixing in Jupiter and other gas giants.Experiments and nu-F .
merical simulations of thermal convection have shown that rotation reduces the convective velocities when the flow has small Rossby number, Ro, the ratio of the rotation period to the convective turnover time (e.g., Stevenson 1979;Fernando et al. 1991;Barker et al. 2014;Dietrich & Wicht 2018;Aurnou et al. 2020;Hindman et al. 2020).Further, numerical experiments wherein rapidly-rotating convection zones grow via entrainment have shown that the vertical transport of buoyancy and kinetic energy across the convective layer, and thus the entrainment rate, are reduced by rotation (Julien et al. 1996a).Convective flows in Jupiter and other gas giants can be significantly slow compared to rotation (Ro 10 −6 , e.g.Guillot et al. 2004), so the effects of rotation on the erosion of primordial composition gradients are likely to be significant.Indeed, Leconte & Chabrier (2012) suggested that rotation would hamper convection and impede the mixing of heavy elements.
In this Letter, we examine how rotation modifies the rate of mixing of a stable fluid layer by a neighboring convection zone.In Section 2, we use mixing-length theory to argue that rotation reduces the convective velocity and the kinetic energy flux available for mixing.Then, in Section 3 we use the scalings from rotating mixing length theory to show that rotation decreases the rate of entrainment and mixing at the boundary between the stable and unstable layers.In Section 4, we use 3D numerical simulations to support the predictions from the previous sections.Finally, in Section 5 we discuss the implications of our results for Jupiter and other gas giants.

Mixing-length theory for rotating convection
In this section, we summarize the mixing-length scalings for the convective velocity and kinetic energy flux in nonrotating and rotating flows.
In the absence of rotation, convective motions are nearly isotropic with a characteristic length scale, ℓ, that is roughly the depth of the convection zone (e.g., in a laboratory experiment) or a pressure scale height (e.g., in stellar and planetary interiors).Further, the velocity scale of the convective flows is approximately the overturning buoyant "freefall" velocity,  NR = √︁  NR ℓ, where the subscript NR means no rotation,  is the acceleration due to gravity, and  = −( ln /)|  is the coefficient of thermal expansion at constant composition , and  NR is the characteristic temperature perturbation of convective parcels.The heat flux carried by the convective motions is   ∼    NR  NR , where  is the mass density, and   is the specific heat capacity at constant pressure.Then, the convective velocity can be estimated as Typical values for Jupiter's deep interior at  ∼ 10 8 yrs after formation, when the convective luminosity is high enough so that a large fraction of the planet is expected to be mixed, ℓ ∼   ∼ 10 9 cm (a pressure scale height),   / ∼ 10 10 cm,  ∼ 1 g cm −3 , and   ∼ 6 × 10 5 cgs (e.g., Cumming et al. 2018;Helled et al. 2022), give  NR ∼ 40 cm s −1 (and accordingly, the convective turnover time is ℓ/ NR ∼ 0.8 yr).
In the limit of rapid rotation, convection is quasigeostrophic, and convective flows are set by a balance between the Coriolis, inertial, and buoyancy (Archimedean) forces (called CIA balance, e.g., Stevenson 1979).In this limit, a new length scale, ℓ ⊥ , associated with motions perpendicular to the rotation vector, , emerges (e.g., Barker et al. 2014;Aurnou et al. 2020;Vasil et al. 2021).From CIA balance, one obtains where the subscript R means rotating and Ω is the angular frequency of the rotation vector.Using the definition of the heat flux   ∼    R  R , and eliminating  R and ℓ ⊥ from Equations (2), it follows that Plugging in the previous values for Jupiter, together with 2/Ω ∼ 10 hours, we obtain  R ∼ 6 cm s −1 .Rotation therefore reduces the convective velocity by a factor of ∼ 6; then, the convective turnover time increases by a factor of ∼ 6.6, giving ℓ/ R ∼ 5 yrs.
Although the effect of rotation on the convective velocity could be considered minor in the sense that the resulting turnover time is still many orders of magnitude smaller than the age of the planet, the important quantity for entrainment is the kinetic energy flux, which is reduced by a much larger factor.This can be seen from the ratio between the kinetic energy flux  K ∼ (1/2)  3 conv associated with the convective motions and the thermal heat flux carried by convection.Without rotation, Equation (1) gives  K,NR /  = (1/2)  3 NR /  ∼ (ℓ/  ).In the limit of rapid rotation, Equation (3) gives where Ro =  R /2Ωℓ 1 is the Rossby number of the flow based on the rotationally-constrained convective velocity.Then, for a given heat flux, rotation reduces the kinetic energy flux by a factor of Ro 1/2 1.Using the values above for Jupiter, we estimate Ro ∼ 10 −5 , and consequently, the kinetic energy flux available to mix the fluid is reduced by a factor of ∼ 4 × 10 −3 .This reduction is important because the energy source for the upward transport and mixing of heavier material is the kinetic energy of the convective motions.We investigate this in more detail in Section 3.

Entrainment in rapidly rotating flows
Mixing of material (entrainment) across convective boundaries (e.g., a core-envelope interface) has been investigated with great detail in laboratory experiments and numerical simulations (e.g., Turner 1968;Fernando 1987;Molemaker & Dĳkstra 1997;Fuentes & Cumming 2020).In those studies, a convection zone advances into a stable region, and the returning flows from overshooting convection entrain fluid with the chemical composition of the stable layer.This material rapidly mixes, so the convection zone grows over time (see Anders & Pedersen 2023, for a review of boundary mixing processes).By reducing the convective velocities, rotation reduces the entrainment (mixing) efficiency, and consequently, decreases the growth rate of the convective layer.
Consider a plane-parallel fluid layer of depth , stablystratified with composition decreasing linearly with height,  0 () = | 0 /|(−).Note that  could represent the spatial variation in mass fraction of heavy elements, or the mean molecular weight of the fluid, as would be induced by salt in water or heavy elements in hydrogen.Then, a heat flux  0 escapes through the top of the stable layer, such that convection is driven from above and mixes the composition gradient from top to bottom.For simplicity, we adopt the Boussinesq approximation (Spiegel & Veronis 1960), wherein the fluid density is constant everywhere except for buoyant perturbations, which is a decent approximation in the vicinity of a convective boundary.
We now estimate the rate of entrainment ℎ/, where ℎ() is the depth of the convection zone at time .Lifting and mixing a quantity Δ of heavy fluid requires an energy expenditure Δ = (1/6)  0 Δℎ 2 (see, e.g., Fuentes & Cumming 2020), where mass conservation dictates that Δ = 0.5| 0 /|ℎ.Here,  0 is the density of the layer, and  = ( ln /)  is the coefficient of compositional expansion (both quantities assumed constant).Differentiating Δ with respect to time, we write the rate of change of potential energy due to mixing as  (Δ)/ = (1/4)  0 | 0 /|ℎ 2 ℎ.We adopt the entrainment hypothesis, which states that the rate of change of potential energy across the layer is set by the kinetic energy flux of the convective motions  KE = (1/2)  0  3 conv (Linden 1975).It follows that Equation ( 5) assumes that all the kinetic energy is used to do work and lift material from below, and that thermal effects do not significantly affect the rate of change in the potential energy across the layer (a good approximation since Δ Δ).
We see that given  conv , Equation ( 5) can be integrated to obtain ℎ().Without rotation, we use  =  0 ,   =  0 , and ℓ = ℎ and assume  conv ≈  NR ∝ ℓ 1/3 (Equation 1), to obtain (within a factor of order unity due to the assumptions above) A full derivation including the mixing efficiency, interfacial heat transport, and thermal effects on the density can be found in Fuentes & Cumming (2020).
where  = √︁ | 0 /| is the buoyancy frequency.From Equation (6), the rate at which the convection zone expands inward depends on how quickly it cools down (which sets the convective velocity magnitude) and on the magnitude of the initial composition gradient.
On the other hand, in the limit of rapid rotation we should use  conv ≈  R ∝ ℓ 1/5 (Equation 3), and integration of Equation (5) gives We can see that rotation slows down the propagation of the convection zone, it reduces the power-law index from 0.5 to approximately 0.41, but the most significant effect comes from the prefactor ( ) −1/12 , which can be very small over long timescales.
It is instructive to compare the mixing timescales in both cases.From Equations ( 6) and ( 7), the time required for convection to mix a stable layer of vertical size  (i.e., ℎ() = ) is in the non-rotating case, and in the rotating case.Therefore, the mixing timescale increases approximately by a factor of (5/6)Ro −1/2 , a commensurate change to the decrease in the kinetic energy flux-compare Equations ( 4) and ( 10).
In a system like Jupiter where Ro ∼ 10 −5 at  ∼ 10 8 yrs after formation, this corresponds to a mixing timescale that is more than two orders of magnitude longer.This difference is substantial, e.g., for  ∼ 10 9 cm,  2 ∼ /  ∼ 10 −6 s −2 , Equations ( 9) and ( 10) give  NR ∼ 10 8 yrs and  R ∼ 10 10 yrs, i.e., for the rotating case the mixing timescale would be more than the age of the Solar System.

Numerical experiments and results
In the following, we use two hydrodynamic simulations (with and without rotation) to test the theory above.We model 3D thermal convection in a plane-parallel binary mixture of fluid (heavy and light elements).The domain initially consists of two layers: a convection zone of uniform composition, on top of a region stabilized by a composition gradient.The top of the convection zone is constantly cooled so that convection propagates downwards.We focus on how rotation affects the mixing of heavy material at the lower convective boundary.
In our simulations, as in our derivation above, we assume that the length scale of the fluid motions is much smaller .The overline denotes that variables are averaged over the horizontal directions.Note that u 2 =  2 +  2 +  3 , where , , and  are the , , and  components of the velocity.The results are shown at at a times where the convection boundary is located at  = 0.75 (i.e., ℎ = 1.25) for both runs.This time corresponds to  = 3750 for the non rotating case, and  = 11500 for the rotating case.The background colors distinguish between the convection zone and the stable region.than a density scale height, so that the effects of curvature can be ignored and the Boussinesq approximation is valid (Spiegel & Veronis 1960).Although the density in giant planets can vary by many orders of magnitude, this approximation fully captures nonlinear mixing near the convective boundary, which is the focus of our study.In the rotating case, we consider uniform rotation throughout the fluid, and assume that gravity and rotation are aligned and point in the z direction,  = Ω ẑ, and g = − ẑ, i.e., the simulation is at a polar latitude.The density of the fluid depends on both composition and temperature fluctuations.In what follows, all results are presented in dimensionless form.For details on the numerical setup, nondimensionalization, and code, we refer the reader to Appendices A and B.
Figure 1 shows 3D snapshots of the  component of the vorticity for both the non-rotating and rotating simulations.Initially, in both cases, the convection zone spans the upper half of the domain,  ∈ [1, 2].After significant evolution ( ∼ 3750 for the non-rotating case, and  ∼ 11500 for the rotating case), the bottom of the convection zone has advanced to  = 0.75, mixing the primordial composition gradient.Note that in the rotating case, it takes approximately 3 times longer for the zone to increase its vertical extent by Δℎ = 0.25.We show later that this difference agrees with Equation (10).Rotation visibly alters the flow morphologies.In the nonrotating case, the characteristic flow length scale is similar in all directions, while in the rotating case, the horizontal length scale of the flow is much smaller than the vertical scale (as expected from the Taylor-Proudman theorem, e.g., Proudman 1916;Taylor 1917).
Figure 2 shows horizontally-averaged profiles of the composition (left), vertical convective composition flux (middle), and vertical kinetic energy flux (right); fluxes are depicted at the same instant in time as displayed in Figure 1.Rotation does not affect the composition profile, which is homogeneous in the convection zone with a transition to the primordial linear distribution in the stable region.This is expected because mass conservation requires that the amount of material transported and mixed in the convection zone depends only on the thickness of the convective layer, Δ = 0.5| 0 /|ℎ.The initial composition in the convection zone is  = 0.5, and with | 0 /| = 1 and Δℎ = 0.25 giving Δ = 0.125, the expected evolved convection zone composition is  = 0.625, which is what we observe (see left panel).In the rotating case, both the convective composition flux and the kinetic energy flux are smaller by a factor of ≈ 3 when compared with the values of the non-rotating case (see middle and right panels), supporting both our hypotheses that rotation decreases the vertical transport of dense material and that this reduction is related to the kinetic energy flux.Although we only focus on the compositional transport and kinetic energy flux, it is worth mentioning that the thermal structure of the fluid across the box is similar to that of the composition.The only difference is that in the non-rotating case, the convection zone tends to be isothermalized by strong turbulence, while in the rotating case, a small temperature gradient tends to be sustained in the convection zone (this has been reported in previous studies of thermal convection in rotating flows, see, e.g., Julien et al. 1996b;Cheng et al. 2020).
We expect from the entrainment rates calculated in Section 3 that the convection zone will grow more rapidly in the nonrotating case than the rotating case.We quantify this difference by measuring the size of the convection zone ℎ as a function of time.We define the bottom of the convection zone as the location where the convective composition flux is maximized (see Figure 2, middle), so that ℎ is the distance between that location and the top of the box.Figure 3 (a) shows ℎ vs. time for the two cases; as expected, ℎ in the non-rotating case grows more rapidly than in the rotating case.We fit a power-law function of the form ℎ() = [ℎ 1/ 0 +  ( −  0 )]  , where ℎ 0 is the size of the convection zone at initial time  0 , and  is a constant.In both simulations, we demand that ℎ 0 = 1, since initially the upper half of the box is mixed.The time  0 is formally the time when the convection zone starts to grow.This only occurs after the convective turbulence is fully established.The initial state in our simulations is one of quiescence.In the non-rotating case, the convection continues to form and equilibrate up to  0 = 225 and only then does entrainment of the stable fluid begin.The convection in the rotating case takes longer to be established,  0 = 325.We find  ≈ 1.61 × 10 −4 (6.62 × 10 −5 ) and  ≈ 0.49 (0.4) in the non-rotating (rotating) case.These values are close to what is predicted by the simplified model in Section 3, where the predicted powerlaw index is  = 1/2 = 0.5 for the nonrotating case, and  = 5/12 ≈ 0.41 for the rotating case.In dimensionless form, the predicted values of the prefactor are  = 4(  /Pe ff ) ∼ 10 −4 for the non-rotating case, and  = (24/5)  (Ro ff /Pe ff 2 ) 3/5 ∼ 2×10 −5 for the rotating case, where the input dimensionless parameters   , Pe ff , and Ro ff , are the density ratio, and the Peclet and Rossby numbers based on the free-fall velocity, respectively.The reader is refered to Appendix A for more details on the nondimensionalization.
Finally, we measure the Rossby number of the flow, Ro sim =  /2Ωℎ, where  and ℎ are the magnitude of the flow velocity averaged within the convection zone and the thickness of the convection zone, respectively.We find Ro sim ≈ 0.065, which decreases slightly as the convection zone evolves and deepens (see Figure 3b).When compared with the ratio between the kinetic energy flux in the rotating and non-rotating case (as predicted by mixing-length theory, see Equation 4), we find Ro sim 1/2 /(  ,R /  ,NR ) ≈ 0.87.Finally, our theory (Equation 10) predicts that mixing a depth Δℎ should take longer in the rotating case by a factor of (5/6)Ro −1/2 .In our simulations this difference is approximately a factor of 3, close to (5/6)Ro sim −1/2 ≈ 3.3.

Summary and Discussion
In this work, we studied the effect of rotation on the compositional transport across the boundary between a convection zone and a stable region (e.g., a core-envelope interface in a gas giant).Mixing length theory predicts that rotation reduces the convective velocity (compare Equations 1 and 3).For Jovian conditions at  ∼ 10 8 yrs since formation, this corresponds to a velocity reduction of approximately a fac-F .
tor of 6.Consequently, the kinetic energy flux available for the entrainment and mixing of heavy material gets reduced by a factor of ( NR / R ) 3 ∼ Ro −1/2 ∼ 6 3 ∼ 200 (see Equation 4).Using the prescriptions for the convective velocity from mixing-length theory, we derived a simple analytical model for the rate of entrainment and the propagation of a convection zone.The thickness of the outer convection zone follows ℎ() ∝  1/2 in the non-rotating case (Equation 6), and ℎ() ∝ Ω 1/4  5/12 in the rotating case (Equation 7), supporting the hypothesis that rotation decreases compositional mixing.From this model, the timescale to mix a stable layer of vertical size  is larger by a factor of (5/6)Ro −1/2 , where Ro is the Rossby number, when rotation is taken into account in the convective velocity (Equation 10).We ran two numerical experiments (including and not including rotation) of convective entrainment and found a good agreement between the analytic predictions and the simulations (e.g., see Figures 2-3).The agreement between our numerical simulations and the mixing-length theory predictions is encouraging.The mixing time can increase by two orders of magnitude for Jupiter-like values at  ∼ 10 8 years after formation.This could significantly reduce the rate of core erosion, preventing strong mixing without the need to invoke a double-diffusive staircase.For example, Moll et al. (2017) derived that the rate at which core mass is transported into the envelope is given by  core =  −1  core /  , where  core is the luminosity at the core-envelope interface, and  −1 is ratio of the buoyancy flux associated with heavy-element transport to the buoyancy flux associated with heat transport.Using 3D simulations of nonrotating incompressible flows, Moll et al. (2017) constrained  −1 0.5 for Jupiter conditions.In the entrainment model investigated in this work, the rate at which mass is transported into the outer convection zone is  core ∝ ℎℎ/.Then, it follows that rotation reduces the erosion rate by a factor of  core,R /  core,NR =  −1 R / −1 NR ≈ (/2Ω) 1/2 (24/5) −1/6 (where we have used Equations 6 and 7).Using values for Jupiter at  ∼ 10 8 years after formation, we estimate  −1 R / −1 NR ≈ 0.01.Although we have derived this ratio with a rather simple model, the fraction is significantly smaller than unity and illustrates the importance of rotation for convective boundary mixing.
Although our numerical experiments are done in a Cartesian box and assume that gravity and rotation are parallel to each other (as expected at the poles of the planet), the reduction in the convective velocity does not depend on the exact geometry or shape of the composition profile.However, numerical simulations with misaligned gravity and rotation have shown that thermal mixing is not uniform across latitude (e.g., Dietrich & Wicht 2018;Currie et al. 2020).A logical next step would be to investigate this problem using spherical geometry and a primordial composition gradient that is more consistent with formation models (e.g., Stevenson et al. 2022).
Our simulations do not include density stratification, as we adopt the Bousinessq approximation which formally limits the vertical scale to be much less than a density scale height.Including density stratification and compressibility will change the mixing rate by a combination of two effects.First, the potential energy required to mix a region should decrease; while composition is uniformly mixed in the convection zone in both cases, the compressible case mixes the density to the adiabatic gradient, whereas Boussinesq convection mixes the density to a constant value, which requires a greater energy expenditure.Second, stratification increases the asymmetry between upflows and downflows and in turn increases the net kinetic energy flux, so it should matter whether a stable entraining region is at the top or bottom of a convection zone (see e.g., Fig. 4 Meakin & Arnett 2007).At this time we cannot speculate on the net result of those effects.To our knowledge, there are no simulations of convective entrainment accounting for both composition gradients and rotation in a compressible flow.We will investigate this problem in a future paper.
We encourage future work on hydrodynamic models of gas giants, as improved mixing prescriptions at convective boundaries provide a promising avenue for better aligning 1D evolution models with observational inferences of gas giant interiors from NASA's missions Juno and Cassini.

Figure 1 .Figure 2 .
Figure 1.Volume renderings of the  component of the vorticity ω = ∇ × u for the non-rotating case (left) and rotating case (right).For both cases, the snapshots are shown at times when the bottom of the convection zone has expanded beyond the initial location (from  = 1 to 0.75).The Rossby number in the rotating case is Ro sim ≈ 0.065, computed using the averaged flow velocity in the convection zone, and thickness of the convection zone.

Figure 3 .
Figure 3. Panel (a): Thickness of the convection zone ℎ as a function of time .The results are shown for both the rotating and non-rotating cases.The solid lines correspond to the best fits of the form ℎ () = [ℎ 1/ 0 +  ( −  0 ) ]  for each case.We obtain  ≈ 0.49 for the non-rotating case (compared to the analytic prediction  = 0.5) and  ≈ 0.40 for the rotating case (compared to the analytic prediction  = 5/12 ≈ 0.41).See text for more details.Panel (b): Ro sim , Ro sim 1/2 , and   ,R /  ,NR , as a function of the thickness of the convection zone, ℎ.Note that Ro sim 1/2 has similar values to   ,R /  ,NR , giving support to Equation (4).