© The Authors 2024

1 Introduction

Unveiling Jupiter’s internal structure is key for studying and constraining its formation and evolution with implications for other giant planets (Vazan et al. 2018; Helled et al. 2022; Miguel & Vazan 2023; Helled & Stevenson 2024). NASA’s Juno mission (Bolton et al. 2017) provided precise measurements of Jupiter’s gravity field, which are essential for constraining its interior (Iess et al. 2018; Durante et al. 2020). Additional constraints come from atmospheric measurements taken by the Galileo entry probe and Juno (von Zahn et al. 1998; Seiff et al. 1998; Wong et al. 2004; Li et al. 2020), and the cloud-level temperature derived from Voyager radio occultations (Gupta et al. 2022).

Early Juno measurements suggested the need for a dilute core in Jupiter to match the hemispherically symmetric gravity field (Wahl et al. 2017), but this came at the expense of incompatibility with atmospheric composition measurements. To address this, an inward decrease in the heavy element mass fraction, Z (or metallicity), was proposed (Debras & Chabrier 2019), which was shown to be unlikely from an evolution modeling perspective (Howard et al. 2023b). Other approaches to fit the atmospheric metallicity involved either lowering internal densities by increasing the temperature at 1 bar (Miguel et al. 2022; Howard et al. 2023a), or arbitrarily modifying the equation of state (EOS) for hydrogen and helium (Nettelmann et al. 2021; Howard et al. 2023a). More recent studies considered the existence of a radiative layer that could separate Jupiter’s upper atmosphere from the deep Z-poor interior, allowing for atmospheric enrichment in metallicity (Howard et al. 2023b; Müller & Helled 2024).

Juno also revealed hemispherical asymmetries in Jupiter’s gravity field, which are attributed to deep winds that significantly affect both the symmetric and asymmetric components of the gravitational signature, adding constraints to the range of possible internal structures (Kaspi et al. 2018, 2023; Guillot et al. 2018). These internal flows were shown to reach a depth of -3000 km within Jupiter (Kaspi et al. 2018, 2023). Most interior models account for the dynamical contribution to the gravity field (ΔJ_2n) by allowing it to have some typical ranges (Debras & Chabrier 2019; Miguel et al. 2022; Howard et al. 2023a), which do not necessarily reflect the plausible solutions for each specific model and its solid-body properties. Militzer et al. (2022) takes an alternative approach, optimizing the cloud-level wind and wind decay profile for each interior model and its background density, allowing for very large ΔJ_2n by arbitrarily varying wind decay depth with latitude.

Modeling the interior of a gas planet constrained by its gravity field is a multidimensional, nonlinear, non-unique problem that is traditionally analyzed statistically through per parameter distributions and pair-wise correlations (e.g., Miguel et al. 2022; Howard et al. 2023a), or by selecting the best-fitting (preferred) models (e.g., Debras & Chabrier 2019; Militzer et al. 2022). Identifying typical, characteristic interior structures (i.e., combinations of interior parameters) can provide further information on the wide multidimensional solution phase space and highlight effective parameters that can simplify the representation of plausible structures, making the search for them more effective, reducing the high computational costs involved (Ziv et al. 2024).

In this study, we propose a fully consistent coupled interior and wind model where each interior model, computed using the accurate concentric Maclaurin spheroid (CMS) method (Hubbard 2013), is uniquely matched to Jupiter’s observed gravity field by using its own solid-body solution and a physically viable wind model. This method enables the generation of a large sample of plausible dilute core models for Jupiter, defining the permissible range for the dynamical contribution to the gravity field. Then, performing clustering analysis across the high-dimensional space of interior models simplifies them into a 2D phase space. This framework can also be adapted to explore the interiors of other giant planets.

In Sect. 2, we describe our modeling setup and the process used to generate a sample of plausible interior structures. Section 3 focuses on analyzing the range of plausible interior models and their associated observables. Section 4 presents a clustering analysis to identify key interior structures of Jupiter and shows how these structures simplify the internal model’s dimensionality. We conclude in Sect. 5.

2 Methods

In this section, we build on the methodology developed by Ziv et al. (2024), but also incorporated a wind model to obtain a sample of feasible solutions. Figure 1 illustrates the interior structure model and the schematic workflow used in this study. Details of the computational steps are outlined below and in Sect. 4.1.

2.1 Jupiter interior structure model

We calculated interior models using the CMS method (Hubbard 2012, 2013), based on a publicly available GitHub code¹ (Movshovitz et al. 2020). The CMS method assembles a rotating fluid planet from N constant density spheroids, constructed by their equatorial radii to compute the planetary mass, shape, moment of inertia, and gravitational moments by solving each spheroid’s gravitational and rotational hydrostatic equilibrium. Moreover, radial profiles of various physical properties are computed (e.g., density, pressure, entropy, and composition). The results of CMS are sensitive to the spacing grid of spheroids (Debras & Chabrier 2018), so we used the same equatorial grid as in Howard et al. (2023a) with N = 1041 spheroids, in similarity to the exponential grid shown to best reproduce Jupiter analytical polytropic solution (Debras & Chabrier 2018). The equatorial radius measured at 1 bar pressure, R_eq = 71 492 km (Lindal 1992) is unaltered in our models, so as Jupiter’s rotation period of Ω = 9.92492 h (Riddle & Warwick 1976).

The gravitational moments, representing the planet’s shape and mass distribution, are precisely measured by Juno (Durante et al. 2020): $$J_n=-\frac{1}{a^{n}M}{\displaystyle \int r^n}{P}_n(\sin \theta )\rho (r,\theta )d^{3}r,$$(1)

where n is the harmonic degree, a = R_eq is Jupiter’s equatorial radius, M is the planetary mass, r is the radial coordinate, θ is the latitude, P_n is the nth Legendre polynomial, and ρ is the density which can be expressed as ρ = ρ_s + ρ′, where ρ_s is the density from solid-body rotation (calculated by the CMS method), and ρ′ is the dynamical density due to the wind (Kaspi et al. 2010). The normalized moments of inertia ($\text{NMoI}\equiv \frac{C^2}{M{R}_{\text{eq}}^2}$, where C is the moment of inertia about the rotation axis) is computed according to Eq. (5) from Hubbard & Militzer (2016).

Following the setup of Miguel et al. (2022) and Howard et al. (2023a), we modeled Jupiter’s internal structure with four layers: an outer envelope, an inner envelope, a dilute core, and a possible compact core (see Fig. 1). The outer envelope is assumed to be isentropic, with its adiabat determined by the temperature at 1 bar, T_1bar. This temperature, acting as a proxy for the deep interior entropy, was measured by the Galileo entry probe as 166.1 ± 0.8 K and reanalyzed from Voyager radio occultations to reach 170.3 ± 3.8 K. We set the hydrogen and helium abundances to be consistent with the Galileo probe atmospheric measurements Y₁/(X₁+ Y₁) = 0.238, where X₁ and Y₁ are the mass fractions of hydrogen and helium in the outer envelope, respectively (von Zahn et al. 1998). The envelope’s metallicity, Z₁, was measured in the atmosphere to exceed the solar abundance (see Fig. 1 in Howard et al. 2023a, and references therein). Recent interior models continue to struggle with meeting this observational constraint (Howard et al. 2023b).

The lower-than-protosolar helium abundance measured in Jupiter’s atmosphere by the Galileo probe (von Zahn et al. 1998), assuming the planet’s overall helium abundance is protosolar, suggests a deep enrichment of helium within the planet. A region where hydrogen and helium are immiscible could separate the He-poor outer envelope from the He-rich inner envelope, with inward helium enrichment potentially caused by “helium rain” (Stevenson & Salpeter 1977; Mankovich & Fortney 2020; Howard et al. 2024). While the phase transition pressure defining this region, P₁₂, remains uncertain, numerical simulations suggest it ranges between 0.8 and 3 Mbar (Morales et al. 2013; Schöttler & Redmer 2018). We set both envelopes to have the same metallicity, Z₁ = Z₂. The helium mass fraction in the inner envelope and dilute core is adjusted to ensure the total planetary mean helium abundance aligns with the protosolar value of Y_proto = 0.278 ± 0.006 (Serenelli & Basu 2010).

The dilute core is modeled as a region with a gradual inward increase in heavy elements, following Miguel et al. (2022): $$Z=Z_2+\frac{Z_{\text{dillute }}-Z_2}{2}\left[1-\mathrm{erf}\left(\frac{m-m_{\text{dilute }}}{\delta m_{\text{dil }}}\right)\right],$$(2)

where Z_di1ute is the maximum mass fraction of heavy materials in the dilute core, m_dilute denotes the extent of the dilute core in normalized mass, and δm_dil controls the steepness of the Z gradient. We set m_dil = 0.075. Formation-evolution models of Jupiter suggest that the dilute core should be relatively small, extending up to 20% of Jupiter’s mass (Müller et al. 2020). Most current models do not feature a sufficiently small dilute core, with only models from Howard et al. (2023a) matching this theoretical constraint. We explored interior models with dilute cores extending between 11% – 60% of Jupiter’s mass, with maximum metallicity ranging from 0.06 to 0.45. Finally, a compact core composed entirely of heavy elements is considered, defined by its normalized equatorial radius, r_core, with its mass, M_core, determined through the CMS calculation. Recent studies resulted in compact cores sized 0–8 M_⊕ (e.g., Nettelmann et al. 2021; Miguel et al. 2022; Howard et al. 2023a). We explored interior models with a compact core ranging between 0 and 12% Jupiter’s radius, equivalent to masses of 0–9.6 M_⊕.

The variety of EOSs for hydrogen and helium being used and uncertainties related to their interpolation, provide further uncertainty in Jupiter’s plausible structures. This has been extensively studied by examining how using various EOSs affects the resulting interior models (Miguel et al. 2022; Howard et al. 2023a). In this study, we adopted the state-of-the-art pure H and He from Chabrier et al. (2019), and accounted for their non-ideal interaction using the tables from Howard & Guillot (2023). For a given pressure and temperature, we calculate densities using the additive volume law (AVL) including the mixing effects (Howard & Guillot 2023; Howard et al. 2023b): $$\frac{1}{\rho (P,T)}=\frac{X}{{\rho }_{\text{H}}(P,T)}+\frac{Y}{{\rho }_{\text{He}}(P,T)}+XY{V}_{\mathrm{mix}}(P,T)+\frac{Z}{{\rho }_Z(P,T)},$$(3)

where ρ_H, ρ_He, ρ_Z and X, Y, Z are the pure species densities and abundances of hydrogen, helium, and heavy elements, respectively. The volume change resulting from interactions between H and He is accounted for by V_mix. The entropy is similarly computed using the AVL: $$S(P,T)=X{S}_{\text{H}}(P,T)+Y{S}_{\text{He}}(P,T)+XY{S}_{\mathrm{mix}}(P,T)+Z{S}_Z(P,T),$$(4)

where S_H, S_He, S_Z are the pure species entropies. The entropy change due to mixing of H and He is represented by S_mix. We used the Sesame water EOS (Lyon & Johnson 1992) to account for the heavy elements. To compute the densities ρ (in g cm⁻³) as a function of pressure P (in Mbar) in the compact core we used the “rock” analytical formula from Hubbard & Marley (1989): $$P={\rho }^{4.406}\exp \left(-6.579-0.176\rho +0.00202{\rho }^2\right).$$(5)

Fig. 1 Schematic view of Jupiter’s interior model used in this study (left) and the exploration workflow (right). The model’s free parameters are shown. For each step of the exploration, we state which method is being used, which observables are used to constrain the solutions, and the number of plausible models resulting from this step. The prediction errors of NeuralCMS are denoted by ϵ (Ziv et al. 2024).

2.2 Exploration of plausible interior models with NeuralCMS

We reduced Jupiter’s interior structure to seven controlling variable parameters discussed in Sect. 2.1: T_1bar, Z₁, P₁₂, Y_proto, m_dilute, Z_dilute, and r_core. To explore the phase space of these parameters and find parameter combinations consistent with Juno gravity data, we used NeuralCMS², a deep learning approach regressing the CMS model (Ziv et al. 2024). Neu-ralCMS is a deep neural network model, predicting the even gravitational moments J₂ to J₈ and the mass given a combination of the seven interior parameters mentioned above. It is trained on a large sample of wide-ranged CMS-computed interior models of Jupiter allowing the performance of an extensive grid search for plausible interior models with a reduced computational runtime by a factor of 10⁵. NeuralCMS was applied to eliminate models inconsistent with Juno and reduce the parameter range. Due to the small measurement uncertainty of J2 and J₄, actual CMS calculations were performed to exclude models falsely predicted to be plausible.

The grid search procedure follows the approach described in Ziv et al. (2024), where all parameter combinations are explored using a uniformly spaced grid for each parameter with m grid points, resulting in m⁷ interior models per grid search iteration, without prior assumptions on the model parameters. Models consistent with Juno’s observations are selected based on the interior criteria in Table 1, added with a prediction error from NeuralCMS, ϵ, to the allowable deviation from Juno’s measurements. The interior criteria allow interior models to deviate significantly from the Juno 3σ uncertainty on J₂ and J₄, which is ~2 × 10⁻⁹ to investigate the permissible range of interior models allowed by a self-consistent wind model.

In this study, we build upon the two grid search iterations conducted by Ziv et al. (2024). The initial iteration considers the maximum prediction errors on the dataset used to train Neural-CMS, to roughly reduce the range of the interior parameters. The second grid search uses a denser grid, the 3σ prediction errors on the training dataset, ϵ_3σ, and the reduced parameter range, followed by CMS calculations to retrieve a sample of plausible interior models. To obtain a larger sample, a third grid search iteration was performed with the further narrowed parameter range from the second iteration, ϵ_3σ, and an even denser grid. The third and final iteration with NeuralCMS yielded approximately 80 000 plausible interior models within ϵ_3σ, which were recalculated with CMS to retrieve a sample of ~20 000 parameter combinations that satisfy the interior criteria (Table 1). The three grid search stages are summarized in Table A.1. We confirmed that NeuralCMS does not add uncertainty to the distribution of interior parameters by verifying that no plausible models were excluded by it. This was done by recalculating (with the CMS method) the third grid search iteration with a larger prediction error margin, close the maximum errors on the validation dataset for J₂, J₄, and M.

2.3 Coupled thermal wind model with modified surface wind

The observed cloud-level wind and gravity data impose constraints on the flow dynamics and, consequently, the dynamical contribution to Jupiter’s gravity field. In fast-rotating giant planets like Jupiter, large-scale flows are predominantly governed by geostrophic balance (Pedlosky 1987; Kaspi et al. 2009). Given the largely zonally symmetric nature of Jupiter’s flow, the leading order dynamical balance is between the flow gradient in the direction parallel to the axis of rotation and the latitudinal density perturbations (Galanti et al. 2017), given by a thermal wind (TW) balance: $$2\text{Ω}r\frac{\partial }{\partial z}\left({\rho }_{\text{s}}u\right)=g_{\text{s}}\frac{\partial {\rho }^{\prime }}{\partial \theta },$$(6)

where u is the wind, ρ_s (r) and g_s (r) are the static (solid-body) density and gravity fields, computed with CMS, z is the direction parallel to the axis of rotation, and ρ′ (r, θ) is the anomalous density field. For a detailed derivation of this form of the TW balance, see Kaspi et al. (2018). The wind-induced gravitational moments are then calculated based on the anomalous density field $$\text{Δ}{J}_n=-\frac{2\pi }{M{R}_{\text{eq}}^n}{\displaystyle \underset{-\pi /2}{\overset{\pi /2}{\int }}\cos }\theta d\theta {\displaystyle \underset{0}{\overset{R_{\text{eq}}}{\int }}{\displaystyle r^{n+2}}{P}_n(\sin \theta ){\rho }^{\prime }(r,\theta )dr,}$$(7)

where n = 2, 3, …,N. This balance was used to determine the depth of Jupiter’s winds using Juno’s gravity data (Kaspi et al. 2018), and similarly for Saturn using Cassini’s gravity data (Galanti et al. 2019).

In this study, we followed the methodology of Galanti et al. (2019), where the cloud-level wind profile and its decay are adjusted to best fit the Juno-measured even and odd gravity harmonics, J₂ to J₁₀. We accounted for uncertainties in the observed cloud-level wind (Tollefson et al. 2017), by selecting plausible models that are within 20 ms⁻¹ of the observed wind (see Fig. 2a). Firstly, we decompose the observed wind into the first N = 99 Legendre polynomials $u^{\text{obs}}(\theta )={\displaystyle {\sum }_{i=0}^N{A}_i^{\text{obs}}{P}_i(\sin \theta )}$, where the coefficients $A_i^{obs}$ set the latitudinal wind profile. During optimization, these coefficients are allowed to vary, resulting in a modified cloud-level wind, $u^{\text{sol}}(\theta )={\displaystyle {\sum }_{i=0}^N{A}_i^{\text{sol}}{P}_i(\sin \theta )}$. Next, the modified wind u^so1 (θ) is projected parallel to the axis of rotation and then decayed radially using a continuous, monotonic decay function that is a linear combination of exponential and normalized hyperbolic tangent functions (see Fig. 2b). The wind depth remains constant across all latitudes. Projecting parallel to the rotation axis and applying radial decay have been demonstrated to provide the most probable flow structure for explaining the gravity measurements (Kaspi et al. 2023). We note that Militzer et al. (2022) also uses this procedure, in addition to another degree of freedom that allows for latitude-dependent wind decay.

The CMS calculation provides radial profiles of each interior model’s equipotential surfaces, density, temperature, and composition. These profiles serve as the background (static) state for each interior model, which is then used in Eq. (6), to couple the interior model with a wind model. This approach allows for the consistent selection of interior models that match with the observed cloud-level wind within its uncertainties and the Juno-measured gravitational moments within Juno’s 3σ uncertainty (Durante et al. 2020). We applied the wind model to the ~20 000 models accepted by the interior criteria discussed in Sect. 2.2 (shown as blue histograms in Fig. 2) and identified 491 plausible interior structures that satisfy both the interior and the wind criteria (see Table 1). These wind-constrained models are shown in red histograms in Fig. 3.

Fig. 2 Optimized wind solutions for all models accepted by the wind criteria (Table 1). Panel a: latitudinal cloud-level wind profiles. Panel b: the radial decay function with depth. The blue profile represents the observed cloud-level wind (Tollefson et al. 2017), and the dashed yellow profile indicates the mean of all decay profiles.

3 Results I: Plausible range of interior structures

Figure 3 presents the distributions of interior model parameters and observables for the plausible interior structures that satisfy the interior criteria (blue histograms) and the subset that additionally meets the wind criteria (red histograms). In this section, we examine the differences between these two sets of models and analyze the resulting ranges of parameters and observables. Figure B.1 shows the pair-wise relations (corner plot) of all 491 plausible interior models.

3.1 Distribution of the observables

We find that the lower gravity harmonics J₂ and J₄ (Figs. 3a, b) are strongly constrained by applying the wind model, indicating that the explored range was sufficiently broad. The wind constraint is different from that used by Miguel et al. (2022) (black lines in Fig. 3), where they examined the plausible dynamical contribution to the even gravity harmonics through random sampling of ΔJ_2n, using a constant background state for the TW optimization described in Sect. 2.3. They selected plausible solutions based on similar wind criteria to those used in this study (Table 1). We note that the results from Miguel et al. (2022) presented here represent only the wind-allowed range for the static J_2n, not the distribution of plausible interior models. The comparison here is therefore intended for these allowable ranges, as also shown in Table 2. Comparing our results to theirs, we observe a different range for ΔJ₂ and a comparable distribution for ΔJ₄ (see Figs. 3a, b and Table 2). The higher-order even gravity harmonics exhibit a narrower range due to the applied interior criteria, meaning that their values are largely determined by the interior structure. For ΔJ₆, the observable which is more challenging to fit (Debras & Chabrier 2019; Militzer et al. 2022), we find a different range than Miguel et al. (2022) with the largest allowed deviation from Juno to be ΔJ₆ × 10⁶ = −0.22 (for an interior model with T_1bar = 171K and Z₁ = 0.009). For J₈ and J₁₀ we observe much tighter distributions compared to Miguel et al. (2022). Table 2 also compares our findings with two previous studies that determine the permissible range of dynamical contribution to the gravity harmonics matching the gravity field observations (Kaspi et al. 2020), and also magnetic field constraints (Galanti & Kaspi 2021). Our results are more consistent with these studies than with Miguel et al. (2022). For J₆, our deviations from the Juno measurements are larger than in the previous studies, consistent with the finding that the best-fit wind solutions favor smaller values for ΔJ₆ (Kaspi et al. 2023).

During our exploration, we allowed the planetary mass to vary within the uncertainty stemming from various values reported for the gravitational constant G (Tiesinga et al. 2021; Ziv et al. 2024). This is similar methodologically to allowing variations in the equatorial radius while keeping the mass constant, as done by Howard et al. (2023a) and Miguel et al. (2022). Figure 3g shows that the mass was not constrained by the wind model. Similarly to the gravity harmonics higher-order than J₄, the NMol is primarily determined by the interior structure, and the wind model further constrains it. We find values of NMoI = 0.26395 ± 0.00002 (3σ uncertainty) for the sample of 491 plausible models. This result is consistent with Militzer & Hubbard (2023), who derived values of 0.26393–0.26398 with abstract models that match Juno’s J₂ – J₆, and the range for the preferred model from Militzer et al. (2022), with 0.26393 ± 0.00001 (shown in green line in Fig. 3), which is on the edge of our allowed range. Neuenschwander et al. (2021) presented empirical structure models of Jupiter matching the measured J₂, J₄, and equatorial radius, but with a lower derived NMoI of 0.263408–0.263874.

The preferred model from Militzer et al. (2022), which matches Jupiter’s gravity field by allowing the wind depth to vary with latitude, is represented by the green lines in Fig. 3. Their study indicated that if the wind depth is kept constant across latitudes, the optimized cloud-level wind deviates by 50m s⁻¹ from the observed values, much more than the observed uncertainties (Tollefson et al. 2017), and well outside our wind criteria. Our modeling setup is different than their preferred model by the exact helium rain and the dilute core regions implementation, the EOS used for H-He mixture, where they use the tables from (Militzer & Hubbard 2013), but most importantly the wind approach. We show in Figs. 3a–f that their model is either on the edge of our allowed range for the gravitational moments and NMoI, or falls outside it. We find that while their model meets the interior criteria for the even gravity harmonics, it is excluded as a plausible structure after applying the coupled wind model. Fitting J6 is challenging with most available EOSs, as demonstrated by Howard et al. (2023a), who showed that achieving a good fit requires EOSs associated with high T_1bar values, above 180 K, much higher than the Galileo probe measurement of 166.1 K imposed by Militzer et al. (2022).

Fig. 3 Distribution of observables (a–g) and interior structure parameters (h–n) for plausible interior models. The gravity harmonics (a–e) are the static components, which are fitted to Juno’s measurements using the wind model. Panel o shows the compact core mass determined by rcore. Blue histograms correspond to models that satisfy the interior criteria (19 982 models), while red histograms represent those that also meet the wind criteria (491 models). The red vertical line marks Juno’s J2n measurements (Durante et al. 2020) and derived mass (Ziv et al. 2024). Black Gaussians represent the allowed range for the static gravity harmonics accounting for differential rotation $\left(J_{2n}^{\text{static }}=J_{2n}^{\text{Juno }}-\text{Δ}{J}_{2n}^{\text{dynamical }}\right)$ from Miguel et al. (2022). For J10 (e), the Juno measurement lies outside the range shown, but the black Gaussian covers the displayed values. The preferred model from Militzer et al. (2022) is shown in green lines (we show the inner edge of the He rain region compared to P12). The gravity harmonics (a-e) distance from the red lines represents ΔJ2n. The histogram color corresponds to Fig. 1.

3.2 Distribution of the interior physical parameters

Next, we present the range of parameters defining Jupiter’s interior, along with the resulting mass and distribution of heavy elements within the planet. First, it is important to note that NeuralCMS significantly narrows the ranges of parameters such as T_1bar, Z₁, m_dilute, Z_dilute, and r_core (Ziv et al. 2024). Figures 3h–o show that the distribution of the interior parameters remains largely unchanged when wind constraints are applied. We note that the peak values in the histograms do not necessarily correspond to a plausible combination of parameters. Our models indicate a preference for high T_1bar, and small compact core mass (M_core) ranging from 0 to 5.3 M_⊕, which is consistent with findings from previous studies (Nettelmann et al. 2021; Miguel et al. 2022; Howard et al. 2023a). Additionally, T_1bar and Z₁ exhibit strong positive correlation (Ziv et al. 2024), resulting in a higher number of models with increased envelope metallicity, reaching up to Z₁ = 0.0342 (2.3 times the solar value), associated with T_{1 bar} = 187 K.

The characteristics of the dilute core in our models are consistent with those calculated using the same EOS from Howard et al. (2023a), though with a different setup. The plausible range of m_dilute is found to be between 0.25 and 0.6 and between 0.072 and 0.204 for Z_dilute. We report values for the masses of heavy elements in the different regions in the interior (also shown in Fig. B.1) defined the same as in Howard et al. (2023a). The mass of heavy elements in both envelopes, corresponding to the envelope metallicity and not the enrichment in the dilute core M_Z,env ranges from 1.6 to 10.9 M_⊕ and is primarily determined by the value of Z₁. The excess mass of heavy elements in the dilute core region M_Z,dil ranges from 7.6 to 19.6 M_⊕. The total mass of heavy elements in the planet M_Z,total varies between 17.9 and 25.8 M_⊕. This means that our models have most of their heavy elements in the dilute core region instead of in the envelope, in agreement with Howard et al. (2023a).

Fig. 4 Normalized maximum variance for increasing number of clusters tested. Shown is the analysis for the full sample of plausible interior structures (thick blue), and the mean and standard deviation of the analysis for 10 randomly selected sub-samples with different sizes. The red circle marks our selection of four clusters for this analysis.

4 Results II: Characteristic interior structures

The range of plausible interior structures for Jupiter, shaped by the complex interplay of properties like composition, compositional gradients, and thermal structure, makes it challenging to identify typical structures. However, determining such representative structures can simplify our characterization of Jupiter’s interior, provide valuable multidimensional insights, and enable meaningful comparisons with other giant planets.

4.1 Clustering analysis based on the interior parameters

The variability among plausible interior structures can be explored and simplified through clustering analysis, which helps identify characteristic interior structures of Jupiter by classifying common combinations of parameters. In this study, we used the k-means method which was successfully used to characterize the behavior of Jupiter’s magnetosphere (Collier et al. 2020). We used the sample of 491 plausible interior models with their seven defining parameters (Fig. 1) for this analysis. To avoid biases, we standardized the parameters by subtracting the mean and dividing by the standard deviation.

The k-means algorithm starts by randomly selecting k cluster centroids and assigning each of the 491 data points to the nearest centroid based on Euclidean distance. The centroids are then updated as the mean of their assigned points. This process repeats until centroids no longer change. To ensure robust results, the algorithm runs multiple times with different initializations, selecting the setup with the lowest sum of squared distances. The analysis was performed using the Matlab “kmeans” function³. The optimal number of clusters was determined by an “elbow analysis”, which shows how variance decreases with an increasing number of clusters (Fig. 4). In our case, variance drops sharply up to four clusters, with only marginal improvement beyond that. Adding more than four clusters yields little additional simplicity despite a minor variance reduction at nine clusters.

To demonstrate that the sample size of 491 is sufficient for statistical analysis, we performed the same analysis on sub-samples of varying sizes, as shown in Fig. 4. Our results indicate that for sub-samples with fewer than N = 150 data points, the variability in the elbow analysis is too high to statistically validate the quality of the clustering. Additionally, for these smaller sub-samples, there is a significant reduction in variance when using more than four clusters, which is not observed in larger sub-samples. In contrast, for the larger sub-samples, the behavior of the elbow plot closely resembles that of the full sample, supporting the robustness of our clustering analysis.

Figure 5 shows the variability of the interior parameters across the four clusters. The parameters controlling the envelope (T_1bar and Z₁) are distinguished between clusters, where two clusters are characterized by high values of these parameters and the other two clusters have low values. Conversely, The transition pressure setting the helium rain region P₁₂, and the mean planetary helium mass fraction set by Y_proto are not classified within the four clusters. Similar to the envelope parameters, the parameters controlling the planet’s core (m_dilute, Z_dilute, and r_core) are also separated between clusters, but with two different clusters characterized by high (low) values for widiiute and r_core (Z_dilute), and vice versa for the other two clusters.

These relationships between the core-controlling parameters are consistent with previous studies (Miguel et al. 2022; Howard et al. 2023a). Additionally, we find that the smallest plausible dilute cores are associated with a relatively hot and heavy envelope. Overall, the clusters represent four unique combinations of two envelope configurations and two core configurations. Again, we note that the parameter combination represented by the mean values of each cluster is not necessarily a plausible combination.

Fig. 5 Interior parameters’ means (points) and standard deviations (error bars) within the four clusters, shown in different colors. Clusters 1 and 2 (red and blue) show high values of T1bar and Z1, while clusters 3 and 4 (green and yellow) show lower values of these parameters. Clusters 1 and 3 feature high values of mdilute and rcore, and low values of Zdilute, whereas clusters 2 and 4 display the opposite trend. We note that the compact core mass, Mcore, was not used in the clustering analysis.

Fig. 6 All 491 identified plausible interior models presented in the phase space defined by the product of the two envelope-controlling parameters and the ratio of the two dilute core-controlling parameters. Colors represent different clusters, consistent with Fig. 5. Black circles show the most observationally constrained interior models (see the observational criteria in Table 1). The legend displays the fraction of models assigned to each cluster. “Heavy” and “light” refer to metallicity.

4.2 Reducing the dimensionality to two effective parameters

As discussed above, the resulting clusters are distinguished by their envelope and core states, indicating that interior structures can be effectively represented in a 2D phase space with two key parameters representing these states. Figure 6 presents all 491 plausible interior structures, color-coded by their assigned cluster, according to two effective parameters: the product of temperature at 1 bar and envelope metallicity (which distinguishes between hot and heavy versus cold and light envelopes), and the ratio of dilute core extent to maximum dilute core metallicity (which differentiates between extended and light versus small and heavy dilute cores). The latter parameter also correlates with the compact core radius and mass (see Fig. 5), where an extended dilute core is associated with a larger (heavier) compact core. The descriptors used (i.e., cold versus hot, small versus extended, light versus heavy) are relative to the derived plausible range (Fig. 3), with “light” and “heavy” referring to metallicity (molecular weight). Figure 6 also presents the fraction of interior models within each cluster, showing a higher density of models with a hot and heavy envelope and a small and heavy dilute core. A schematic of the four characteristic Jupiter interior structures, including the full parameter range of each cluster, is provided in Fig. 7. Importantly, we show that Jupiter’s interior structure can be effectively explained using a 2D phase space, clearly representing the four distinguishable characteristic structures. Moreover, the overlap at the edges of clusters suggests they represent characteristic structures rather than distinct end members.

Additional observational and theoretical constraints can be applied to the sample of plausible interior structures to refine the selection. Models consistent with atmospheric measurements and theoretically bound for the helium rain region, meeting the observational criteria listed in Table 1, are marked with black circles in Fig. 6. These criteria include a temperature at 1 bar lower than 178 K (close to the upper bound from reanalysis of the Voyager radio occultations derived by Gupta et al. 2022), higher than solar atmospheric (envelope) metallicity, and a helium rain transition pressure lower than 3 Mbar (Morales et al. 2013). Most highlighted structures fall within a cluster characterized by a cold and light envelope and a small, heavy dilute core, suggesting that Jupiter’s interior structure might be confined to this specific characteristic structure or at least a narrow region in the phase space (Fig. 6). This structure is also highlighted with a thick frame in Fig. 7. Under these constraints, the minimal T_1bar required to accommodate at least a solar abundance of heavy elements in the envelope exceeds 175 K, which is higher than the upper limit of 174.1 K derived from Voyager radio occultations (Gupta et al. 2022), suggesting that the deep interior entropy might not be represented by the measured T _1bar.

Fig. 7 Schematic diagram of the four characteristic interior structures of Jupiter identified in this study. The interior parameters mean values and full range within each cluster are presented. The fraction of plausible interior models within each cluster is shown in the same color coding as in Figs. 5 and 6. The structure with the most observationally constrained interior models is shown with a thicker frame. The diagrams are not to scale.

5 Conclusion

In this study, we explored a wide range of Jupiter’s structural models with self-consistent interior and wind models constrained by Juno’s gravity field measurements, observed cloud-level wind, and atmospheric observations from both Juno and the Galileo entry probe. These models feature four layers: an outer envelope defined by the temperature at 1 bar and atmospheric metallicity, an inner envelope with the same metallicity but enriched in helium due to helium rain, a dilute core characterized by compositional gradients, and a possible compact core.

We present a robust framework for characterizing Jupiter’s interior using machine learning techniques, which can also be applied to study the interiors of other giant fluid (exo)planets. Our approach uses NeuralCMS, a deep neural network model, to efficiently reduce the range of plausible interior structures and explore the full parameter space. Subsequently, we compute interior models using the accurate CMS method, integrated with a wind model that has been demonstrated to best explain the observed asymmetric gravity field (Kaspi et al. 2023). We used a state-of-the-art EOS for hydrogen and helium; however, using alternative available EOS could yield different plausible ranges and distributions of interior models, potentially impacting the clustering analysis. This methodology ensures a consistent selection of plausible models.

The wind constraints result in a different range of plausible dynamical contributions to the gravity field (ΔJ_2n) compared to the values found by Miguel et al. (2022), who used random sampling of ΔJ_2n combinations to account for the wind in their interior models, as well as by Howard et al. (2023a). These differences may arise from variability between the background density profiles used to derive wind-induced gravity harmonics and differences in the EOS. We show that the even gravity harmonics J₆ to J₁₀ are predominantly determined by the interior structure, thus suggesting that our derived range more accurately reflects the wind effects relevant to interior modeling.

The application of wind constraints has further implications: while the distribution of parameters controlling the interior does not change significantly by applying the wind model, the range of plausible interior parameters is more narrowly defined. This is particularly important for addressing physical features that recent interior models have struggled to reconcile with observations and theory, such as the temperature at 1 bar, atmospheric metallicity, and the extent of the dilute core.

Clustering analysis reveals combinations of seven interior parameters with similar characteristics, offering deeper insights than traditional pair-wise relations or preferred models. We identify four characteristic interior structures of Jupiter that differ in their envelope and planetary core configurations. Specifically, the envelope is categorized as either hot and heavy or cold and light, and the core configuration varies between a small, heavy dilute core with a small compact core or the opposite (relative to the derived plausible interior range). These four different configurations are expected from a mass balance perspective.

The classification process enabled us to simplify Jupiter’s interior characteristics into two effective parameters: one representing the envelope state and the other representing the core state. This reduced the dimensionality of the analysis from a 7D space to a 2D phase space, clearly highlighting the four characteristic structures. Within the reduced 2D phase space, we highlight the models most consistent with the combination of atmospheric measurements (T_1bar and Z₁) and show that they mostly fall within a specific characteristic cluster, thus providing a further important reduction towards one key internal model. These models are characterized by a small and heavy dilute core with mean values of m_dilute = 0.36 and Z_dilute = 0.16. Our refined results on the extent of the dilute core are consistent with previous studies by Debras & Chabrier (2019) and Militzer et al. (2022), but diverge from the formation-evolution models of Müller et al. (2020), which suggest a less extended dilute core. However, even these constrained models face challenges in matching with both supersolar atmospheric metallicity and the observationally consistent 1 bar temperature, indicating that these measurements might not fully represent Jupiter’s entire envelope within this modeling setup.

Acknowledgements

We thank the referee for valuable comments and suggestions, which have improved the quality of this manuscript. We thank the Juno Interior Working Group for useful discussions and comments, and Moria Abu for her graphic assistance with the schematic figures. We acknowledge support from the Israeli Space Agency and the Helen Kimmel Center for Planetary Science at the Weizmann Institute.

Appendix A Grid search iterations with NeuralCMS

Appendix B Corner plot of plausible interior models

Figure B.1 shows the pair-wise relations of the sample of 491 plausible interior models.

Fig. B.1 Resulting pair-wise relations and histograms of the observables and the interior parameters for the plausible models. The gravitational moments are scaled by 106. The parameters for the distribution of the mass of heavy elements (MZ) are defined the same as in Howard et al. (2023a). Each point represents a specific interior model. Error bars are not shown for better readability.

References

All Tables

All Figures

Fig. 1 Schematic view of Jupiter’s interior model used in this study (left) and the exploration workflow (right). The model’s free parameters are shown. For each step of the exploration, we state which method is being used, which observables are used to constrain the solutions, and the number of plausible models resulting from this step. The prediction errors of NeuralCMS are denoted by ϵ (Ziv et al. 2024).

In the text

	Fig. 2 Optimized wind solutions for all models accepted by the wind criteria (Table 1). Panel a: latitudinal cloud-level wind profiles. Panel b: the radial decay function with depth. The blue profile represents the observed cloud-level wind (Tollefson et al. 2017), and the dashed yellow profile indicates the mean of all decay profiles.
In the text

Fig. 3 Distribution of observables (a–g) and interior structure parameters (h–n) for plausible interior models. The gravity harmonics (a–e) are the static components, which are fitted to Juno’s measurements using the wind model. Panel o shows the compact core mass determined by rcore. Blue histograms correspond to models that satisfy the interior criteria (19 982 models), while red histograms represent those that also meet the wind criteria (491 models). The red vertical line marks Juno’s J2n measurements (Durante et al. 2020) and derived mass (Ziv et al. 2024). Black Gaussians represent the allowed range for the static gravity harmonics accounting for differential rotation $\left(J_{2n}^{\text{static }}=J_{2n}^{\text{Juno }}-\text{Δ}{J}_{2n}^{\text{dynamical }}\right)$ from Miguel et al. (2022). For J10 (e), the Juno measurement lies outside the range shown, but the black Gaussian covers the displayed values. The preferred model from Militzer et al. (2022) is shown in green lines (we show the inner edge of the He rain region compared to P12). The gravity harmonics (a-e) distance from the red lines represents ΔJ2n. The histogram color corresponds to Fig. 1.

In the text

	Fig. 4 Normalized maximum variance for increasing number of clusters tested. Shown is the analysis for the full sample of plausible interior structures (thick blue), and the mean and standard deviation of the analysis for 10 randomly selected sub-samples with different sizes. The red circle marks our selection of four clusters for this analysis.
In the text

Fig. 5 Interior parameters’ means (points) and standard deviations (error bars) within the four clusters, shown in different colors. Clusters 1 and 2 (red and blue) show high values of T1bar and Z1, while clusters 3 and 4 (green and yellow) show lower values of these parameters. Clusters 1 and 3 feature high values of mdilute and rcore, and low values of Zdilute, whereas clusters 2 and 4 display the opposite trend. We note that the compact core mass, Mcore, was not used in the clustering analysis.

In the text

Fig. 6 All 491 identified plausible interior models presented in the phase space defined by the product of the two envelope-controlling parameters and the ratio of the two dilute core-controlling parameters. Colors represent different clusters, consistent with Fig. 5. Black circles show the most observationally constrained interior models (see the observational criteria in Table 1). The legend displays the fraction of models assigned to each cluster. “Heavy” and “light” refer to metallicity.

In the text

Fig. 7 Schematic diagram of the four characteristic interior structures of Jupiter identified in this study. The interior parameters mean values and full range within each cluster are presented. The fraction of plausible interior models within each cluster is shown in the same color coding as in Figs. 5 and 6. The structure with the most observationally constrained interior models is shown with a thicker frame. The diagrams are not to scale.

In the text

	Fig. B.1 Resulting pair-wise relations and histograms of the observables and the interior parameters for the plausible models. The gravitational moments are scaled by 106. The parameters for the distribution of the mass of heavy elements (MZ) are defined the same as in Howard et al. (2023a). Each point represents a specific interior model. Error bars are not shown for better readability.
In the text