Relation of Gravity, Winds, and the Moment of Inertia of Jupiter and Saturn

We study the relationship of zonal gravity coefficients, J_2n, zonal winds, and axial moment of inertia (MoI) by constructing models for the interiors of giant planets. We employ the nonperturbative concentric Maclaurin spheroid (CMS) method to construct both physical (realistic equation of state and barotropes) and abstract (small number of constant-density spheroids) interior models. We find that accurate gravity measurements of Jupiter's and Saturn's J_2, J_4, and J_6 by Juno and Cassini spacecrafts do not uniquely determine the MoI of either planet but do constrain it to better than 1%. Zonal winds (or differential rotation, DR) then emerge as the leading source of uncertainty. For Saturn, they are predicted to decrease the MoI by 0.4% because they reach a depth of ~9000 km while on Jupiter, they appear to reach only ~3000 km. We thus predict DR to affect Jupiter's MoI by only 0.01%, too small by one order of magnitude to be detectable by the Juno spacecraft. We find winds primarily affect the MoI indirectly via the gravity harmonic J_6 while direct contributions are much smaller because the effects of pro- and retrograde winds cancel. DR contributes +6% and -0.8% to Saturn's and Jupiter's J_6 value, respectively. This changes the J_6 contribution that comes from the uniformly rotating bulk of the planet that correlates most strongly with the predicted MoI. With our physical models, we predict Jupiter's MoI to be 0.26393+-0.00001. For Saturn, we predict 0.2181+-0.0002, assuming a rotation period of 10:33:34 h that matches the observed polar radius.


INTRODUCTION
The angular momentum of a giant planet must be accurately known to calculate the planet's precession rate, which is the crucial quantity to determine whether it is in a spin-orbit resonance. Such resonances have been invoked, along with additional assumptions to explain the obliquities of Saturn, 27 • (Saillenfest et al. 2021;Wisdom et al. 2022), Jupiter, 3 • (Ward & Canup 2006), and Uranus, 98 • (Saillenfest et al. 2022). The planetary spin angular momentum contributes 99% of the angular momenta of the Jovian or Saturnian systems, the rest coming from the most massive satellites. To high order, the total angular momentum of a planetary system is conserved over billions of years while the planet's moment of inertia C changes (Helled 2012;Nettelmann et al. 2012a) due to secular cooling and other processes like helium rain (Wilson & Militzer 2012a) and core erosion (Wilson & Militzer 2012b), and satellite orbits exchange angular momentum with the planet through tidal interactions (Fuller et al. 2016).
The space missions Juno (Bolton et al. 2017) and Cassini (Spilker 2019) have provided us with a wealth of new data for Jupiter and Saturn. Multiple close flybys have mapped the gravity fields of these planets with a high level of precision (Durante et al. 2020;Iess et al. 2019) that far exceeds the earlier measurements by the Pioneer and Voyager missions (Campbell & Synnott 1985;Campbell & Anderson 1989). The new measurements have also led to a revision of the assumptions that are employed when interior models are constructed. Traditionally, the interiors of Saturn and Jupiter were represented by three layer models (Guillot et al. 2004;Saumon & Guillot 2004a;Nettelmann et al. 2012b;Hubbard & Militzer 2016a) that start with an outer layer that is predominantly composed of molecular hydrogen, a deeper layer where hydrogen is metallic, and compact core that is composed of up to 100% of elements heavier than helium. There was sufficient flexibility in choosing the layer thicknesses and the mass fractions of helium, Y , and heavier elements, Z, to match the earlier spacecraft measurements.
Still, the predictions from various types of three layer models were not always found to be in perfect agreement for two main reasons. Early interior models relied on the theory of figures (ToF) (Zharkov & Trubitsyn 1978), a perturbative approach, to capture the gravitational and rotational effects in a planet's interior. Most calculations employed the third and fourth order version of the ToF but this technique has recently been extended to seventh order (Nettelmann et al. 2021). With the development of the concentric Maclaurin spheroid method (CMS), it became possible to construct giant planet interior models nonperturbatively (Hubbard 2013).
The second source of uncertainty is the equation of state (EOS) of hydrogen-helium mixtures at high pressure (Vorberger et al. 2007; Morales et al. 2010). While shock wave measurements (Zeldovich & Raizer 1968) now routinely reach the relevant regime of megabar pressures (Da Silva et al. 1997;Collins et al. 1998;Knudson et al. 2001;Celliers et al. 2010), the temperatures in these experiments are much higher than those in giant planet interiors (Militzer et al. 2016). For this reason, interior models invoke theoretical methods (Saumon et al. 1995a) and ab initio simulations (Militzer et al. 2008;Nettelmann et al. 2008) to construct an EOS for hydrogen-helium mixtures and then add heavy elements within the linear mixing approximation (Soubiran & Militzer 2015;Ni 2018). A direct experimental confirmation of the prediction from ab initio simulations of hydrogen-helium mixtures under giant planet interior conditions would be highly valuable even though simulation results for other materials were found to be in good agreement with shock experiments (French et al. 2009;Millot et al. 2020).
For Jupiter, the Juno spacecraft obtained smaller magnitudes for the harmonics J 4 and J 6 than interior models had predicted (Hubbard & Militzer 2016a). Matching and interpreting these measurements has led authors to introduce a number of novel assumptions into interior models. One can adopt a 1-bar temperature that is higher (Wahl et al. 2017b;Miguel et al. 2022) than the Galileo value of 166.1 K or invoke a less-than-protosolar abundance of heavy Z elements (Hubbard & Militzer 2016a;Nettelmann 2017;Wahl et al. 2017b). Both modifications reduce the density of the molecular outer layer, which makes it easier to match J 4 and J 6 . Wahl et al. (2017b) introduced the concept of a dilute core, which partially addressed the J 4 -J 6 challenge. Debras & Chabrier (2019) adopted the dilute core concept and then decreased the heavy Z element fraction at an intermediate layer. Most recently Militzer et al. (2022) matched all observed J n values exactly by simultaneously optimizing parameters of the dilute core and models for the zonal winds.
The high-precision values from the Juno and Cassini missions for Jupiter's and Saturn's zonal gravitational harmonics, J n , provide important constraints on the interior mass distributions and thereby also constrain the moment of inertia as we will demonstrate in this article. A different constraint, the value of the spin angular momentum, J , comes from measurement of forced precession of the planet's rotation axis. As the precession periods are very long, respectively ∼0.5×10 6 years for Jupiter and ∼2×10 6 years for Saturn (Ward & Canup 2006), high-precision pole-position measurements over a long time baseline are necessary to measure J to better than 1%. In principle, if the planet rotates uniformly and its spin rate, ω, is known, one can obtain the axial moment of inertia, C, via C = J /ω, which would provide an independent constraint on the interior mass distribution.
For convenience, a planet's momentum of inertia is typically reported in normalized form, MoI≡ C/M R 2 e . While we normalize by the planet's mass, M , and the present-day equatorial radius, R e at a pressure of 1 bar, one should note that other authors have used the volumetric radius (Ni 2018) or made the radius age-dependent (Helled 2012). In this paper, we systematically investigate how much MoI can vary for models which have exact fixed values of ω and zonal gravitational harmonics J n up to some limiting degree n, and thus illustrate the role of MoI as an independent constraint. Note that the approximate Radau-Darwin formula (e.g., Bourda & Capitaine (2004)), which posits a unique relation between J 2 , ω, and MoI, is too inaccurate to be relevant to this investigation because Jupiter and Saturn rotate rapidly and the density varies significantly throughout their interiors (Wahl et al. 2021). When we construct models for giant planet interiors, we assume hydrostatic equilibrium and that the density increases with pressure. Since this concentrates mass in the planet's center, we expect the inferred MoI to be substantially less than 2/5, the value for a single constant-density Maclaurin spheroid independent of its rotation rate.
The article is organized as follows. In Sec. 2, we show how a planet's moment of inertia and angular momentum are calculated with the CMS method. We introduce physical and abstract models for giant planet interiors. We also explain that differential rotation (DR) in a planet has direct and indirect effects. The direct effect is introduced when the observed zonal winds, that move at different angular velocities, are projected into the interior and thereby cause a planet's angular momentum to deviate from the value of an object that rotates uniformly. However, the zonal winds also make dynamical contributions to a planet's gravitational harmonics. They thereby reduce the static contributions slightly that come from the mass distribution in the planet's interior when models are constructed to match a spacecraft's gravity measurements. This change in the mass distribution also affects the resulting moment of inertia, which we call the indirect effect of DR.
In Sec. 3, we first discuss our predictions for Saturn's momentum of inertia and illustrate how sensitively it depends on the gravity harmonics J 4 and J 6 . We find that the dynamical contributions to J 6 play a critical role. Then we derive the angular momentum for arbitrary giant planets, for which the mass, equatorial radius, J 2 , and rotational period have been measured. We present results from different models for Jupiter's interior, which includes CMS calculations that we performed for Jupiter models of other authors. Finally we compare our momentum of inertia values with earlier predictions in the literature before we conclude in Sec. 4.

METHODS
The normalized moment of inertia of an axially symmetric body can be derived from this integral over all fluid parcels as function of radius and µ = cos(θ) with θ being the colatitude, where l = r 1 − µ 2 is the distance from the axis of rotation and R(µ) marks the outer boundary of the planet. In the CMS method, one represents the mass in the planet's interior by a series of nested constant-density spheroids each adding a small density contribution, δ j , that lets the combined density increase with depth. After carrying out the radius integration, the MoI can be written as a sum over spheroids, where r j (µ) marks the outer boundary of the spheroid with index j. For a uniformly rotating (UR) body, the normalized spin angular momentum is given by J UR norm = √ q rot C/M R 2 e with q rot being the dimensionless rotational parameter, that compares the magnitudes of the centrifugal and gravitational potentials. If the body is rotating differentially, one needs to revert to the 2D integral, where v(r, µ) is the fluid velocity andv is that of the uniformly rotating background,v = l * ω. For convenience, one may choose to define an effective or average moment of inertia for a differentially rotating body, and compare it with predictions of Eq. 2. We call this difference the direct effect of DR on the predicted MoI, to be compared with the indirect effect that emerged because DR affects the interior density structure and thus the calculated gravity harmonics, in particular J 6 . In Tab. 1, we quantify the indirect DR effect by comparing the MoI values, C (DR) and C (UR) , derived from Eq. 2 for a model that invokes DR effects and a model that does not when they both match the observed J n . We find that the magnitude of the direct DR effect is much smaller than the indirect one (Tab. 1) because contributions from pro-and retrograde jets to the direct effect partially cancel. Direct DR effects increase Jupiter's MoI by 0.0015% because the prograde winds in the equatorial region dominate. For Saturn, we find that the retrograde winds at a latitude of ∼35 • dominate over the prograde equatorial jet, which implies that direct DR effects lower the planet's angular momentum by −0.13%.

CMS Technique
The spheroid surfaces r j (µ) are contours of constant pressure, temperature, composition, and potential. The potential combines centrifugal and gravitational contributions, Q + V . According to Zharkov & Trubitsyn (1978), the gravitational potential can be expanded in the following form, where P n are the Legendre polynomials of order n and the J n are the gravity harmonics given by According to Hubbard (2013), the gravitational potential V j of a point (r j , µ) on spheroid j is decomposed into contributions from interior spheroids (j = i . . . N − 1), and exterior spheroids (j = 0 . . . i − 1), Following the derivation in Hubbard (2013), we define the interior harmonics and the exterior harmonics with a special case for n = 2, and finally, where M is the total mass of the planet. One should note that during the numerical evaluation of these expressions, it is recommended to work with harmonics that have been renormalized by the powers of the equatorial spheroid radii, λ i . These equatorial points (r j = λ j , µ = 0) serve as anchors for all spheroid shapes. This is where the reference value of the potential is computed that one uses to adjust the spheroid shape until a self-consistent solution emerges for which all spheroids are equipotential surfaces. It is important to choose the λ i grid points wisely in order to minimize the discretization error that is inherent to the CMS approach. We recommended choosing them so that a logarithmic grid in density emerges, ρ(λ i+1 )/ρ(λ i )=constant ). This grid choice allows us to obtain converged results when we construct our physical models with N S = 2048 spheroids.
In addition to gravity, one needs to consider the centrifugal potential, which takes the following simple form for a uniformly rotating body, Q(l) = 1 2 l 2 ω 2 . We employ this formula when we construct models for Jupiter's interior and then introduce DR effects by solving the thermal wind equation (Kaspi et al. 2016) to derive the density perturbation, for a rotating, oblate planet (Cao & Stevenson 2017) in geostrophic balance. z is the vertical coordinate that is parallel to the axis of rotation. ρ is static background density that we derive with the CMS method. u is the differential flow velocity with respect to the uniform rotation rate, ω. g is the acceleration that we derive from the gravitationalcentrifugal potential, V + Q, in our CMS calculations. s is the distance from the equatorial plane along a path on an equipotential. We represent the flow field u as a product of the surface winds, u s , from Tollefson et al. (2017) and a decay function of sin 2 (x) form from Militzer et al. (2019). This function facilitates a rather sharp drop similar to functions employed in Galanti & Kaspi (2020) and Dietrich et al. (2021).
Since the winds on Saturn reach much deeper, we treat them nonperturbatively by introducing DR on cylinders directly into the CMS calculations by modifying the centrifugal potential, Since we assume potential theory, a cylinder's angular velocity, ω(l), cannot decay with depth, which means we are only able to include the prograde equatorial jet and first retrograde jet at ∼35 • that were characterized by tracking the cloud motion in Saturn's visible atmosphere (Sanchez-Lavega et al. 2000;García-Melendo et al. 2011).

Physical Interior Models
In Fig. 1, we illustrate our physical interior models for Jupiter and Saturn. Since planets cool by convection, we assume most layers in their interiors are isentropic and of constant composition. We represent their outer envelope where hydrogen is molecular by the parameters (S mol ,Ỹ mol , Z mol ) for entropy, helium mass fraction and the fraction of heavy elements. We defineỸ ≡ Y /(X + Y ) with X and Y being the mass fractions of hydrogen and helium so that X + Y + Z = 1. We require Z mol to be at least protosolar, Z PS = 1.53% (Lodders 2010). The entropy is chosen to match the temperature at 1 bar: 142.7 K for Saturn (Lindal et al. 1981) and 166.1 K for Jupiter (Seiff et al. 1997) that was measured in situ by the Galileo entry probe. For Jupiter, we also consider an alternate, slightly higher temperature of 170 K from a recent reassessment of the Voyager radio occultation measurements (Gupta et al. 2022).
To construct EOSs for models in this article, we start from the ab initio EOS that Militzer & Hubbard (2013) computed for one hydrogen-helium mixing ratio. With these calculations, absolute entropies (Militzer 2013) were derived that implicitly set the temperature profiles in our models. We use our helium EOS from Militzer (2006Militzer ( , 2009 to perturb helium fraction in our H-He EOS as we detailed in Hubbard & Militzer (2016b). We also follow this article when we introduce heavily elements into our models. Their detailed composition is not important as long as they are substantially more dense than hydrogen and helium. Ice, rocky materials and iron are all sufficiently dense so that they add mass but do not increase the volume of the mixture too much. At low pressure where the ab initio simulations do not work, we revert back to the semi-analytical EOS by Saumon et al. (1995b).
When hydrogen assumes an atomic/metallic state at approximately 80-100 GPa (Morales et al. 2009), helium remains an insulator and the two fluids are predicted to become immiscible (Stevenson & Salpeter 1977;Brygoo et al. 2021). There is indeed good evidence that helium rain has occurred in Jupiter because the Galileo entry probe measured a helium mass fraction ofỸ = 0.238 ± 0.005 (von Zahn et al. 1998) that is well below the protosolar value of 0.2777 (Lodders 2010). Furthermore, neon in Jupiter's atmosphere was measured to be nine-fold depleted relative to solar, and this can be attributed to efficient dissolution in helium droplets (Roulston & Stevenson 1995;Wilson & Militzer 2010). So for our Jupiter models, we adopt the value from the Galileo entry probe forỸ mol and for Saturn, we make it a free parameter but constrain it to be no higher than the protosolar value because we have no information on how much helium rain has occurred in this planet.
For both planets, we chose values for the beginning and ending pressures of the helium rain layer that are compatible with the immiscibility region that Morales et al. (2013) derived with ab initio computer simulations (see Militzer et al. (2019) for details). Across this layer, we assume (S,Ỹ , Z) vary gradually with increasing pressure until they reach the values of the metallic layer (S met ,Ỹ met , Z met ) where they are again constant since we assume this layer to be homogeneous and convective.Ỹ met is adjusted iteratively so that the planet as a whole assumes a protosolar helium abundance. This also assuresỸ met >Ỹ mol . We prevent the heavy element abundance from decreasing with depth, Jupiter Saturn Figure 1. Models of Jupiter and Saturn based on CMS calculations that match the gravity measurements of both planets. The zonal winds on Saturn are predicted to reach a depth of ∼9000 km, involving ∼7% of the planet's mass. On Jupiter, they are predicted to reach ∼3000 km and thus involve only 1% of the planet's mass. Because Jupiter is more massive, the pressure rises more rapidly with depth. Therefore the helium rain layer, predicted to start approximately at 80-100 GPa, is located closer to the surface. While the gravity measurements for Jupiter imply that the planet has a dilute core, the state of Saturn's core is less certain. Here we show a model with a compact core constructed to match the gravity measurements.
Z met ≥ Z mol . Every layer is either homogeneous and convective or Ledoux stable (Ledoux 1947). This sets our models apart from those constructed by Debras et al. (2021) who introduced a layer where Z decreases with depth in order to match Jupiter's J 4 . Instead our Jupiter models all have a dilute core with Z ≈ 0.18 (see Fig. 1) because this key restriction allows us to match the entire set of gravity measurements of the Juno spacecraft under one set of physical assumptions ). For our Monte Carlo calculations of Jupiter's interior (Militzer 2023), we vary the beginning and end pressure of the helium rain layer but apply constraints so that they remain compatible with H-He phase diagram as derived by Morales et al. (2009). We also vary a parameter α that controls the shape of the helium profile in this layer, as we explain in Militzer et al. (2022). During the Monte Carlo calculations, we also vary the beginning and end pressure of the core transition layer, which we assume to be stably stratified since the abundance of heavy elements increases from Z mol to Z met in this layer. We also allow Z mol and Z met to vary as long as they meet the constraint we discussed in the previous paragraph. More details of our Monte Carlo approach are given in Militzer et al. (2022).
For our Saturn models, we assume a traditional compact core that is composed up to 100% of heavy elements because this assumption was sufficient to match the gravity measurements by the spacecraft (Iess et al. 2019), but there are alternate core models constructed to match ring seismological data (Mankovich & Fuller 2021).

Abstract N Spheroid Models
In the previous section, we described physical interior models in hydrostatic equilibrium that rely on a realistic EOS for H-He mixtures. To explore more general behavior, we now investigate simplified models with N S spheroids. We  (4) 0.078826 ± 0.000003 0.08655 ± 0.00008 Direct DR effect, ( still require each spheroid surface to be an equipotential but spheroid densities, ρ i , are arbitrary as long as the densities monotonically increase toward the planet's interior, ρ i+1 ≥ ρ i . We can set the density of the outermost spheroid to zero, ρ 0 = 0, because in realistic interior models, the density of the outermost layer is typically much lower than that of deeper layers. (We also construct models in which ρ 0 is a free parameter, but they behave similarly, and in the limit of large N S , the difference becomes negligible.) We initialize the equatorial radii of all spheroids, starting from i = 0 . . . N S − 1, to fall in a linear grid, λ i = 1 − i/N S . While we keep the outermost spheroid anchored at λ 0 = 1, we repeatedly scale all interior λ i>0 points uniformly to obtain a model that matches the planet's mass and J 2 exactly. We add a penalty term to the Monte Carlo (MC) cost function if λ i > λ 0 .
Since matching M and J 2 requires two free parameters, we also scale all density values, ρ i , uniformly. So after every update of the spheroid shapes, we employ a Newton-Raphson step to scale ρ i and λ i grids simultaneously. We also institute a maximum density of 10 PU (planetary unit of density, M/R 3 e ) to prevent pathological situations in which the radius of the innermost spheroid becomes very small while its density becomes extremely large. Movshovitz et al. (2020) and Neuenschwander et al. (2021) also introduced upper limits on density. We consider 10 PU to be a reasonable choice because for Jupiter, it corresponds to a density of 52 g cm −3 , which exceeds the density of iron that is ∼27 g cm −3 at Jupiter's core conditions (Wilson & Militzer 2014). The described set of assumptions lead to a stable procedure with N S -1 free input parameters (ρ i>0 ) that is amendable for MC sampling.
Since we do not employ a physical EOS or make specific assumptions about the planet's composition or temperature profile, our abstract models share similarities with the empirical models by Helled et al. (2009) andNeuenschwander et al. (2021) or the composition-free models by Movshovitz et al. (2020) who represented the Saturn interior density profile by three quadratic functions before conducting MC calculations to match the Cassini gravity measurements.

Saturn
In Fig. 2, we show MoI values computed for the physical models of Saturn's interior in Fig. 1, as well as for the abstract N S spheroid models. The dominant source of uncertainty in the computed MoI is the planet's period of rotation, which cannot be derived from the planet's virtually axisymmetric magnetic field. This is not the case for Jupiter, whose rotation period is known to a fraction of a second (see Tab. 1). Without any constraints on the rotation period, the predictions for Saturn's MoI vary by ∼2%. Still all values predict that Saturn is not currently in a spin-orbit resonance with Neptune today (Wisdom et al. 2022). For all rotation periods shown in Fig. 2, we can construct interior models that match the entire set of gravity coefficients that the Cassini spacecraft measured during its ultimate set of orbits (Iess et al. 2019), so gravity measurements alone are insufficient to constrain the rotation period. Only if we match the planet's polar radius as measured by the Voyager spacecraft using radio occultation, the now-preferred period of 10:33:34 h ± 55 s emerges . This rotation period is in remarkably good agreement with the value of 10:33:38 h +112s −89s inferred from waves observed in Saturn's rings (Mankovich et al. 2019 Once a rotation period has been selected, the remaining uncertainty is dominated by effects of differential rotation (DR), which amount to about 0.4%. Without DR effects, we are only able to match the gravity harmonics J 2 -J 6 , and already matching J 6 requires us to introduce one additional adjustable parameter, so we add an artificial density jump (Iess et al. 2019). The comparison of predictions from model with and without DR in Fig. 4 illustrates that DR effects are much more important for Saturn than for Jupiter. When we include DR effects in our Saturn models, we are able to match the entire set of gravity harmonics J 2 -J 12 without an artificial density jump. We find that resulting MoI drops 0.4% below predictions from models that match J 2 -J 6 without DR.
To better understand this drop, we constructed MC ensembles of abstract models of Saturn's interior that match q rot and J 2 without invoking DR. In Fig. 3c, we plot the posterior distribution of the computed MoI in J 4 -J 6 space. We also show the Cassini measurements (Iess et al. 2019) and the model from Fig. 4 without DR nor artificial density jump, matching the observed J 2 and J 4 . We estimate DR effects increase Saturn's J 6 from ∼81 × 10 −6 to the observed value of 86.340 ×10 −6 . Fig. 3c shows that the Cassini measurements place Saturn in a regime where an increase in J 6 (or in J 4 ) leads to an increase in the MoI: ∂MoI ∂J int 6 > 0. At the same time, models without DR in Fig. 2 predict a larger MoI than models with DR. This lets us conclude that when models with DR are constructed to match the Cassini measurements, DR effects reduce the contribution to J 6 that comes from the uniformly rotating bulk of the interior, J int 6 . So when models with and without DR are    (Zharkov & Trubitsyn 1978). The variation of the ratios with the order n is due to the contribution of higher-order terms in the ToF, captured to high precision in the nonperturbative CMS method. The figure also illustrates that the effects of differential rotation are much stronger for Saturn (blue shaded area) than for Jupiter because Saturn's winds extend to a greater depth of ∼9000 km (Iess et al. 2018). Saturn models with differential rotation fit the observed moments up to J12 while uniform rotation models can only fit coefficients up to J4, or up to J6 if a density jump is included. Here we compare the Juno measurements with a uniform rotation model for Jupiter's interior while the model with differential rotation in Militzer et al. (2022) matches the entire set of gravity coefficients. compared, both matching the spacecraft data, models with DR predict a smaller MoI because their J int 6 is reduced by contributions to J 6 from DR. It is primarily this change to the J 6 term that affects the MoI while the DR contributions to J 2 and J 4 are too small to matter. On the other hand, DR effects dominate the higher order J n starting with J 8 (see Fig. 4) but their values are controlled by the outer layers of the planet (Guillot 2005;Nettelmann et al. 2013;Fortney et al. 2016;Militzer et al. 2016) where the density is comparatively low, and therefore they do not contribute much to the MoI. We conclude that DR effects couple to the MoI mostly via J 6 .
While the models in Fig. 3 only match J 2 , we compare MC ensembles of Saturn models in Fig. 5 that either match J 2 and J 4 or all three J 2 -J 6 . The posterior distribution of MoI value narrows substantially with every additional constraint.
Abstract models that match J 2 -J 6 yield a MoI range from ∼0.2180 until a sharp drop off at 0.2189. Our physical models yield a MoI value of 0.2181 with a 1-σ error bar of 0.0002. Broadly speaking the predictions from the two ensembles are compatible. However, with increasing spheroid number, our abstract models cluster around the most likely value of 0.2188, which is a bit higher than our physical models predict. This difference is a consequence of the   Fig. 3). The third panel shows the same trend in four-spheroid models of Saturn. The lowest panel shows predictions from models that match Saturn's J2 through J6. As the number of spheroids is increased, models tend to cluster in a narrower MoI interval. With 50 spheroid models that match J2 through J6, we obtained a range from 0.26393-0.26398 for Jupiter's MoI. (Between 5 × 10 5 and 6 × 10 7 models were constructed to compute every individual MoI histogram.)   16 that becomes exact in the limit of small qrot and large J2. In this limit, the 50-spheroid models are forced to approach the uniform-density limit of MoI= 2 5 . One can also approach this limit with two-spheroid calculations. So in panel d), we show how the volume fraction of the inner spheroid changes as a function of qrot. As the fraction approaches unity, the choices for qrot and J2 imply a constant planet.
way the two ensembles are constructed. In one case, we apply a number of physical assumptions. In the other, we do not and let the Monte Carlo procedure gravitate towards the most likely parameter space as long as the spacecraft measurements are reproduced. So one may expect to see small deviations in the predictions of the two ensembles.

Giant planets in general
The results in Fig. 5 show that the MoI of a giant planet can already be constrained reasonably well even if only q rot and J 2 are known. We therefore derive the MoI for a set of hypothetical giant planets by performing MC calculations with N S = 50 spheroids on a grid of q rot and J 2 points, which will help us to understand why Jupiter's and Saturn's MoI differ by ∼20%.
The ensemble averages of the computed MoI are shown in Fig. 6. One finds in Fig. 6a that for a given q rot , the MoI rises rapidly with increasing J 2 . To first approximation, J 2 is a measure of the planet's oblateness. So if J 2 is increased, while the equatorial radius and the rotation period are kept constant, more mass is moved towards the equatorial region, increasing the MoI. In Fig. 6b, we also show the predictions of the Radau-Darwin approximation, While there exist slightly different formulations of this approximation (Zharkov & Trubitsyn 1978), they all become exact in the limit of small q rot and large J 2 . In this limit, the planet's density becomes more and more uniform throughout its interior. Eventually the MoI approaches 2 5 , the value for a uniform-density fluid planet (Maclaurin spheroid) regardless of rotation rate. The 2 5 value cannot be exceeded unless one permits the density in the interior to be less than that of the exterior, which we exclude from consideration.
The uniform-density limit is also approached by models that have just two spheroids. While we fix the parameters of the outer spheroid, ρ 0 = 0 and λ 0 = 1, the two parameters of the inner spheroids, ρ 1 ≥ 0 and λ 1 ≤ 1, are just sufficient to match a pair of prescribed q rot and J 2 values. In Fig. 6d, we plot the volume fraction of the inner spheroid as function of q rot . When this fraction approaches 1 for small q rot , the density of the planet becomes uniform. For a given J 2 , this occurs at the same q rot value that leads to a MoI value of 2 5 in Fig. 6b. The two-spheroid calculations in Fig. 6d also confirm the trends that we see in the N S spheroid calculations in the other figure panels: With increasing q rot , more and more mass needs to be concentrated in the planet's center to satisfy the J 2 constraint. This leads to a decrease in the MoI if q rot is increased for a given J 2 , explaining the trends in Fig. 6b.
Finally we performed calculations for our two-spheroid models for Saturn's and Jupiter's q rot and J 2 values. While such models are crude, they show that the volume fraction of the inner spheroid is ∼54% for Jupiter and only ∼42% for Saturn. This implies that a higher fraction of Saturn's mass is concentrated near the the center, consistent with the fact that typical Jupiter models have a dilute core, while Saturn models matching the gravity measurements typically do not require one.

Jupiter
In Fig. 3a, we compare the posterior distributions of abstract Jupiter models with different numbers of spheroids. All models were constructed to match Jupiter mass, equatorial radius, and J 2 exactly. Models with fewer spheroids tend to show a wider range of J 4 and J 6 values, which is counterintuitive because, e.g., the entire space of 10 spheroid models is included in that of the 20 spheroid models. (In an 20 spheroid model, one only needs to set ρ 2i = ρ 2i+1 to obtain a valid 10 spheroid model.) However, the available space of 20 spheroid models is much bigger and in most models, the magnitude of the density steps, ρ 2i ≤ ρ 2i+1 , is smaller than that between two densities in a 10 spheroid model. In most 20 spheroid models, the density varies slightly more gradually than in the coarser 10 spheroid models. As a result, a representative set of 20 spheroid models occupies a smaller area in J 4 -J 6 space than a set of 10 spheroid models. Despite this reduction with increasing N S , the range of every model ensemble includes the J 4 and J 6 values from the Juno measurements (Durante et al. 2020) as well as the predictions from the static gravity terms (no DR) according to the dilute core models from Militzer et al. (2022). We will refer to them as five layer models throughout this article.
In Fig. 3b, we compare the average MoI as function of J 4 and J 6 . In general, small J 6 and J 4 , that are less negative, lead to larger MoI values. One also notices that as J 6 is increased for a given J 4 , the MoI goes through a maximum and the Juno measurements place Jupiter in the regime where ∂MoI ∂J int 6 < 0 while the opposite is true for Saturn. From the shape of contour lines, we can infer that Jupiter's MoI is almost independent of J 4 .
The five layer models from Militzer et al. (2022) predict DR contributions to Jupiter's J 6 to be negative: -0.27 ×10 −6 or -0.8%. They are much smaller in magnitude than for Saturn (it was +6%) and have the opposite sign. However, since ∂MoI ∂J int 6 also has the opposite sign, we are again in a situation where models matching the gravity data with DR effects predict a smaller MoI than models without DR. The magnitude of the MoI difference between the two types of models is, at -0.01%, much smaller for Jupiter while it was -0.4% for Saturn.
While a -0.01% correction was derived from our more recent five layer models , one may also ask whether the DR effect could make a larger contribution to J 6 . Our preliminary Jupiter model (Hubbard & Militzer 2016b), put together before Juno data became available, differs in J 6 by -0.8 ×10 −6 from the now-available gravity data. Even if such a large discrepancy came from DR effects, the MoI would only decrease by -0.04%, still smaller than the 0.1% precision that Juno is expected to ultimately achieve for the MoI measurements.
In Fig. 7, we compare the MoI of two and three layer models for Jupiter's interior (Saumon & Guillot 2004b;Guillot et al. 2004;Militzer et al. 2008) that are based on a physical EOS for the hydrogen-helium mixture but do not contain sufficient flexibility to match all observations. The predicted MoI values range from 0.26385-0.26400. In panel 7b, the temperature of Jupiter's interior was increased by raising the 1 bar temperature step by step from the value of the Galileo entry probe, 166.1 K, up to the extreme value of 185 K (Miguel et al. 2022). Raising 1 bar temperature lowers the density of the hydrogen-helium mixture, which enables one to add more heavy elements and thereby produce models that have at least a protosolar heavy element abundance, Z PS = 1.53%. An increase of 10 K allows one to approximately add one Z PS worth of heavy elements to an existing model. Still most models require the transition Figure 7. Predictions from two and three layer models for Jupiter's interior that were constructed under different assumptions and match J2 exactly. Colors label models consistently across all panels but not all curves are shown in every panel for clarity. The numbers specify the assumed transition pressure in GPa between molecular and metallic layers so that models can be traced across different panels. Panel (b) compare the abundance of heavy elements in the outer molecular layer, Z1, with the value of the protosolar nebula, 1.53% (Lodders 2010). No winds were included except for the last models shown in light blue color. Only the brown curves refer to models with a compact core (6 Earth masses, rocky composition). The light green star shows the preferred model from Hubbard & Militzer (2016b) while the pentagons and triangles indicate other compact core models based on the EOSs by Militzer & Hubbard (2013) and Saumon et al. (1995a). pressure to be 400 GPa or higher, which is not compatible with predictions for the metallization of hydrogen and for the hydrogen-helium immiscibility. Both are assumed to occur at approximately 80-100 GPa (Morales et al. 2010).
Like the abstract models in Fig. 3, all physical models in Fig. 7 match J 2 exactly but the fact that the equation of hydrostatic equilibrium is satisfied and that a physical EOS is employed means that J 4 and J 6 are now much more tightly correlated. While abstract models permitted a wide interval of J 6 values from 32.5-36.5×10 −6 for J 4 = −587 × 10 −6 , the more physical assumptions narrow this range to 34.2-34.5×10 −6 in Fig. 7a.
In Fig. 8, we compared the MoI from ensembles of interior+wind models that match the entire set of Juno's even and odd gravity coefficients up to J 10 (Durante et al. 2020). The posterior distribution of our five-layer reference ensemble is centered around the MoI value of 0.26393, which we consider to be our most plausible prediction for Jupiter's MoI. If we increased the 1 bar temperature to 170 K, the resulting ensemble of MoI shifted to higher MoI values by a modest amount of ∼ 7 × 10 −6 . Slightly larger shifts were obtained when we changed the H-He EOS by reducing the density by 3% over a selected pressure interval ). The largest positive shift was obtained for a density reduction from 10-100 GPa and the largest negative was seen if the density was reduced from 50-100 GPa. Both MoI shifts were on the order to 10 −5 , which is why we report 0.26393±0.00001 for Jupiter's MoI.
In Fig. 9 we plot results from an ensemble of five layer models in order to show how the computed MoI correlates with different gravity harmonics. The MoI correlates positively with J 2 , negatively with J 4 , and not in a significant Probability density Reference model with five layers (N m = 5 × 10 5 ) 3% density reduction from 10-100 GPa (N m = 2 × 10 5 ) 3% density reduction from 50-500 GPa (N m = 6 × 10 4 ) Increased T 1 bar =170 K (N m = 1 × 10 6 ) Figure 8. Probability density distributions of normalized moments of inertia derived from different ensembles of interior models matching even and odd gravity harmonics up to J10. All models include a dilute core and contributions from winds. The red circles represent our reference models with five layers. When four layer models are constructed by removing either the helium rain layer or the core transition layer, the MoI distribution hardly changes. If the 1-bar temperature is increased from 166.1 to 170 K, the MoI increases by only ∼10 −5 . Slightly larger changes are seen when the density of the hydrogen-helium mixture is reduced over pressure intervals from 10-100 or from 50-500 GPa. In the caption, we specify the size of the ensemble that was used for each histogram.
way with J 6 . (The correlations differ from predictions of two and three layer models in Fig. 7 because they only match the Jupiter's mass and J 2 .) While the sign and slopes of the correlation of the MoI with J 2 and J 4 in Fig. 9 differ, one needs to consider that the sign and the magnitude of J 2 and J 4 differ as well (see Tab. 1). If one removes that dependence by evaluating J 2 ∂MoI ∂J2 = 10 −7 and J 4 ∂MoI ∂J4 = 8 × 10 −8 , one finds the correlations between the MoI and both gravity coefficients are rather similar. The small magnitudes of ∼10 −7 illustrates that an individual gravity coefficient would need to change a lot to alter the MoI significantly. Fig. 9 also shows that the posterior distributions of J 2 and J 4 are centered at the Juno gravity measurements as expected.
In Fig. 10, we demonstrate fairly good agreement between the density profiles of our abstract and physical models for Jupiter's interior. For a fractional radius of 0.2 and larger, the density of our physical five layer reference model falls within one standard deviation from the mean of the abstract ensemble that matches the planet's mass, equatorial radius and the gravity coefficients J 2 and J 4 . Both gravity coefficients do not constrain the core region very well and the abstract models can thus yield larger density values there. As expected, models that are only constrained by J 2 show a wider range of density values for given radius. Larger density values favored for r < 0.3 and smaller values for r > 0.4. Still for most radii, we find that the predictions from the J 2 and J 4 constrained models fall within one standard deviation of the J 2 constrained models.
In Tab. 2, we compare our result with different predictions for Jupiter's MoI in the literature. Early determinations based on Pioneer and Voyager measurements by Hubbard & Marley (1989) and Wisdom (1996), who assumed uniform rotation, predicted Jupiter's MoI to be 0.2640, which is very close to the 0.26393 ± 0.00001 value that we derived when we match the Juno measurements with models that included DR effects. This now preferred value is also included in the ranges from earlier CMS calculations by Hubbard & Militzer (2016b) and Wahl et al. (2017a). With a low-order ToF, Helled et al. (2011) predicted smaller MoI values. Nettelmann et al. (2012a) predicted a very wide range of MoI values because not all models were constructed to match J 2 , J 4 , and J 6 . Ni (2018) adopted the approach from Anderson & Schubert (2007) when he adjusted coefficients of a polynomial function for the density profile in Jupiter's interior in order to match the first gravity measurements of the Juno spacecraft. With the theory of figures, he obtained  Abstract models that match J 2 Abstract models that match J 2 and J 4 Physical interior model Figure 10. Comparison of the density profiles of abstract models with that of our physical reference model with five layers. The shaded regions represent the standard deviations in density among the abstract models for a given radius.

Planet/reference
Rotation period J2 × 10 6 J4 × 10 6 J6 × 10 6 C/M R 2 e Jupiter 9:55:29.711 h Our 5 layer model Militzer et al. (2022) 14696. a range for Jupiter's MoI that includes our most reliable value. In the lower part of the table, we show how the range of predicted MoI shrinks when more and more of Juno's gravity harmonics are reproduced. In Tab. 3, we compare our predictions for Jupiter and Saturn with results from CMS calculations that we performed for models that other authors had constructed with the theory of figures. The central quantity of this approach is the volumetric radius, s, of different interior layers. When we read in the model files from other authors, we construct a with 50 spheroids for Jupiter's interior that match the measured harmonics J 2 , J 4 and J 6 , we derived a narrower range of possible MoI values from 0.26393-0.26398. Finally with our most plausible five layer models for Jupiter's interior, we predict the planet's MoI to be 0.26393 ± 0.00001, which is about ∼10% above the critical value of C/M R 2 = 0.236 for the planet to be in spin-orbit resonance with Uranus today (Ward & Canup 2006). Wisdom et al. (2022) argue that available high-precision measurements of Saturn's zonal harmonics suffice to infer a tight MoI range that rules out a current Saturn precession resonance with Neptune. By the same token, our predicted range for Jupiter's MoI needs to lie within the range constrained by Juno's extended mission measurement of MoI.