Dwindling Surface Cooling of a Rotating Jovian Planet Leads to a Convection Zone That Grows to a Finite Depth

Recent measurements of Jupiter's gravitational field (by Juno) and seismology of Saturn's rings (by Cassini) strongly suggest that both planets have a stably stratified core that still possesses a primordial gradient in the concentration of heavy elements. The existence of such a “diffusely” stratified core has been a surprise as it was long expected that the Jovian planets should be fully convective and hence fully mixed. A vigorous zone of convection, driven by surface cooling, forms at the surface and deepens through entrainment of fluid from underneath. In fact, it was believed that this convection zone should grow so rapidly that the entire planet would be consumed in less than a million years. Here we suggest that two processes, acting in concert, present a solution to this puzzle. All of the giant planets are rapidly rotating and have a cooling rate that declines with time. Both of these effects reduce the rate of fluid entrainment into the convection zone. Through the use of an analytic prescription of entrainment in giant planets, we demonstrate that these two effects, rotation and dwindling surface cooling, result in a convection zone that initially grows but eventually stalls. The depth to which the convective interface asymptotes depends on the rotation rate and on the stratification of the stable interior. Conversely, in a nonrotating planet, or in a planet that maintains a higher level of cooling than current models suggest, the convection zone deepens forever, eventually spanning the entire planet.


Introduction
Recent observations of the gravitational field of Jupiter by the Juno spacecraft and the seismology in Saturn's rings by Cassini suggest that neither planet is fully mixed.Instead of consisting of a compact segregated core of heavy elements surrounded by a deep well-mixed zone of convection (as is assumed in conventional models of gas giant interiors, e.g.Pollack et al. 1996), the convection zone of each planet is shallower, occupying a half or less of the planet's radius.Below the convection zone, there likely exists a stably-stratified, diffuse core with a smoother radial gradient of heavy elements (Bolton et al. 2017a,b;Wahl et al. 2017;Debras & Chabrier 2019;Mankovich & Fuller 2021;Militzer et al. 2022;Howard et al. 2023).
Despite the fact that composition gradients can be a natural outcome of formation models (e.g., Stevenson et al. 2022), the survival of such gradients over evolutionary timescales is not well understood.The expectation had been that at each planet's birth the high surface temperatures and the associated rapid cooling led to vigorous convection that quickly burrowed its way through the entirety of the planet, hence hindman@colorado.edufully mixing the planet's interior on timescales as short as 1 Myr (e.g., Müller et al. 2020).
Traditionally, thermo-compositional layers (also known as "staircases", resembling the ones observed in the artic sea on Earth) have been proposed as a mechanism to stop the growth of the outer convection zone and prevent mixing in the deeper layers of gas giants (Chabrier & Baraffe 2007;Vazan et al. 2015;Moll et al. 2017;Vazan et al. 2018).Although a convective staircase is a plausible phenomena to occur in the interior of a gas giant, hydrodynamical simulations have shown that staircases do not persist over evolutionary timescales, as multiple layers tend to merge over short timescales until a single well-mixed convective layer is left (Mirouh et al. 2012;Wood et al. 2013;Garaud 2018;Fuentes et al. 2022;Garaud 2021).
Recently, Fuentes et al. (2023) studied how a convection zone cooled from above at a fixed rate, mixes a primordial compositional gradient.In particular, they focused on how rotation modifies convective mixing at the boundary between the convection zone and the stable region below.Utilizing both 3D numerical simulations and recent scaling theory (Barker et al. 2014;Aurnou et al. 2020) that provides estimates for the speed of turbulent convective motions in a rapidly rotating fluid, they showed that rotation significantly retards the advance of the base of the convection zone, thus reducing the mixing and entrainment of heavy elements.If confirmed with more realistic simulations, rotation would provide a simple alternative mechanism to prevent mixing in gas giants.
Another possibility, that we will explore here, is that the luminosity of the Jovian planets has dwindled over time, starting at formation with a luminosity that is five orders of magnitude larger than at present (e.g., Marleau & Cumming 2014).Since convection in gas giants is driven by the fast cooling from the outer surface, this diminution of the luminous flux results in an ebbing of the convective entrainment.In this paper, we investigate the effect of convection driven by a cooling flux that decreases over time.Adapting the model of Fuentes et al. (2023), we demonstrate that if the cooling flux decreases with sufficient rapidity, the growth of the convection zone can stall and the depth of the convecting layer asymptotes to a fixed value as time advances.Further, since rotation also diminishes the rate of entrainment, for a rotating planet the rate of decay of the cooling flux can be more sedate and still lead to a stall.In section 2 we present an analytic model for entrainment of heavy fluid from the underlying region of stable stratification and in section 3 we explore how a dwindling cooling rate modifies the depth of the planet's convection zone as a function of time.In section 4, we summarize our findings and discuss the implications of our results for the Jovian planets.

The Entrainment Model
We build an analytic entrainment model which describes the depth of a gas giant's convection zone as a function of time.We do so for both a rotating planet and a nonrotating planet by following the prescription of Fuentes et al. (2023).However, here we allow the surface cooling rate to vary with time under the implicit assumption that the time scale for change in the cooling rate is on a long evolutionary time scale that is much longer than the convective overturning time.Hence, the convection is always in a state of statistical quasi-equilibrium, where the heat flux and other properties of the convection are allowed to equilibrate as the cooling flux evolves.This condition is easily met; for example, in Jupiter, the convective turnover time   = ℎ/ has a typical value of a year (Fuentes et al. 2023) and the cooling rate changes on a time scale of a million years or longer (Marleau & Cumming 2014).

Initial Atmosphere
For simplicity, consider a plane-parallel atmosphere for which the mass density  is a linearly increasing function of depth.A portion of the density gradient is due to vertical variation in the atmosphere's composition of heavy elements and the remainder arises from thermal stratification.We write this linear relation in the form, where  is the height within the atmosphere (with  = 0 corresponding to the top of the atmosphere) and  0 and Γ are positive constants that represent characteristic values of the density and the reciprocal of the density scale height.

Entrainment
Surface cooling will cause a convection zone to form at the upper surface and this zone will deepen with time as convection scours the interface between the convection zone and the stably stratified fluid below.Heavy fluid will be dredged upward and mixed into the convection zone.If we assume that the convection zone is well-mixed to an adiabatic density gradient,  0 Γ ad , and at a given time has a depth of ℎ(), the change in the gravitational potential energy from the initially unmixed state, Δ, is given by where  is the buoyancy frequency, which for small density fluctuations about the fiducial density  0 is a constant value given by  2 =  (Γ − Γ ad ) where  is the gravitational acceleration (assumed constant).
The entrainment hypothesis states that the rate of change of potential energy is proportional to the kinetic energy flux within the convective motions (Linden 1975).Dimensional analysis dictates that the kinetic energy flux is proportional to the convective flow speed, , times the kinetic energy density,  0  2 /2.Hence, we find a relationship between the speed of the convective motions and the rate at which the convective interface descends, where  is a constant of proportionality called the mixing efficiency.Equation (3) has been well validated by both laboratory experiments and numerical simulations (e.g., Turner 1968;Fernando 1987;Molemaker & Dĳkstra 1997;Fuentes & Cumming 2020;Fuentes et al. 2023).The mixing efficiency  depends on the strength of turbulence of the flow, ranging from 0.1 in experiments with salty water, to approximately 1 in astrophysical flows (Fuentes & Cumming 2020).Since fluid motions in gas giants are highly turbulent (low viscosity), we adopt  = 1.In a nonrotating system, the convective flow speed, , can be estimated by using mixing length arguments.We start by making three assumptions: 1) the convective kinetic energy arises from buoyant acceleration over the entire depth of the convecting layer, 2) the density fluctuations  within the convection are proportional to the thermal perturbations , and 3) the convective heat flux equals the rate of surface cooling .These three assumptions lead to the following three relations where the constant   and  are the the specific heat capacity at constant pressure and the coefficient of thermal expansion, respectively.When these three equations are combined, one finds that the convective velocity in the nonrotating system,  NR , scales with the cube root of the cooling rate, When convection occurs in a rotating system, the importance of rotation is quantified by the Rossby number Ro, defined as the ratio of the rotational period to the convective turnover time.In a rapidly rotating system, one with Ro ≪ 1, mixing-length theory leads to a very different scaling.Instead of a pure balance between inertia and buoyancy, one expects CIA balance (e.g., Stevenson 1979;Barker et al. 2014;Aurnou et al. 2020), which is a three way balance between the Coriolis, inertial, and buoyancy (Archimedean) forces.This leads to a convective velocity,  R , that is significantly reduced compared to a nonrotating system, In both Jupiter and Saturn Ro ∼ 10 −6 , so we expect CIA balance to hold and convection to be strongly constrained by rotation.
By inserting these expressions for the convective velocity into Equation (3), we obtain ODEs that relate the depth of the convective layer to the cooling rate.For a nonrotating planet, we derive whereas, for a rotating planet, we obtain

Evolution of the Cooling Rate
When first formed, a gas giant is extremely hot and cools rapidly through radiation.As the planet ages, its surface temperature falls and the cooling rate slows.Figure 1𝑎 illustrates the radiant luminosity of Jupiter as a function of time, as calculated by Marleau & Cumming (2014)1 using the MESA stellar evolutionary code (Paxton et al. 2011(Paxton et al. , 2013)).As noted by Marleau & Cumming (2014), the luminosity dwindles at a rate roughly consistent with the reciprocal of time, i.e.,  ∝  −1 .The red dashed curve in Figure 1𝑎 illustrates this power law dependence.We have chosen the constant of proportionality such that  = 8.7 × 10 −10  ⊙ at the planet's current age  = 4.6 × 10 9 years, with  ⊙ = 3.846 × 10 33 ergs s −1 .
Over the same evolutionary timescales, the rotation rate of a gas giant changes only moderately.When young and luminous, magnetic braking spins down the planet (e.g.Takata & Stevenson 1996;Batygin 2018) gravitational contraction due to cooling slowly spins it up.However, after about a million years the gas giant has contracted to a density where degeneracy pressure opposes further contraction (Stevenson & Salpeter 1977), thus, the planet's radius and rotation rate stop changing.Since, the planet's rotation rate varies by factors of order unity over the entire period after planetary formation (e.g., see Batygin 2018), we assume a constant rotation rate.Similarly, we ignore changes in the planetary radius.

Constant Cooling Rate
For a point of comparison, first consider a cooling rate  that is temporally steady.For such cooling, Equations ( 9) and (10) can both be integrated analytically to provide the depth of the convection zone as a function of time (Fuentes et al. 2023), , (12) where  is an overshooting length scale that characterizes the depth to which convective plumes penetrate locally into the stable region.The integration constant ℎ 0 provides the depth of the convection zone at time  = 0.
In the limit of long times, for the nonrotating fluid, one recovers the well-known result that the layer grows with a square-root of time dependence.However, in a rotating fluid, the growth rate is slower; the layer advances as a power law with an index of 5/12 (Fuentes et al. 2023), We emphasize that in both the nonrotating and rotating planet, the convection zone continual deepens and never reaches a finite asymptotic value.

Dwindling Cooling Rate
Now consider a cooling rate that dwindles like the reciprocal of time,  =  0 ( 0 /).We adopt values of  0 and  0 to match the red-dashed line in Figure 1𝑎:  0 = 6.0 × 10 8 ergs cm −2 s −1 and  0 = 4.0 × 10 4 such that Jupiter has a luminous flux  = 4 2   that is equal to its current intrinsic luminosity ( = 8.7×10 −10  ⊙ ) at its current age ( = 4.6×10 9 years) and a luminosity of 10 −4  ⊙ at  =  0 .We adopt a constant Jovian radius of   = 7.15 × 10 9 cm.For this power-law form for the energy flux, the ODEs, Equations ( 9) and ( 10), can be integrated analytically giving,  16) and ( 17) and using the parameter values that appear in Table 1.The convection zone in a nonrotating Jupiter deepens rapidly engulfing the entire interior of the planet in less than 10 5 years.Conversely, the convection zone in a rotating giant planet initially grows, but asymptotes to a finite depth.For parameter values appropriate for Jupiter, the convection zone stops growing after it occupies the outer 37% of the planetary radius.
Figure 1𝑏 illustrates the depth of the convection zone as a function of time for both a nonrotating planet (red curve) and for a rotating planet (blue curve).We use parameter values that are appropriate for the interior of Jupiter, see Table 1.For the nonrotating planet, the convection zone quickly grows, engulfing the entire Jovian interior in less than 10 5 years.The gray region of the diagram indicates states for which the planet is fully mixed.In the rotating planet on other hand, the convection zone initially grows rapidly but then tapers off approaching a constant asymptotic value (which is marked by the horizontal dashed blue line).If we assume that the planet starts with a very shallow convection zone, ℎ 0 ≪   , we find that for large times,  ≫  0 , the asymptotic depth has the following value, For the parameter values indicated in Table 1, one obtains an asymptotic depth that is only a fraction of Jupiter's radius, ℎ ∞ ≈ 0.33   .Thus, due to a combination of Jupiter's rapid rotation and its dwindling luminous flux, its convection zone will never completely engulf the planet.A central core remains stably stratified and unmixed for all time.
We note that in a realistic model of a gas giant, the buoyancy frequency  will be a function of depth determined by the net stratification of the dilute core.In the analytic model that we illustrate in Figure 1𝑏, we adopted a constant value of  ∼ 10 −4 s −1 within the stably stratified fluid below the convection zone.Within the convection zone itself,  2 ≤ 0 by definition.This value corresponds to  ∼ 0.2  dyn , where is the natural frequency of the planet.The fraction of 0.2 results by averaging the buoyancy frequency in recent models of dynamical tides in Jupiter (Idini & Stevenson 2022;Lin 2023) over the interior below the convective interface at 0.7  .Models using recent equations of state for hydrogen-helium mixtures at Jovian conditions (e.g., Militzer et al. 2022) give values that are 3-5 times larger (interestingly, similar values can be estimated for Saturn, e.g., Figure 1b in Mankovich & Fuller 2021).
Table 1.Thermodynamic and buoyancy properties of Jupiter and Saturn (order of magnitude estimates).
Given the current uncertainty in the value of the buoyancy frequency in the outer regions of a gas giant's primordial atmosphere, in Figure 2 we present isocontours of the asymptotic depth as a function of the buoyancy frequency  of the interior stable region and the rotation rate of the planet Ω.For comparison the current rotation rate of Jupiter is indicated by the horizontal dotted line.The dashed blue curves indicate isocontours for values of the asymptotic depth that are a hundredth and a tenth of Jupiter's radius: 0.01  and 0.1  .The solid gray curve indicates where the convection zone asymptotically reaches a depth equal to Jupiter's radius and the gray region marks the parameter regime where the interior of the planet eventually becomes completely mixed.The region that is shaded red corresponds to asymptotic depths that fall within the range ℎ ∞ ∈ [0.3   , 0.5   ].
It is important to note that recent studies of Jupiter's gravitational moments by the Juno mission (e.g., Militzer et al. 2022) have constrained the upper convection zone to span the outer 30% of the planet's radius (being consistent with our calculations).However, Militzer et al. (2022) also found that inner regions of the planet ( < 0.4  ) can be convective and well-mixed (presumably due to dissolution of core material into the hydrogen-rich envelope).Although we do not investigate the inclusion of an inner convective zone, the gravity data does not rule out the possibility of a core that is stable to convection.

Summary and Discussion
In a gas giant, surface cooling by radiation drives a convection zone that deepens due to the entrainment of fluid from the stably-stratified fluid that lies beneath.Given that we expect a gas giant to be formed with a radial gradient in the concentration of elements heavier than helium (with heavy elements more heavily concentrated in the center), the growing depth of the convective interface should homogenize the planet's composition by dredging up fluid rich in heavy elements and mixing it throughout the convection zone.We made a simple analytic model that describes convective boundary mixing and the downward propagation of a convective interface.This model was designed to illustrate the effects of rotation and temporally varying surface cooling, hence to avoid unnecessary complications, we assume that the fluid is nearly incompressible.Using Jupiter-like values, we have shown that in a nonrotating planet, the convection zone grows rapidly and, within less than a million years, the convection zone should have ingested the entirety of the planet's interior.In a rapidly rotating planet with surface cooling that dwindles with time, the convection zone grows and saturates to a finite value, mixing only the upper ≈ 33% of the planet's radius.However, we make clear that final size of the outer convection zone is sensitive to parameters that are not well-constrained in the planet's interior, particularly to the value of the buoyancy frequency (ℎ ∞ ∝  −5/6 ).Therefore, our models should be viewed as illustrative instead of quantitative.

Why Does the Convection Zone Stop Growing?
A simple energy argument can be used to understand how the convection zone can asymptote to a finite depth.The energy required to entrain and mix the fluid ultimately comes from the surface cooling.Cooling at the surface causes contraction and generation of gravitational potential energy in a surface boundary layer.The fluid in this boundary layer continually rains downward converting the potential energy into kinetic energy and this kinetic energy is used to perform the work needed to entrain and mix the dense fluid from the underlying layer into the convection zone.
In a nonrotating star, the flux of kinetic energy generated by falling cool fluid is directly proportional to the total amount of thermal energy that has been radiated away, i.e.,  0  3 ∝ see Equation ( 7).Hence, if the cooling rate is integrable, a finite amount of energy is available for mixing and we should expect the deepening of the convection zone to stall at a finite depth.If we consider a general power law for the temporal behavior of the cooling rate,  ∝  − , we see that the power law index  must be greater than 1 for the total radiated energy to be finite, Figure 2. Isocountours of the asymptotic depth of the convection zone within a rotating planet with a surface cooling rate that varies like ∼  −1 .The asymptotic depth is shown as a function of the buoyancy frequency  of the stably stratified interior and the planet's rotation rate Ω.To ease interpretation, we have indicating the rotation rate of Jupiter with the dotted horizontal line.The two blue dashed lines provide isocontours for ℎ ∞ = 10 −2   and for ℎ ∞ = 10 −1   .The red shaded region indicates all values consistent with current constraints for Jupiter's outer convection zone ℎ ∞ ∈ [0.3   , 0.5   ] (e.g., Militzer et al. 2022).All isocontours were calculated using Equation ( 17) and Table 1.
A power law of  = 1 matches evolutionary models of Jupiter's luminous flux (see Figure 1𝑎).Hence, the depth of the convection zone in a nonrotating Jupiter is likely to be unbounded.Even though the value of the asymptotic depth depends on the stratification , the mixing efficiency , and other properties of the fluid ( 0 , ,   , etc.), whether or not an asymptotic depth exists only depends on the temporal behavior of the cooling rate.If the cooling flux has a finite integral over all time, the convection zone will stall at a finite depth.Since rotation inhibits the generation of kinetic energy from surface cooling, the condition in a rotating planet is subtly different.The kinetic energy flux is proportional to  6/5 instead of simply the flux .Thus, the convection zone will achieve a finite asymptotic depth if the integral of  6/5 exists and is finite, The criterion for a finite amount of lifting work for the rotating case is that  > 5/6, which is satisfied by our previous choice of  = 1.Once again the stratification and the fluid properties do not matter.Surprisingly, beyond the fact that the planet is rotating rapidly (Ro ≪ 1), the rotation rate also does not matter.All that matters is whether  6/5 is temporally integrable.

Implications for Convective Staircases
Here we have concentrated on how rotation and temporally decreasing surface cooling reduce the entrainment of fluid across a convective interface, thus leading to a slowing of the downward advance of that interface.Both of these mechanisms (rotation and dwindling cooling rate) will reduce the speed and penetration depth of convective overshoot into the stable region below the outer convection zone.Undoubtedly, these mechanism should also be important for the existence and longevity of a convective staircase that might reside underneath the outer convection zone.Rotation has been shown to hinder double-diffusive instabilities, suppressing the formation of a convective staircase (Moll & Garaud 2017).However, staircases can also form by different mechanisms, such as thermohaline-shear instabilities (Radko 2016;Garaud 2017) and convective destabilization of the diffusive boundary layer below the interface (as observed in laboratory experiments and numerical simulations Turner 1968;Fernando 1987;Fuentes et al. 2022).In such cases where a staircase forms despite the rotation, rotation and dwindling cooling should make staircases more stable by reducing the chance of disruption by overshooting motions from the vigorous outer convection zone.

Saturn and the Ice Giants (Uranus and Neptune)
Although we have primarily focused on Jupiter, our findings may also be relevant for Saturn, Uranus, and Neptune.For Saturn, seismology of its rings has revealed internal gravity modes trapped within a large scale composition gradient within a region that occupies the inner half of the planet's radius or larger (Fuller 2014;Mankovich & Fuller 2021).Hence, Saturn too is not fully mixed.Using Equations ( 16) and ( 17), and planetary parameters appropriate for Saturn (see Table 1), we obtain an asymptotic convection-zone depth that is 12% of Saturn's radius, ℎ ∞ = 0.12   .The shallower depth that we derive for Saturn is largely a consequence of the fact that we employ a larger buoyancy frequency for Saturn.The density stratification (and hence the buoyancy frequency) is more constrained in Saturn than it is in Jupiter because ring seismology provides reliable and independent estimates for the internal density profile (Fuller 2014;Mankovich & Fuller 2021).Even when the final size of outer convection zone is smaller in our model (versus 0.3-0.4 , see Mankovich & Fuller 2021), we remind the reader that our analytic model is qualitative.
For Uranus and Neptune, fully-mixed adiabatic models have difficulties explaining their observed luminosities (Fortney et al. 2011;Nettelmann et al. 2013;Scheibe et al. 2019), suggesting that their interiors may possess stably-stratified layers with a non-uniform distribution of heavy elements (Vazan & Helled 2020;Helled et al. 2011;Scheibe et al. 2021).Moreover, incomplete radial mixing can also explain the complex magnetic fields of the two ice giants (Stanley & Bloxham 2004, 2006).If the convection zones of either Uranus or Neptune spanned the entire interior of the planet, we would expect primarily dipolar planetary magnetic fields.Thus, the fact that observed magnetic fields for the ice giants are dominated by the quadrupole contributions (e.g., Connerney et al. 1987;Connerney et al. 1991), suggests that the convection zone only occupies an outer shell.Since all of the giant planets are rapidly rotating, with thermal luminosities that decrease in time in a way qualitatively similar to Jupiter (Fortney et al. 2011), the same mechanisms that we explore here for Jupiter may result in a convection zone that occupies only a fraction of the interior of the other giant planets.This work was supported by NASA through grants 80NSSC18K1125, 80NSSC19K0267, and 80NSSC20K0193.

Figure 1 .
Figure 1.Panel (a): Radiant luminosity  as a function of time for a MESA model of a Jupiter-like planet (Marleau & Cumming 2014) that undergoes cooling for 5 Gyr (black solid line) and for an analytic model where  ∝  −1 (dashed red line).Panel (b): Depth of the convection zone of Jupiter as a function of time for the analytic models described by Equations (16) and (17) and using the parameter values that appear in Table1.The convection zone in a nonrotating Jupiter deepens rapidly engulfing the entire interior of the planet in less than 10 5 years.Conversely, the convection zone in a rotating giant planet initially grows, but asymptotes to a finite depth.For parameter values appropriate for Jupiter, the convection zone stops growing after it occupies the outer 37% of the planetary radius.