Draft version May 14, 2024

Typeset using LATEX twocolumn style in AASTeX631

Constraints on the long-term existence of dilute cores in giant planets

A. Tulekeyev,1 P.Garaud,1 B. Idini,2 and J. J. Fortney2
1 Department of Applied Mathematics, Baskin School of Engineering, University of California, Santa Cruz
1156 High St, Santa Cruz, CA 95064 USA
2 Department of Astronomy and Astrophysics, University of California, Santa Cruz

1156 High St, Santa Cruz, CA 95064 USA

arXiv:2405.06790v1 [astro-ph.EP] 10 May 2024

ABSTRACT
Ring seismology has recently revealed the presence of internal gravity waves inside Saturn that extend
up to 60% of Saturn’s radius starting from the center, in what is recognized today as Saturn’s stably-
stratified dilute core. Similarly, gravity measurements on Jupiter suggest the existence of a dilute core
of still poorly constrained radial extent. These cores are likely in a double-diffusive regime, which
prompt the question of their long-term stability. Indeed, previous DNS (Direct Numerical Simula-
tions) studies in triply-periodic domains have shown that, in some regimes, double-diffusive convection
tends to spontaneously form shallow convective layers, which coarsen until the region becomes fully
convective. In this letter, we study the conditions for layering in double-diffusive convection using
different boundary conditions, in which temperature and composition fluxes are fixed at the domain
boundaries. We run a suite of DNS varying microscopic diffusivities of the fluid and the strength of the
initial stratification. We find that convective layers still form as a result of the previously discovered γ-
instability which takes place whenever the local stratification drops below a critical threshold that only
depends on the fluid diffusivities. We also find that the layers grow once formed, eventually occupying
the entire domain. Our work thus recovers the results of previous studies, despite the new boundary
conditions, suggesting that this behavior is universal. The existence of Saturn’s stably-stratified core,
today, therefore suggests that this threshold has never been reached, which places a new constraint on
scenarios for the planet’s formation and evolution.

Keywords: hydrodynamics – instabilities – planets and satellites: gaseous planets – planets and satel-
lites: interiors – turbulence

1. INTRODUCTION lute or fuzzy core that extends up to at least ∼50% of

Simple formation scenarios via the “core accretion” the planetary radius in the case of Saturn (Fuller 2014;
mechanism predict the existence of a compact mass of Mankovich & Fuller 2021), and up to 30% to 70% of the
heavy elements concentrated near the center of Jupiter radius in the case of Jupiter (Wahl et al. 2017; Debras
and Saturn, in what is known as the traditional view & Chabrier 2019; Miguel et al. 2022; Idini & Stevenson
of gas giant planet interiors (Safronov 1969; Perri & 2022; Militzer et al. 2022; Howard et al. 2023). These
Cameron 1974; Mizuno et al. 1978; Mizuno 1980; Pol- cores are composed of a mixture of heavy elements and
lack et al. 1996; Guillot 2005). Over the past several H/He fluid, with an abundance of heavy elements that
years, new high-precision constraints on the gravity field increases towards the planetary center, thus creating a
of Jupiter from the Juno Mission (Bolton et al. 2017; strongly stabilizing composition gradient.
Folkner et al. 2017; Durante et al. 2020), and the grav- The case for Saturn is especially compelling, as Sat-
ity field and seismic oscillation frequencies of Saturn urn’s interior is constrained not just by the gravity
from the Cassini Mission (Hedman & Nicholson 2013, field but also by ring seismology (Marley & Porco
2014; Hedman et al. 2018; French et al. 2019; Iess et al. 1993; Fuller 2014; Mankovich 2020; Mankovich & Fuller
2019), have revolutionized our understanding of the in- 2021). Mankovich & Fuller (2021) have recently found
terior structure of these planets, providing observational strong evidence for a dilute core structure that is stably-
evidence that puts this traditional view into question. A stratified, meaning that it admits internal gravity wave
wealth of evidence suggests that these planets have a di- normal modes, and is not freely convecting. Similar
2

seismic constraints are not available for Jupiter, but we where ν is the kinematic viscosity of the fluid, κC and
have instead more precise measurements of the gravity κT are compositional and thermal diffusivities, respec-
field thanks to the Juno mission. All interior models tively, and NT2 and NC2 are the Brunt-Väisälä frequen-
published to date that fit this gravity field, by differ- cies due to the thermal and compositional stratification,
ent groups using several different high pressure equa- respectively. While the Prandtl number and the diffu-
tions of state for H/He mixtures, require a dilute core sivity ratio characterize properties of the fluid, the in-
(Miguel et al. 2022; Idini & Stevenson 2022; Militzer verse density ratio, R−1 , defines the strength of strati-
et al. 2022; Howard et al. 2023). However, the true ex- fication, with more strongly stratified regions having a
tent of Jupiter’s core, and whether it is statically sta- larger R−1 . The criterion for linear instability to ODDC
ble today, remains unclear due to the relatively large is (Walin 1964; Kato 1966):
uncertainties in our knowledge of the equation of state
Pr + 1 −1
for H-He and how they complicate the interpretation of 1 ≤ R−1 ≤ ≡RC . (2)
Jupiter’s gravity field (Howard et al. 2023). Pr + τ
Interior models that aim to reconstruct the planet’s in- The limit R−1 = 1 corresponds to the Ledoux criterion,
terior structure from the gravity field or ring seismology and a region of fluid with R−1 ≤ 1 is unstable to over-
are inherently static models. Furthermore, such mod- −1
turning convection. A region with R−1 ≥ RC on the
els typically have little power to constrain the thermal other hand is linearly stable to ODDC.
structure of the planetary interior. However, additional In general ODDC-induced transport is weak, charac-
information and constraints can be obtained by taking terized by a total1 temperature flux FT that is at most
into account the stability of the fluid to turbulent mo- ten times larger than the radiative flux, and a total
tions and the effect that mixing caused by these motions compositional flux FC that is at most a hundred times
would have on the heavy element distribution. For in- larger than pure diffusion (see figure 4 of Mirouh et al.
stance, while static models may allow a dilute core to 2012). However, direct numerical simulations (DNS)
support a strong compositional gradient, that solution have shown that ODDC can also transition into layered
could be unphysical from a fluid dynamical perspective convection under some circumstance (Rosenblum et al.
depending on the strength of the temperature gradient. 2011), and the fluxes in the layered state can be orders
Indeed, should the latter be strong enough to drive some of magnitude larger than in the ODDC state (Wood
form of convection despite the stabilizing compositional et al. 2013). The layering transition is caused by the γ-
gradient, the turbulent motions could rapidly homog- instability, first discovered by Radko (2003) in context
enize the heavy elements throughout the core. In that of fingering convection in the Earth’s ocean, and applied
respect, it is therefore crucial to establish the stability of to ODDC by Rosenblum et al. (2011) and Mirouh et al.
dilute cores to convective mixing, either via overturning (2012). The γ-instability can cause the growth of large-
convection or double-diffusive convection. scale density anomalies in double-diffusive convection
In what follows, we therefore study the properties of a (ODDC or fingering convection) because temperature
Saturn-like dilute core, with a destabilizing temperature and composition are not transported at the same rate
gradient, but a strongly stabilizing heavy-element gradi- in these small-scale instabilities, but both contribute to
ent (so that the core overall is stable to the Ledoux crite- the total density.
rion). This type of stratification can nevertheless be un- To model this instability, we usually define the so-
stable to oscillatory double-diffusive convection (ODDC called ‘inverse flux ratio’ γ −1 as 2 :
hereafter) in some cases, as well as secondary convec-
tive instabilities via layering (see more on this below). FC
γ −1 = . (3)
ODDC, which is often also called semi-convection, is a FT
form of weak convection triggered by diffusively desta-
In homogeneous ODDC, γ −1 is only a function of P r,
bilized gravity waves (Walin 1964; Kato 1966; Baines &
τ , and R−1 . Rosenblum et al. (2011) and Mirouh et al.
Gill 1969). It is primarily governed by three nondimen-
(2012) showed analytically and numerically that the γ-
sional parameters: the Prandtl number P r, the diffusiv-
instability only takes place when γ −1 is a decreasing
ity ratio τ , and the inverse density ratio R−1 (Rosen-
function of R−1 . For a given set of fluid parameters, the
blum et al. 2011), defined respectively as:

1 By total here we mean the sum of the turbulent flux and the
diffusive flux.
2 This corresponds to γ −1 in Mirouh et al. (2012) but we omit the
ν κC −1 NC2 tot
subscript ‘tot’ to simplify the notation.
Pr = ,τ = ,R = , (1)
κT κT NT2
 3

function γ −1 (R−1 ) is always concave (see for instance time, i.e. so-called ‘fixed-flux’ boundary conditions, and
figure 5 of Mirouh et al. 2012), decreasing from R−1 = 1 are the same at both boundaries. While this is not neces-
−1
to a global minimum at RL , and then increasing again sarily more realistic than periodic boundary conditions,
at larger R−1 . The position of this minimum depends the comparison with previous works will establish which
on P r and τ . This implies that convective layers are results are robust and which are model-dependent.
expected to form by the γ-instability whenever In section 2, we present the model setup and the nu-
merical method of solution. In section 3, we present a
−1
1 ≤ R−1 < RL (P r, τ ). (4) suite of selected simulations that span a range of initial
stratifications and analyze the results in the context of
As shown by Mirouh et al. (2012) and Wood et al.
the γ-instability theory. Finally, in the section 4, we dis-
(2013), these convective layers are initially shallow but
cuss the implications of our results for planet formation
seem to always rapidly merge to occupy the entire ther-
and evolution models.
mally unstable region. The implication of these results,
as discussed in Wood et al. (2013), is that ODDC-
2. MODEL SETUP
unstable regions are not necessarily long-lived: if R−1
−1
drops below RL , the γ-instability rapidly triggers the For the purpose of our study, we consider a Cartesian
formation of shallow layers that then gradually merge domain of height H under uniform gravity g = −gez .
into a large-scale, chemically homogeneous convection We assume that the domain is bounded by imperme-
zone. able parallel plates at z = −H/2 and z = H/2, and is
However, it is important to note that all of these pre- horizontally periodic in both x and y directions. The
vious studies were very idealized. In particular, they equations are the usual Boussinesq equation for com-
assumed triply periodic boundary conditions, in which pressible gases (Spiegel & Veronis 1960) which are valid
the thermal and compositional fluxes in the domain can as long as fluid velocities are much smaller than the lo-
increase freely in response to layer mergers (Wood et al. cal sound speed and H is much smaller than the density,
2013). By contrast, in a real planet these fluxes would be pressure, or temperature scale heights3 . These are:
set by processes that take place on the much larger plan-
etary lengthscale and on the evolutionary timescale. As ∇ · u = 0, (5)
a result, an important question is whether the shallow ∂u 1
convective layer formation and gradual mergers happen + u · ∇u = − ∇p + (αT T − αC C)gez + ν∇2 u,
(6)
∂t ρm
as described above if the fluxes are externally fixed in- ∂T dTad
stead. + u · ∇T − w = κT ∇2 T, (7)
∂t dz
In giant planets, for instance, the temperature flux ∂C
at a given radius is due to the loss of thermal energy, + u · ∇C = κC ∇2 C, (8)
∂t
initially generated from the planet’s formation and sub-
sequent additional gravitational contraction. The lo- where u = (u, v, w) is the velocity field, p is the pres-
cal temperature flux, up to the surface that radiates sure perturbation away from hydrostatic equilibrium, T
to space, therefore depends on the integrated local heat is the temperature perturbation away from the mean
capacity and contraction rate, such that values nears temperature Tm of the domain, and C is the composi-
the planet’s center will necessarily differ from that at tion perturbation away from the mean composition Cm .
the surface.(e.g., Nettelmann et al. 2015, their Figure Note that the composition C is simply the mass frac-
17). The compositional flux in the center arises by con- tion of heavy elements (usually called Z), but we use C
vective erosion of the compact core while the one at here to avoid any confusion with the coordinate z. We
the surface is zero. These realistic boundary conditions assume that the fluid has a constant mean density ρm ,
could, in principle, be modelled in a full-planet simu- kinematic viscosity ν, thermal diffusivity κT , and com-
lation. However, this is not computationally possible positional diffusivity κC . It also has constant coefficients
owing to the high resolution required to capture the dy- of thermal expansion and compositional contraction, αT
namics of ODDC, so we are forced to model a small and αC respectively, defined as:
region within the dilute core instead, of extent at most
a few hundred meters in each direction (Moll & Garaud 1 ∂ρ 1 ∂ρ
αT = − , αC = . (9)
2017). The question of which boundary conditions to ρm ∂T ρm ,Cm ρm ∂C ρm ,Tm
apply at the edge of the computational domain is much
less clear. So, in this study, we assume for simplicity 3 Note that this implies that our domain only spans a small portion
that fluxes of temperature and composition are fixed in of the core.
4

Finally, dTad /dz is the adiabatic temperature gradient at z = ±Ĥ/2, where hats are used to denote nondimen-
and is also assumed to be constant. sional dependent variables.
Both plates are assumed to be no-slip. We impose The nondimensional equations are:
fixed flux boundary conditions, namely, a positive tem-
perature flux FT , and a positive composition flux FC at ∇ · û = 0, (20)
each plate so that: ∂ û
+ û · ∇û = −∇p̂ + P r(T̂ − Ĉ)ez + P r∇2 û,(21)
∂t
∂T FT ∂C FC
=− and =− (10) ∂ T̂
∂z κT ∂z κC + û · ∇T̂ − ŵ = ∇2 T̂ , (22)
∂t
at z = ±H/2. In the absence of fluid flow, the steady- ∂ Ĉ
state solutions of (7) and (8) with boundary condi- + û · ∇Ĉ − R0−1 ŵ = τ ∇2 Ĉ, (23)
∂t
tion (10) have constant temperature and composition
where the independent variables (x, y, z, t) are now im-
gradients dTb /dz = −FT /κT and dCb /dz = −FC /κC
plicitly nondimensional as well, but the hats are omitted
throughout the domain. We then define the fluctua-
to avoid crowding the equations. The usual parameters
tions away from that steady background, T̃ and C̃, such
for double-diffusive convection emerge, namely
that :
dTb dCb ν κC αC |dCb /dz|
T =z + T̃ , C = z + C̃. (11) Pr = , τ= , R0−1 = .
dz dz κT κT αT |dTb /dz − dTad /dz|
The new boundary conditions are then: (24)
The Prandtl number, P r, is typically of order 10−2 or
∂ T̃ ∂ C̃
= 0, and =0 (12) 10−1 , while the diffusivity ratio τ is often an order of
∂z ∂z
magnitude smaller than P r (French et al. 2019). Both
at each plate, and the perturbation equations are are fluid properties. By contrast, R0−1 depends on the
thermal and compositional stratification of the fluid, and
∇ · u = 0,
for the simulations presented here, is the inverse density
∂u 1
+ u · ∇u = − ∇p̃ + (αT T̃ − αC C̃)gez + ν∇2 u, ratio at the start of the simulation.
∂t ρm
The initial conditions are:
∂ T̃ dTb dTad
+ u · ∇T̃ + w −w = κT ∇2 T̃ , û = 0, T̂ = 0, Ĉ = ϵ, (25)
∂t dz dz
∂ C̃ dCb where ϵ is a small amplitude random noise. To study
+ u · ∇C̃ + w = κC ∇2 C̃. (13)
∂t dz this model numerically, we use ’Dedalus’, an open-source
Next, we use the standard nondimensionalization for Python library for spectral methods (Burns et al. 2020).
double-diffusive convection (Radko 2013), where the The computational domain is Cartesian, and configured
unit length [l], the unit time [t], the unit velocity [u], to use Fourier bases in both horizontal x and y direc-
the unit temperature [T ] and the unit composition [C] tions, and a Chebyshev basis in the vertical direction
are given by z. The nondimensional domain size is 100 × 100 in the
 1/4 horizontal plane, and is of height Ĥ = 200 vertically.
κT ν
[l] = d = , (14) 3. SIMULATIONS RESULTS
αT g|dTb /dz − dTad /dz|
d2 We ran a number of simulations with two different val-
[t] = , (15)
κT ues of the Prandtl number P r and diffusivity ratio τ . In
κT this section, we present a number of sample simulations
[u] = , (16)
d at P r = τ = 0.1, that demonstrate typical behaviors and
[T ] = d|dTb /dz − dTad /dz|, (17) dynamics observed in all performed simulations listed in
αT Table 1. In each case, we show the evolution of the total
[C] = d|dTb /dz − dTad /dz|. (18)
αC nondimensional mean density profile, ρ̄(z, t), defined as
In these new units, the top boundary is located at the :
nondimensional height z = H/2d = Ĥ/2 , and the bot- Z t+∆t Z Lx Z Ly
1
tom boundary at z = −H/2d = −Ĥ/2. The nondimen- ρ̄(z, t) = ρ̂(x, y, z, t)dxdydt
2∆tLx Ly t−∆t 0 0
sional boundary conditions become: (26)
where ρ̂ = −(T̂ − z) + (Ĉ − R0−1 z) is the total nondimen-
∂ T̂ ∂ Ĉ
û = v̂ = ŵ = 0, = 0, and = 0, (19) sional density, and ∆t ≃ 50. The time integral is useful
∂z ∂z
 5

to filter out fast gravity waves dynamics. We also com- eventually overturns into convection. Later, this con-
pute the local inverse density ratio, R−1 (z, t), defined vective layer grows to fill the entire domain around
as: t ≈ 75, 000, and the system is now in a statistically-
∂ C̄(z,t)
stationary fully-convective state. The mechanism by
R−1 (z, t) = ∂z
, (27)
∂ T̄ (z,t) which the convective layer grows will be explored in a
∂z
future publication.
where C̄ and T̄ are defined similar to ρ̄.
We find that our simulations fall into 4 categories de- 3.3. Intermediate R0−1
scribed below.
As we continue to weaken the initial stratification
3.1. Highest R0−1 (i.e. lower R0−1 ), we find that the long-term evolution
We begin with the simplest case which takes place of the system is very different, and notably, that con-
at R0−1 close to, but still smaller than RC −1
. In this vective layers form from the boundaries of the domain,
case, the fluid becomes unstable to ODDC, and the sys- rather than from the middle. This is illustrated in Fig-
tem eventually reaches a statistically-stationary ODDC ure 3 where we present results for a simulation with
state. Figure 1 (left) represents a time series of den- P r = τ = 0.1 and R0−1 = 1.75. As usual, the simu-
sity profiles for the simulation with P r = τ = 0.1 and lation starts with a linear density profile at t = 0 and
R0−1 = 4.0 which is illustrative of this regime, i.e. a rapidly becomes unstable to ODDC. This causes weak
strong stratification, but still unstable to ODDC by the deviations from the linear profile observable around
criterion given in equation (2). The vertical axis corre- t ≈ 1000. However, in this case, we can see convective
sponds to the z-coordinate. The simulation starts with layers emerge near the boundaries of the domain around
a linear density profile at t = 0. Throughout the evo- t ≈ 2000. Inspection of R−1 (±75, t) in Figure 3 (right)
−1
lution, we observe weak perturbations from the linear reveals that this time R−1 (z, t) drops below RL near
−1
background, but overall the shape and slope of the den- the domain boundaries first. Meanwhile, R (0, t) re-
−1
sity profiles remain approximately the same. In Figure 1 mains above RL , so the intermediate region remains in
(right), we show the time evolution of R−1 (0, t), namely the state of ODDC, at least temporarily. These bound-
the instantaneous inverse density ratio in the middle of ary convective layers grow in time between t ≈ 2000 and
the domain, as well as R−1 (±75, t), the inverse density t ≈ 6000 through an erosion process that appears to be
ratio closer to the boundaries. We see that both curves consistent (at least qualitatively) with that observed by
start at the initial value R0−1 , increase slightly at early Fuentes et al. (2022) in related simulations. After some
times as a result of ODDC, and later converge asymp- time, the intermediate ODDC region completely disap-
totically to the same constant value as the simulation pears, and is replaced by a thin, very stably stratified
reaches a statistically-stationary state. interface.
The convection on both sides continues to gradually
3.2. High R0−1 erode the remaining interface which disappears around
Next we look at a slightly lower R0−1 . Figure 2 t ≈ 17000. This erosion mechanism will be explored
shows the same diagnostics as Figure 1, for parameters in more detail in a future publication. In the final,
P r = τ = 0.1, R0−1 = 2.5. We see that the domain statistically-stationary stage of the simulation, the do-
also becomes unstable to ODDC at early times, and main is again fully-convective, and has a constant den-
stays in that state for a very long time. In contrast sity except for thin boundary layers.
with Figure 1, however, it does not reach a statistically-
stationary ODDC state. We see that as the ODDC de-
velops, R−1 (0, t) decreases with time, which means that 3.4. Low R0−1
the background density stratification becomes weaker. Finally, we present a simulation at very low R0−1 (but
Eventually around t ≈ 60, 000, we see a clear convec- still R0−1 > 1). We see yet another completely different
tive layer appear in the middle of the domain, char- behavior in the evolution of the system. In Figure 4, we
acterized by a constant density. To understand why present a simulation with P r = τ = 0.1 and R0−1 = 1.25.
a convective layer form, we note that R−1 (0, t) crosses Similarly to other cases, the simulation starts with a lin-
−1
the γ-instability threshold RL , shown as the horizon- ear density profile at t = 0 and evolves into a state of
tal dashed red line in Figure 2 (which is approximately ODDC. This initial phase only lasts until t ≃ 800. After
equal to 1.5 at these parameters, see Mirouh et al. 2012) a short period of time, we observe the formation of five
at t ≈ 60, 000. Shortly thereafter, the γ-instability convective layers in the domain at t ≈ 1000. The con-
causes an inversion in the mean density profile, which vective layers are separated by thin strongly-stratified
6

Figure 1. Left: Time-series of nondimensional mean vertical density profiles (with ρ̄ on the horizontal axis and the z-coordinate
on the vertical axis) for P r = τ = 0.1 and R0−1 = 4.0. Each line represents a mean density profile at a specific time, and each
density profile is horizontally offset (shifted) from earlier time steps for ease of visualization. Lines are plotted every 2500 time
units increments, in 10 different colors. Each blue line corresponds to ∆t = 25, 000 time unit increments. Right: Inverse density
ratios (vertical axis) vs time (horizontal axis): R−1 (0, t) (green), and R−1 (±75, t) (orange). The red horizontal dashed line
−1
shows RL ≃ 1.5, which is the layering threshold for P r = τ = 0.1.

Figure 2. As in Figure 1, for P r = τ = 0.1 and R0−1 = 2.5. Left: Time-series of vertical density profiles. Each blue line
corresponds to ∆t = 12, 500 time unit increments. Right: Inverse density ratios vs time: R−1 (0, t)(green), and R−1 (±75, t)
−1
(orange). The red horizontal dashed line shows RL ≃ 1.5, which is the layering threshold for P r = τ = 0.1, and the black
dashed line corresponds to t = 60, 000 when convection is triggered.

interfaces which are identified by sharp gradients in the t ≈ 3300 and t ≈ 7500, the interface is eroded over time
density profiles. These layers are clearly formed by the by the two convective layers and eventually disappears
−1
γ-instability. Indeed, in this case, R−1 (z, t) < RL ev- at t ≈ 7800, leading once again to a fully convective
erywhere at the start of the simulation, meaning that chemically homogeneous domain.
the γ-instability may drive the formation of convective To summarize, we have found that the long-term evo-
layers everywhere in the domain very early on. Later, lution of this simple model system depends on whether
the number of convective layers is reduced from five to the γ-instability is triggered or not. In the first case,
−1
two, as a result of a combination of interface erosion (e.g. R−1 (z, t) never decreases below RL , and the domain
lowest interface) and merger (e.g. top three interfaces). ultimately achieves a statistically-stationary state of
These processes happen quickly between t ≈ 1000 and ODDC. In all the other cases, the γ-instability drives
t ≈ 3300, resulting in a single interface in the middle of the formation of convective layers as soon as R−1 (z, t)
−1
the domain separating two convection zones. Between drops below the threshold RL somewhere in the do-
 7

Figure 3. As in Figure 1, for P r = τ = 0.1 and R0−1 = 1.75. Left: Time-series of vertical density profiles. Each blue line
corresponds to ∆t = 2000 time unit increments. The various evolutionary stages are marked with arrows. Right: Inverse
−1
density ratios vs time: R−1 (0, t)(green), and R−1 (±75, t) (orange). The red horizontal dashed line shows RL ≃ 1.5, which is
the layering threshold for P r = τ = 0.1, and the black dashed line corresponds to t = 2000 when the convection is triggered.

Figure 4. As in Figure 1, for P r = τ = 0.1 and R0−1 = 1.25. Left: Time-series of vertical density profiles. Each blue line
corresponds to ∆t = 1, 000 time unit increments. The various evolutionary stages are marked with arrows. Right: Inverse
−1
density ratios vs time: R−1 (0, t)(green), and R−1 (±75, t) (orange). The red horizontal dashed line shows RL ≃ 1.5, which is
the layering threshold for P r = τ = 0.1.

main. As in Wood et al. (2013), these convective layers of this conclusion. However, we believe it is quite ro-
grow until they fill the entire domain, either by merging, bust. First, note that while we have presented results for
or by entraining material from the neighboring ODDC P r = τ = 0.1, very similar results are obtained at other
regions. In other words, our numerical experiment has input parameters, in as much as we have observed the
demonstrated that stably stratified ODDC regions can same sequence of possible cases at P r = τ = 0.3, with
−1
only survive for long periods of time provided R−1 (z, t) the same conclusion regarding the role of RL (see Table
−1
never drops below RL . 1). Second, Radko (2007) has demonstrated that a den-
sity staircase (formed of convective layers separated by
4. DISCUSSION & CONCLUSIONS thin stable interfaces) necessarily coarsens with time ei-
Given the idealized nature of the DNS conducted in ther when the density flux increases with layer height or
the previous section, one may question the generality when it decreases with interface stratification. Both of
8

R0−1 Resolution −1
R−1 (z, t) < RL Outcome function of the rotation rate for rapidly-rotating plan-
−1
P r = τ = 0.1 (RL ≃ 1.5) ets. A thorough study of the effects of rotation on the
4.0 64×64×128 no ODDC flux ratio γ −1 is therefore needed. In addition, Fuentes
3.0 64×64×128 no ODDC et al. (2023) showed that rotation can slow down the ero-
2.5 64×64×128 yes, locally CCL sion rate of stably-stratified region adjacent to convec-
2.0 64×64×128 yes, locally BCL tive ones. This could slow down, but cannot completely
1.75 128×128×256 yes, locally BCL halt, the growth of the fully convective zones once they
1.4 128×128×256 yes, locally BCL have been created.
1.25 128×128×256 yes γ-CL Of course, one should not directly compare the de-
P r = τ = 0.3 (RL−1
≃ 1.25) tailed evolution of the model system studied in this pa-
1.75 64×64×128 no ODDC per to the evolution of a planet’s diffuse core, because
1.55 64×64×128 yes, locally CCL
the boundary conditions we have used impose fixed and
equal fluxes of heat and composition on both boundaries,
1.4 128×128×256 yes, locally BCL
while in reality this would certainly not be the case. In
1.3 128×128×256 yes, locally BCL
particular, depending on where we assume these bound-
1.25 128×128×256 yes, locally BCL
aries are, the top and bottom fluxes are likely to differ
1.15 128×128×256 yes, locally BCL
from one another and evolve with time on the planetary
1.10 128×128×256 yes γ-CL
evolution timescale (Hindman & Fuentes 2023). How-
Table 1. List of simulations for P r = τ = 0.1 and ever, we have demonstrated that the development of the
P r = τ = 0.3, and their outcomes. The first column shows γ-instability is a local process, which is independent of
the initial inverse density ratio. The second column shows what happens at the boundaries, and only depends on
the adopted resolution in terms of numbers of modes, see sec- the local stratification of that core at any point in time.
tion 2. The third column shows whether R−1 (z, t) ever drops Because of this, the following conclusions hold:
−1
below the threshold, RL , at any point in the simulation. Fi-
nally, the fourth column indicates the outcome: ”ODDC” is 1. We confirm the result of Wood et al. (2013) and
a statistically stationary ODDC state, ‘CCL’ means that a Mirouh et al. (2012) that the γ-instability plays
layer forms in the center of the domain first, ‘BCL’ means
a fundamental role in the dynamics of stably-
that layers form near the boundaries first, ‘γ-CL’indicates a
simulation with a multitude of layers that form simultane- stratified regions with a super-adiabatic temper-
ously. ature gradient. More specifically, we have found
that the γ-instability triggers convective layer for-
mation locally as soon as the local inverse den-
these conditions are satisfied in ODDC staircases formed sity ratio, R−1 (z, t) (see eq. 27), drops be-
−1
by the γ-instability (Wood et al. 2013; Moll et al. 2017), low RL (P r, τ ) (see eq. 4). Only regions with
−1 −1
suggesting that the coarsening process observed in figure R (z, t) > RL , ∀z, t can remain in a long-term
4 (and more generally in any weakly-stratified region) is state of ODDC. This, incidentally, shows the im-
universal and will not stop until the entire region is fully portance of having good estimates for both P r and
convective. The robustness of cases illustrated in figures τ (cf. French et al. 2019; Preising et al. 2023) as
−1
2 and 3, where a convective layer grows by eroding the well as as a good model for RL (P r, τ ) (Mirouh
−1
nearby ODDC layer(s), is somewhat less clear but re- et al. 2012). If the planet is rapidly rotating, RL
mains likely because the entrainment process takes place may also depend on the rotation rate.
on dynamical timescales, while the restratification of the
2. Our results, combined with those of Radko (2007),
ODDC layer (necessary to preserve it) would take place
challenge the standard paradigm of composition-
on much longer evolutionary timescales. Finally, giant
ally stratified regions in giant planets which pic-
planets in our solar system are rapid rotators, and so
tures them as a stack of long-lived, well-mixed,
one may wonder how rotation affects the γ-instability
shallow convective layers, separated by thin dif-
and the evolution of convective layers. Moll & Garaud
fusive interfaces (Stevenson 1982). This has not
(2017) showed that the γ-instability is not directly af-
been seen as an outcome of our simulations. At the
fected by rotation, but the flux ratio γ −1 itself, and its
moderate stratifications where they spontaneously
dependence on R−1 , may both depend on the planet’s
form, convective layers always seem to grow and
rotation rate. As such, while our general conclusion con-
merge, until the convection zone occupies the en-
tinues to hold in the presence of rotation, the critical
−1 tire domain. These processes happen rapidly and
inverse density ratio for layer formation RL could be a
irreversibly. We believe that rapid rotation might
 9

slow down, but cannot completely hinder the pro- the satellite Io (Idini & Stevenson 2022). Alternatively,
cess. perhaps the higher initial thermal energy of the more
massive Jupiter, and corresponding higher initial central
3. The fact that Saturn is observed to have temperature, have put the planet in a different regime,
an extended stably-stratified core suggests that having triggered convection sometimes in the past so
R−1 (z, t) has never decreased below the thresh- that the dilute core today is fully or partially homoge-
−1
old RL in that core throughout the evolution of nized.
the planet. This might help constrain formation Further analysis is required to translate the conditions
and evolution models of Saturn. of instability presented here to the conditions expected
in the interiors of Jupiter and Saturn. In the giant plan-
There are a number of further investigations in planetary ets, such work must consider the nontrivial calculation
modeling that should be explored, given these results. of how material properties change with depth, including
We have strong evidence that Saturn’s dilute core to- the adiabatic or super-adiabatic behavior of the mix-
day is stably stratified (Fuller 2014; Mankovich & Fuller ture of H-He fluid with heavy elements and the uncer-
2021; Fortney et al. 2023). Our work here suggests that tainties therein. Future research will aim to place new
this stable stratification must have already been estab- constraints on the range of core temperatures, strength
lished in Saturn during the planet’s formation. This of compositional gradients, and radial extent of dilute
places constraints on the central temperature (and the cores in Jupiter and Saturn, to assess what present and
corresponding temperature gradient in the dilute core) past interior structures are possible given the mixing ar-
at that time, as well as at any time between formation guments developed here.
until today, to ensure that the γ-instability is never trig-
A.T. and P.G. gratefully acknowledge help from the
gered. These temperature constraints directly inform
’Dedalus’ development team and support from National
static structure models today as well as thermal evolu-
Science Foundation grant AST-1908338. B.I. is sup-
tion models that aim to reproduce the planet’s measured
ported by the University of California President’s Post-
luminosity at its current age.
doctoral Fellowship Program. Simulations were per-
Without seismic constraints for Jupiter, we cannot as
formed using the Extreme Science and Engineering Dis-
easily constrain the extent and thermal state of its di-
covery Environment (XSEDE) Expanse supercomputer
lute core. One could imagine that like Saturn, Jupiter’s
at the San Diego Supercomputer Center through alloca-
dilute core is also stably stratified. This scenario is
tion TG-PHY210050, as well as the Lux supercomputer
supported by the suggestion that internal gravity waves
at UC Santa Cruz, funded by NSF MRI grant AST-
trapped in such a core produce enough of a gravitational
1828315.
signature to perturb the gravity field of tides raised by

REFERENCES
Baines, P. G., & Gill, A. E. 1969, Journal of Fluid Fortney, J. J., Militzer, B., Mankovich, C. R., et al. 2023,
Mechanics, 37, 289–306, doi: 10.1017/S0022112069000553 Saturn’s Interior After the Cassini Grand Finale,
Bolton, S. J., Lunine, J., Stevenson, D., et al. 2017, Space doi: https://ui.adsabs.harvard.edu/link gateway/
Science Reviews, 213, 5, 2023arXiv230409215F/doi:10.48550/arXiv.2304.09215
doi: https://doi.org/10.1007/s11214-017-0429-6 French, R. G., McGhee-French, C. A., Nicholson, P. D., &
Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D., & Hedman, M. M. 2019, Icarus, 319, 599,
Brown, B. P. 2020, Physical Review Research, 2, 023068,
doi: https://doi.org/10.1016/j.icarus.2018.10.013
doi: 10.1103/PhysRevResearch.2.023068
Fuentes, J. R., Anders, E. H., Cumming, A., & Hindman,
Debras, F., & Chabrier, G. 2019, The Astrophysical
B. W. 2023, The Astrophysical Journal Letters, 950, L4,
Journal, 872, 100, doi: 10.3847/1538-4357/aaff65
doi: 10.3847/2041-8213/acd774
Durante, D., Parisi, M., Serra, D., et al. 2020, Geophysical
Research Letters, 47, e2019GL086572, Fuentes, J. R., Cumming, A., & Anders, E. H. 2022,
doi: https://doi.org/10.1029/2019GL086572 Physical Review Fluids, 7, 124501,
Folkner, W., Iess, L., Anderson, J., et al. 2017, Geophysical doi: 10.1103/PhysRevFluids.7.124501
Research Letters, 44, 4694, Fuller, J. 2014, icarus, 242, 283,
doi: https://doi.org/10.1002/2017GL073140 doi: https://doi.org/10.1016/j.icarus.2014.08.006
10

Guillot, T. 2005, Annu. Rev. Earth Planet. Sci., 33, 493, Mizuno, H., Nakazawa, K., & Hayashi, C. 1978, Progress of
doi: 10.1146/annurev.earth.32.101802.120325 Theoretical Physics, 60, 699, doi: 10.1143/PTP.60.699
Hedman, M., & Nicholson, P. 2013, The Astronomical Moll, R., & Garaud, P. 2017, Astrophys. J., 834, 44,
Journal, 146, 12, doi: 10.1088/0004-6256/146/1/12 doi: 10.3847/1538-4357/834/1/44
Moll, R., Garaud, P., Mankovich, C., & Fortney, J. J. 2017,
—. 2014, Monthly Notices of the Royal Astronomical
ApJ, 849, 24, doi: 10.3847/1538-4357/aa8d74
Society, 444, 1369, doi: 10.1093/mnras/stu1503
Nettelmann, N., Fortney, J. J., Moore, K., & Mankovich,
Hedman, M., Nicholson, P., & French, R. 2018, The
C. 2015, MNRAS, 447, 3422, doi: 10.1093/mnras/stu2634
Astronomical Journal, 157, 18, Perri, F., & Cameron, A. G. 1974, Icarus, 22, 416,
doi: 10.3847/1538-3881/aaf0a6 doi: https://doi.org/10.1016/0019-1035(74)90074-8
Hindman, B. W., & Fuentes, J. R. 2023, Dwindling Surface Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996,
Cooling of a Rotating Jovian Planet Leads to a icarus, 124, 62,
Convection Zone that Grows to a Finite Depth, doi: https://doi.org/10.1006/icar.1996.0190
doi: 10.3847/2041-8213/ad0642 Preising, M., French, M., Mankovich, C., Soubiran, F., &
Howard, S., Guillot, T., Bazot, M., et al. 2023, arXiv Redmer, R. 2023, ApJS, 269, 47,
preprint arXiv:2302.09082, doi: 10.3847/1538-4365/ad0293
doi: 10.1051/0004-6361/202245625 Radko, T. 2003, J. Fluid Mech., 497, 365,
doi: 10.1017/S0022112003006785
Idini, B., & Stevenson, D. J. 2022, The Planetary Science
Radko, T. 2007, Journal of Fluid Mechanics, 577, 251,
Journal, 3, 89, doi: 10.3847/PSJ/ac6179
doi: 10.1017/S0022112007004703
Iess, L., Militzer, B., Kaspi, Y., et al. 2019, Science, 364,
Radko, T. 2013, Double-diffusive convection (Cambridge
eaat2965, doi: 10.1126/science.aat2965 University Press),
Kato, S. 1966, Proc. Astro. Soc. Japan, 18, 374 doi: https://doi.org/10.1017/CBO9781139034173
Mankovich, C. R. 2020, AGU Advances, 1, e2019AV000142, Rosenblum, E., Garaud, P., Traxler, A., & Stellmach, S.
doi: https://doi.org/10.1029/2019AV000142 2011, Astrophys. J., 731, 66,
Mankovich, C. R., & Fuller, J. 2021, Nature Astronomy, 5, doi: 10.1088/0004-637X/731/1/66
1103, doi: 10.1038/s41550-021-01448-3 Safronov, V. 1969, NASA Tech. Transl. F-677; Moscow,
Nauka
Marley, M. S., & Porco, C. C. 1993, Icarus, 106, 508,
Spiegel, E. A., & Veronis, G. 1960, Astrophys. J., 131, 442,
doi: https://doi.org/10.1006/icar.1993.1189
doi: https://ui.adsabs.harvard.edu/link gateway/
Miguel, Y., Bazot, M., Guillot, T., et al. 2022, A&A, 662,
1960ApJ...131..442S/doi:10.1086/146849
A18, doi: 10.1051/0004-6361/202243207
Stevenson, D. J. 1982, Planet. Space Sci., 30, 755,
Militzer, B., Hubbard, W. B., Wahl, S., et al. 2022, The doi: 10.1016/0032-0633(82)90108-8
planetary science journal, 3, 185, Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017,
doi: 10.3847/PSJ/ac7ec8 Geophysical Research Letters, 44, 4649,
Mirouh, G. M., Garaud, P., Stellmach, S., Traxler, A. L., & doi: https://doi.org/10.1002/2017GL073160
Wood, T. S. 2012, Astrophys. J., 750, 61, Walin, G. 1964, Tellus, 16, 389,
doi: 10.1088/0004-637X/750/1/61 doi: 10.3402/tellusa.v16i3.8930
Mizuno, H. 1980, Progress of Theoretical Physics, 64, 544, Wood, T. S., Garaud, P., & Stellmach, S. 2013, Astrophys.
J., 768, 157, doi: 10.1088/0004-637X/768/2/157
doi: 10.1143/PTP.64.544