© The Authors 2026

1 Introduction

Inferring a planet’s bulk composition and internal structure provides a powerful diagnostic of its formation pathway (Helled & Bodenheimer 2011; Mordasini et al. 2016; Turrini et al. 2021). Gas giant planets represent an ideal case study in this context: as the largest planets, they must have formed early, when ample gas and solids were still available for accretion. Their high masses and rapid growth mean that they preserve valuable information about the composition of the protoplanetary disk (Pollack et al. 1996; Helled & Morbidelli 2021; Mordasini et al. 2016). Hence, constraining their composition is key to probing the conditions under which planets form.

The observational landscape now includes over 6000 confirmed exoplanets (NASA Exoplanet Archive 2025), with a large, well-characterised subset of giants having measured masses and radii from space and ground-based surveys (e.g. Kepler, TESS, WASP, and HAT/HATSouth). The large dataset allows inferences of heavy-element content in individual planets and statistical analyses across populations, enabling the identification of more general trends constrained by bulk properties.

Interior models are designed to reproduce the observed bulk properties by assuming specific distributions of heavy elements. These models often explore two simple end-member configurations: a core+envelope (CE) structure, where a distinct heavy-element core is surrounded by a H/He envelope, and a fully mixed (FM) structure, where the heavy elements are homogeneously distributed throughout the planet (Thorngren et al. 2016; Howard et al. 2025b). However, gravity-field measurements and ring seismology for Jupiter and Saturn suggest more complex interiors. Interior models tailored to match these data suggest that the interiors of the two gas giants consist of compositional gradients of heavy elements or ‘dilute (fuzzy) cores’ (Wahl et al. 2017; Debras & Chabrier 2019; Mankovich & Fuller 2021; Miguel et al. 2022; Militzer et al. 2022; Howard et al. 2023; Ziv et al. 2024). Nevertheless, efforts in applying such a fuzzy or dilute core (DC) structure to exoplanetary studies are still embryonic. Recent retrieval studies have begun testing such structural prescriptions in individual exoplanets. For instance, Bloot et al. (2023) applied homogeneous and inhomogeneous models within a Bayesian framework to a limited sample, highlighting strong degeneracies that prevent quantitative discrimination between the two. Similarly, van Dijk & Miguel (2025) extended this approach by incorporating hot Jupiters, showing that current observations allow only qualitative assessments for a few well-characterised cases. In this study, we extend interior modelling beyond the traditional end-member cases to include composition gradients and DC structures that capture gradual heavy-element distributions.

In this work, we assess how the inferred bulk metallicity (Z) depends on the assumed interior structure for a sample of gas giant exoplanets. We further examine how different structural assumptions affect the resulting population-level trends. This study is important for the interpretation of JWST data (Carone et al. 2025; Deka et al. 2025) as well as future observations of giant planet atmospheres. Atmospheric measurements are currently being used to further constrain the planetary interiors (Thorngren & Fortney 2019; Wilkinson et al. 2024; Knierim & Helled 2025). Future missions such as PLATO and Ariel will provide accurate measurements of planetary radii, masses, ages, and atmospheric compositions, allowing tighter constraints on exoplanet interiors and formation histories. (Rauer et al. 2014; Tinetti et al. 2018).

2 Methods

2.1 Model parameters

We used the Code d’Evolution Planetaire Adaptif et Modulaire (CEPAM) (Guillot & Morel 1995; Howard et al. 2024) to simulate the planetary evolution. CEPAM numerically solves the 1D planetary structure and thermal-evolution equations (mass conservation, hydrostatic equilibrium, energy transport, and energy conservation) coupled to EoSs, which set density and entropy as functions of P, T, and composition. The initial catalogue of giant planets from Otegi et al. (2020); Parc et al. (2024) was first restricted to those with all four required observational parameters (mass, radius, age, and T_eq) available. This requirement alone substantially reduced the sample size. From this subset, we applied additional filtering, requiring relative uncertainties smaller than 5% in radius and 10% in mass, following Howard et al. (2025b), which yelded a final sample of 44 planets with masses between 0.12 and 5.98 M_J. Exoplanets with equilibrium temperatures T_eq > 1000 K were excluded from this study, since additional energy sources can significantly inflate their radii (Guillot et al. 2006; Thorngren & Fortney 2018).

The ice-to-rock ratio was kept constant for all the tested structures and was assumed to be a 1:1 ratio. For hydrogen and helium, the ratio between the two elements was assumed to be protosolar (Y/(X + Y) = 0.27) (Asplund et al. 2021).

The equations of state (EoSs) used are important and are known to affect the inferred composition (Baraffe et al. 2008; Miguel et al. 2016; Müller et al. 2020; Howard et al. 2025b). The EoSs used in this study are the SESAME water and dry sand for ices and rocks (Lyon & Johnson 1992), and the CMS19+HG23 EoS for H–He (Chabrier et al. 2019; Howard & Guillot 2023).

The atmosphere was treated as a non-grey, fully radiative atmosphere. Opacities were treated as two bands (visible and IR), and the mean opacities were tabulated accordingly (Freedman et al. 2008, 2014; Parmentier et al. 2015). Opacity enhancement caused by heavy elements follows the methodology described in Valencia et al. (2013).

For all internal structures, the P–T profile was assumed to be adiabatic, and we did not consider compositional mixing of heavy elements. The initial entropy was not prescribed directly but instead was determined by the assumed initial planetary radius. The initial radii are typically in the range of 1–1.5 R_J, depending on the planetary mass. The choice of the initial radius is not expected to significantly influence the evolution when considering sufficiently old planets and adiabatic cooling (Marley et al. 2007; Müller & Helled 2023). For completeness, we also considered a few cases where we assumed a non-adiabatic interior (see Section 4.1).

2.2 Tested interior structures

We considered different interior structures for the planetary dataset. The three main tested structures include the following:

• Core+envelope: this structure assumes a well-defined core consisting of heavy elements, with the envelope being composed of H-He. As the metallicity in the core Z_core is unity, we used the core mass M_core as a free parameter to infer bulk metallicities. For the basic CE scenario, there are no heavy elements in the envelope: Z_env = Z_atm = 0, where Z_env is the metallicity in the envelope and Z_atm in the atmosphere.

• Core+envelope with stellar-abundance atmosphere (CEA): to assess the contribution of metals mixed into the envelope, we allowed heavy elements in the upper planetary layers Z_env = Z_atm > 0.

• Fully mixed: this structure assumes a homogeneous composition throughout the planetary interior. No core is present. The bulk metallicity equals to Z_env = Z_atm.

• Dilute core: this structure represents a gradual distribution of heavy elements. Following Miguel (2022); Bloot et al. (2023); Howard et al. (2023), we describe the structure using the following function: $$Z(m)=Z_{\mathrm{e}\mathrm{n}\mathrm{v}}+\frac{Z_{\mathrm{d}\mathrm{i}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{e}}-Z_{\mathrm{e}\mathrm{n}\mathrm{v}}}{2}[1-\mathrm{erf}\left(\frac{m-m_{\mathrm{d}\mathrm{i}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{e}}}{\delta m}\right)],$$(1)

  where Z_dilute is the maximum heavy–element fraction in the DC, the parameter δm controls the steepness of the gradient, and m_dilute is the mass coordinate at which the gradient is steepest (see Fig. 1). The value δm was fixed at 0.125 (as variations in this parameter are not believed to substantially affect planetary Z estimates (Bloot et al. 2023)). In the basic DC structure, Z_env = Z_atm = 0.

• Dilute cores with stellar- and super-stellar abundance atmospheres (DCA and DCA3): this structure is the same as that for DC but allows for a non-zero atmospheric metallicity. The atmospheric abundances tested match the host star’s (Z_⋆) metal abundance (Z_atm = Z_⋆). Similar to Jupiter’s atmospheric enrichment (Helled & Stevenson 2024), a more enriched case with Z_atm = 3 × Z_⋆ has also been tested (DCA3). Note that in these cases, m_dilute was fixed to a value of 0.5.

For each structure that is consistent with the observational constraints for a given planet, we report M_Z (heavy-element mass) and Z (bulk metallicity), related by Z = M_Z/M, where M is the planetary mass.

2.3 Root-find Z retrieval

In this section, we describe the statistical tool designed to retrieve the possible ranges of Z that reproduce each planetary radius for every structural assumption. The retrieval of planetary bulk metallicities Z was designed as an inversion of the evolution model: instead of using M and Z to predict (R, Age), we seek the metallicity that reproduces the observed (M, R, Age) within their uncertainties. Our approach incorporates observational uncertainties by treating the mass and radius as normally distributed around their observed values (M ~ 𝒩(M_obs, σ_M) and (R ~ 𝒩(R_obs, σ_R)), while the age is assumed to follow a uniform distribution (Age ~ 𝒰(Age_min, Age_max)). The atmospheric metallicity was not used as an additional constraint, since this measurement is not available for the majority of the planets in our sample. The method proceeds in three steps:

• Step one – evolution runs: we first generated a set of cooling tracks with CEPAM across a broad grid of input parameters. Metallicities were sampled in the range 0 ≤ Z ≤ 1 with steps of 0.025, while masses were varied around the observed mean value, and sampled at the central value, ±1σ, ±2σ, and ±3σ. This ensemble of runs ensures that the parameter space captures both the observational uncertainties in mass and the possible degeneracies between metallicity and mass in shaping planetary radii. A dense coverage of Z values is particularly important for detecting degeneracies in the M–Z–R relation.

• Step two – grid interpolation: the model outputs from CEPAM (describable as a radius/age relation for different structures and chemical abundances) were interpolated to build a continuous parameter grid. While linear interpolation is generally sufficient, we employed cubic spline interpolation (via scipy.interpolate.interp1d (Virtanen et al. 2020)) since it better resolves curved behaviours and possible degeneracies in the cooling tracks. The interpolated grid was refined to a resolution of ΔZ = 0.0025, corresponding to 400 metallicity steps¹. This step provides smooth predictions of planetary radii across the 3D space (M, Z, Age), which are later required for the inversion procedure.

  Fig. 1 Three internal structures tested in this study, with their respective heavy element structure profiles: CE, all the heavy elements situated in a well-defined core; DC, a gradient applied to the planet’s core; and FM: a homogeneous envelope throughout the planet.

• Step three – root–finding and sampling: To invert the grid and estimate Z, we repeatedly sampled observational values of mass, radius, and age. For each planet, we drew n = 5000 tuples (M, R, Age). For each tuple, Brent’s method was applied to solve: $$f(Z)=R_{\text{pred}}(Z,M,\text{Age})-R_{\text{obs}}=0,$$(2)

  where R_pred is the interpolated model prediction and R_obs is the sampled observational value. This procedure yields the metallicity Z that reproduces the observed radius. Repeating this process n times provides a posterior distribution for a prior of Z ~ 𝒰(0, 1), capturing both measurement uncertainties and intrinsic model degeneracies. A similar approach was implemented in Müller & Helled (2023).

Root-finding retrieval provides the posterior distributions of metallicity Z for each assumed interior structure and for every exoplanet analysed. These distributions were smoothed using a kernel density estimate (KDE), and the corresponding mean and standard deviation of Z were extracted and used to investigate planetary trends.

3 Results

3.1 Structural properties

Figure 2 shows the retrieved curves for CE, FM, and three DC structures that assume a different extent of the DC (m_dilute = 0.25, 0.5 or 0.75). Since m_dilute remained fixed in the three different cases, Z_dilute was used as the adjusting parameter. This setup allowed us to explore the impact of changing the extent of the dilution curve. Although the analysis was performed for the full sample, we first focus on three planets with representative masses: K2-287b (0.315 M_J), TOI-1811b (0.972 M_J), and TOI-5153b (3.26 M_J). In all cases, CE yields a lower bulk metallicity than FM. We verified this offset by reproducing the approach of Howard et al. (2025b); the results are consistent with our pipeline (see Appendix B). The physical driver is that, in FM, heavy elements mixed into the envelope and atmosphere increase the mean molecular weight and, critically, the Rosseland mean opacity in the radiative zone. The higher opacity slows cooling and keeps radii larger at fixed age; matching the observed radius then requires higher heavy-element mass than in the CE case. Howard et al. (2025b) also found that FM evolution tracks are systematically warmer, especially at younger ages. The opacity effect is the main contributor in FM structures for higher expected bulk metallicities.

When examining the DC curves for different m_dilute values (Fig. 2), we find that DC models generally have lower values of Z compared to CE and FM models, and their solutions more closely resemble CE than FM. No systematic distinction or trend is evident across the three different m_dilute values tested, as the model-to-model differences (and their intrinsic uncertainties) remain smaller than the observational errors. Since these DC structures exclude atmospheric metallicity (Z_atm = 0), the observed effects are driven solely by the internal distribution of heavy elements. Introducing a gradient, rather than a sharp core boundary, redistributes heavy elements into outer layers, increasing the local density and pressure. This results in a lower radius over time, because a smaller heavy-element mass is required to reproduce the observed radius and age. Nonetheless, such an effect appears to be of lower magnitude compared to the difference in FM-CE estimates.

After testing different m_dilute structures, we next consider DC models with atmospheric metallicities set to one and three times the stellar value (DCA and DCA3, respectively). The results for the three representative planets are shown in Fig. 3. The complete dataset, including mean values and standard deviations for all structures, is presented in Appendix A. Since Z_atm (and thus, the opacity effect) is the varying parameter, DCA and DCA3 directly probe its effect on the Z estimates. As shown in Fig. 3, increasing the atmospheric metallicity systematically increases the inferred Z.

In the three cases shown, both DCA and DCA3 yield higher Z than CE. For some planets (such as HAT-P-12b, TOI-2010b, and TOI-2180b), however, DCA results converge with DC models, suggesting that the opacity effects are minor. Overall, Fig. 3 confirms that the atmospheric opacity significantly affects the retrieved heavy-element abundance, whereas density redistribution alone produces comparatively minor shifts. The CEA (bright-red) curves closely track the DCA estimates (especially for TOI-1811b and TOI-5153b), indicating that both models have similar Z. Because they use identical opacities and differ only in interior structure (distinct core in CEA vs. composition gradient in DCA), this agreement reinforces that opacity is the dominant lever over the density distribution.

However, while we find that FM models generally lead to higher inferred bulk metallicities, we note that DCA3 models can sometimes lead to even higher estimates, especially for low-Z planets (for example in the cases of TOI-4127b and TOI-2589b). Furthermore, the results are affected by the choice of Z_atm since our method does not allow for solutions where the inferred bulk metallicity is lower than the assumed atmospheric metallicity. This limitation particularly affects low-Z planets, as illustrated by the cutoff in the DCA and DCA3 posterior distribution of TOI-5153b shown in Fig. 3. Knowledge of the planetary atmospheric composition is therefore crucial for narrowing the parameter space of possible solutions.

Fig. 2 Retrieved curves for CE, FM, and DC0.25–0.5–0.75 structures (varying mdilute values) across three different planetary cases for different masses. The CE curves are shown in purple, FM in yellow, and DC in black, and can be distinguished through the different line style.

Fig. 3 Retrieved curves for CE, FM, CEA, DCA, and DCA3. The dashed lines (in burgundy and orange) show the atmospheric abundances assumed for DCA and DCA3 models, respectively. Such lines act as a lower threshold for these models as the assumed Zatm was kept fixed for a given planet.

3.2 Mass–metallicity relation

After analysing the overall behaviour of the metallicity estimates across different structural models, we next examined the full dataset of 44 planets for possible population-level trends in the inferred composition. In particular, we investigated the heavy-element mass (M_Z), and the bulk metallicity (Z), and their dependency on the planetary mass (M). To assess these relations, we performed a Bayesian regression. The parameters of the power law M_Z [M_⊕] = β M [M_J]^α were estimated using a Markov chain Monte Carlo (MCMC, Python package emcee (Foreman-Mackey et al. 2013), following Müller & Helled (2023) and Howard et al. (2025b)). Posterior distributions of the slope (α) and intercept (β) were sampled from 5000 chains, taking into account 1σ uncertainties in both axes. The outputs for each tested structure is reported in Table 1.

All models reveal a positive slope between M_Z and M, with similar slopes but different intercepts depending on the assumed structure (Fig. 4). The FM model systematically predicts the highest M_Z, while CE predicts the lowest (with β ≈ 36). This scaling is broadly consistent with the relation reported by Thorngren et al. (2016), Howard et al. (2025b), and Chachan et al. (2025). Our results suggest a steeper dependence, with DCA3 yielding α > 0.7 ± 0.07, only marginally inconsistent with Thorngren’s relation of $M_Z\propto \sqrt{M}$ when uncertainties are considered.

For this reason, we inspected the reliability of such curves by testing the correlation between α and β. The covariance analysis shows a weak but positive coupling (ρ ~ 0.2 for all the models), indicating that the slope and intercept are not fully independent, and that small shifts in one parameter affect the other. We then also evaluated the fits with a weighted least-squares (WLS) χ² test. In all cases, χ² ≫ 1, suggesting that the power-law relations do not formally explain the data within the quoted errors; therefore, the reliance on such curves should be treated with caution. Not meeting such reliability limits with the χ² statistic could be caused by (i) noise arising from the scatter in the data, (ii) a large diversity in planetary compositions, or both. In Appendix D, we further discuss the WSL goodness-of-fit approach.

Since the WSL test revealed the limited reliability of the retrieved power-law fits, it is essential to verify whether the apparent slopes of these relations genuinely reflect an underlying correlation within the data. We further tested the correlations of the dataset for each structural assumption using the non-parametric Kendall’s τ rank test. To account for observational uncertainties, instead of computing the correlation on the mean values of the data, a bootstrap resampling approach was employed. Here, 5000 synthetic datasets were generated by drawing from the reported error distributions, and τ was calculated for each realisation. The median τ values and associated p-values are listed in Table 2. For the M_Z–M relation, nearly all runs returned p-values well below the adopted significance threshold of p = 0.05, confirming the robustness of the positive correlation.

Next, we examined the relation between Z and M (Fig. 5). In contrast to the M_Z–M trend, the Z–M relation shows a clear negative correlation. The ranking of models remains similar: FM predicts the highest average Z, followed by DCA3, DCA, and CE. Compared to the M_Z–M relation, the data are more scattered, particularly at lower planetary masses. Kendall’s τ again indicates statistically significant negative correlations, although the correlation strength is weaker. A portion of the reduced significance may reflect mathematical coupling from dividing by the x-axis variable (Z = M_Z/M), which weakens the apparent correlation and inflates p-values. Median p-values remain generally below 0.05, but with a modest fraction of realisations exceeding the threshold (~2% for FM and DCA3, and ~10–11% for DCA and CE). Nevertheless, the correlations remain statistically compelling across all models.

A similar relationship was tested using the normalisation to the host star (Z/Z_⋆ vs M). The results show similar trends (both in best-fit and correlation tests) to those obtained for Z/M. Such an investigation would have been essential if we had artificially imposed an anti-correlation driven by the stellar metallicity (especially since one of the model assumptions is that the stellar composition reflects that the planet’s atmospheric structure). Although recent studies (Buchhave et al. 2014; Müller & Helled 2025) have found trends between planetary heavy-element mass and host star chemical composition, the dataset used in this study does not show any correlation between Z_⋆ and M.

Fig. 4 Relation between the retrieved heavy-element mass MZ and planetary mass M, shown for the four main structures tested (CE, FM, DCA, and DCA3). Shown as dashed lines in a grey gradient are curves indicating the trends for fixed planetary compositions (0%, 50%, 90%, and 99% gas composition, respectively). Bayesian fits are shown using a power-law function MZ = β Mα. A summary of the fitted parameters and correlation coefficients is presented in Tables 1 and 2.

Fig. 5 Relation between retrieved metallicity Z and planetary mass M, shown for the four main structures tested (CE, FM, DCA, and DCA3). Bayesian fits are shown using a power-law function, Z = β Mα. For additional quantitative results, we refer to Tables 1 and 2.

4 Discussion

4.1 Non-adiabatic structure

In our models with a DC (DC, DCA, and DCA3), we assumed an adiabatic temperature gradient. However, the presence of a composition gradient implies that the temperature gradient could differ from the adiabatic one. To assess the effect of a non-adiabatic temperature gradient on the inferred bulk metallicity, we parameterised the deviation from adiabaticity (likely caused by double-diffusive convection (Debras et al. 2021; Leconte & Chabrier 2012)) as in Howard et al. (2025a), using $${\mathrm{\nabla }}_T={\mathrm{\nabla }}_{\mathrm{a}\mathrm{d}}+R_{\rho }B,$$(3)

where $B=\frac{{\chi }_{\rho }}{{\chi }_{\mathrm{T}}}{\left(\frac{d\ln \rho }{d\ln Z}\right)}_{P,T}{\mathrm{\nabla }}_Z,{\chi }_{\rho }={\left(\frac{d\ln P}{d\ln \rho }\right)}_{T,Z},{\chi }_{\mathrm{T}}={\left(\frac{d\ln P}{d\ln T}\right)}_{\rho ,Z}$ and ${\mathrm{\nabla }}_Z=\frac{d\ln Z}{d\ln P}$.

For a smaller selection of planets, we simulated the planetary evolution using R_ρ = 0.1 and 0.4 and compared it to the adiabatic case (R_ρ = 0). The inferred bulk metallicities Z are shown in Fig. 6. We find that R_ρ = 0.1 and R_ρ = 0.4 lead to increased values of Z by up to 7% and 35%, respectively. The super-adiabatic temperature gradients obtained with R_ρ = 0.1 and R_ρ = 0.4 yield hotter interiors compared to the adiabatic case (see bottom panel of Fig. 6) and naturally lead to an increased amount of heavy elements required to match the planets’ radii. We note that efforts to constrain R_ρ for the Solar System’s gas giants have been recently presented (Sur et al. 2025), although the strong degeneracy in this parameter resulted in a broad range of admissible values. Previous studies estimated R_ρ for Jupiter and Saturn, suggesting R_ρ < 0.1 at around 0.05 (Mankovich et al. 2016; Mankovich & Fortney 2020). This led us to test values of 0.1 and an extreme case of 0.4. The actual value is model-dependent, and efforts to constrain it in giant planets still remains ongoing.

A clear limitation of our models is that convective mixing is not included. This process could affect the thermal evolution and the distribution of heavy elements. For example, Knierim & Helled (2025) showed that for a Jupiter-mass planet neglecting mixing could overestimate the radius by 0.08 R_J. The values of R_ρ we chose were arbitrary, and it is necessary to assess whether DCs with such a temperature structure can remain stable against convection. Future studies should further assess the effect of non-adiabaticity on the mass–metallicity relation and the thermal evolution of gas giants.

4.2 Connection with formation

Extent of the DC. For some high-metallicity cases (Z ≳ 0.5; e.g., TOI-4010d and WASP-196b), no solutions were found for compact gradients with m_dilute = 0.25 or 0.5: the interior could not accommodate the required heavy-element fraction unless the gradient was shifted outwards to higher m_dilute. Although this behaviour might suggest a more extended dilute region in these planets, this interpretation is not uniquely supported by the data. The inferred extent of dilution is strongly model-dependent and can arise simply from our assumption of a fixed Z_atm, which permits additional heavy elements to be stored in the envelope without altering the internal composition gradient. In more massive planets, a wide non-adiabatic DC could substantially increase the heat retained in the deep interior, altering the radius evolution. Constraining m_dilute is therefore crucial for massive planets. A well-confined DC may have little impact on inflation (the thermal profile remains close to adiabatic), whereas high m_dilute values can stabilise significantly more heat. In non-adiabatic cases, the steepness of the gradient may also become important.

Mass–metallicity relation. Our results, obtained across different structural configurations, are consistent with previous studies (Thorngren et al. 2016; Teske et al. 2019; Howard et al. 2025b; Chachan et al. 2025). We recover a clear positive correlation between the total heavy-element mass (M_Z) and planetary mass (M), and a weaker, yet significant, negative correlation between metallicity (Z) and M, also when normalised by Z_⋆.

The overall Z–M trend indicates a gradual decrease in planetary metallicity with increasing mass. There is a sharp decrease in Z for low planetary masses (0.2–0.3 M_J), followed by a flatter behaviour for higher masses. Other studies on mass–metallicity relations have found similar results (Chachan et al. 2025). At the high-mass end, we find hints of renewed enrichment. This enrichment, together with the significant dispersion across the mass–metallicity relation, likely reflects the diversity of formation pathways and interior structures of giant planets. The presence of high-Z, low-M objects (potentially with compositions more akin to Neptune-like planets) is consistent with planet formation theory (Helled 2023). Interestingly, we find that more massive planets (M > 4 M_J) are also metal-rich. Even if discerning any conclusion from two planets is premature, this may suggest a secondary accretion of heavy elements through late planetesimal capture or giant impacts after disk dispersal, leading to enhanced metal enrichment in their deep envelopes (Mordasini et al. 2016; Helled & Stevenson 2017; Venturini & Helled 2020; Shibata et al. 2022).

To assess the robustness of the inferred anti-correlation, we examined how the correlation significance changes when progressively removing planets from the mass extremes. We quantified this sensitivity by recalculating the Kendall τ correlation and the corresponding p-values after each removal for both the M_Z/M and Z/M relations (Fig. 7).

We find that removing low-mass planets has a small effect on the inferred M_Z/M relation (top row in Fig. 7): the p-values remain nearly unchanged for all the models. In contrast, removing high-mass objects systematically increases the p-values. These constitute the ‘anchor points’ that most strongly support the monotonic trend found with the correlation test. Even after the removal of six planets from either end, the correlation remains statistically significant, confirming that the observed trend is robust.

For the Z/M relation, the behaviour is different. Because Z already scales inversely with mass (Z = M_Z/M), plotting Z against M effectively introduces an additional mass dependence in the denominator, which leads to a higher inferred p-value (Kronmal 1993; Brett 2004). The removal of two low-mass planets is sufficient for the CE and DCA models to lose significance at the p = 0.05 level. After the removal of six planets, only the FM model retains a marginally significant p-value. In contrast, removing the most massive planets has a rather small impact. We find that removing high-mass planets from the sample increases the p-values for M_Z versus M, while removing low-mass planets increases the p-values for Z versus M. This is because Z = M_Z/M.

Our findings raise a classification question. Planets with Z ≳ 0.5 and low masses likely reside near the transition between Neptunian and Jovian regimes, where electron degeneracy pressure no longer dominates the M–R relation (Müller et al. 2024). The contrast in how the p-values increase when removing either low-M or high-M planets could hold a physical meaning. This may indicate the diversity of the sample and the potential existence of sub-populations that do not follow a monotonic trend. In this context, the assumption of a universal, continuous mass–metallicity trend is clearly oversimplified and cannot reflect the large diversity in planetary composition.

Fig. 6 Top: metallicity estimates for adiabatic (Rρ = 0) and super-adiabatic (Rρ = 0.1, 0.4) models. This is shown for four different planetary cases. Bottom: corresponding temperature profiles over the planet’s normalised mass, using the structures from the median metallicity of the DCA model at a fixed age of 4 Gyr.

Fig. 7 Kendall’s τ correlation p-value for MZ/M and Z/M with n-planets removed from the initial dataset. Planets are removed from the lowest-mass (left) and highest-mass (right) planets. The number of planets being iteratively removed can be seen on the x-axis. The four different models are respectively colour-coded. The horizontal line shows the p = 0.05 acceptance threshold.

5 Conclusions

We quantified bulk metallicities for transiting gas giants under multiple interior structures using CEPAM with a root–finding retrieval. Our key conclusions can be summarised as follows:

• We confirm that the inferred bulk metallicity for FM structures is higher than for CE structures, due to the opacity effect that prevails over the density redistribution effect. This, in turn, leads to DC structures that yield similar bulk metallicities to CE, independently of the extent of the composition gradient;

• When increasing the atmospheric metallicity in the DC structure to one (DCA) or three times stellar (DCA3), the inferred bulk metallicity becomes closer to the one yielded by FM. This demonstrates the importance of calibrating it with future observational data to better constrain internal metallicity estimates;

• When considering non-adiabatic temperature gradients in DC structures, the inferred bulk metallicity can increase by up to 35% for a used value of R_ρ = 0.4. However, future studies should explore a larger parameter space and assess the long-term stability of composition gradients;

• Using a sample of 44 giant exoplanets, we confirm that heavy-element mass increases with planetary mass. The slope of the M_Z–M relationship does not change with the assumed structure, but only the intercept;

• The planetary bulk metallicity Z generally decreases with planetary mass. Planets with M < 0.2 M_J have metallicities between 0.4 and 0.8. These low M – high Z planets are required to hold the negative correlation between Z and M. Planets with 0.2 < Mp < 4 M_J have a strong scatter in metallicities, which can range from less than 0.01–0.5, showing the diverse composition of giant exoplanets. For the two high planetary masses (M > 4M_J) in our sample, the inferred metallicities remain relatively high (between 0.1 and 0.3). This suggests that giant planets continue to accrete heavy elements after their formation (Shibata et al. 2022).

Future studies that account for additional physical processes such as convective mixing as well as the observed atmospheric compositions can improve giant planet characterisation. In addition, further investigations of breakpoints via multi-linear or orthogonal-distance regressions (Müller et al. 2024) in the inferred relations are desirable. Accurate measurements of the atmospheric chemistry of giant planets by JWST can also further constrain the planetary bulk composition and internal structure. In addition, forthcoming measurements from PLATO and Ariel (expected to be launched in 2026 and 2029, respectively) will sharpen bulk metallicity estimates by providing tighter constraints on mass, radius, age, and atmospheric composition for many planets. This will allow a better understanding of the interiors of giant planets and of their formation and evolution histories.

Acknowledgements

We thank Stefano Wirth, Simon Müller, Henrik Knierim, and Shay Zucker for insightful discussions.

References

Appendix A Planetary data

Appendix B Result comparison with Howard et al. (2025)

In order to assess the strength and accuracy of the retrieved Z for the sample used, we took advantage of the fact that the used dataset is the same as the one used by Howard et al. (2025b). We thus took the initial sample of 21 exoplanets with the most precise data on Radius (R_err<3%), and compared it with the same planets processed with the method explained in Howard et al. (2025b). Results are shown in Fig. B.1, where a visual comparison between the trends for the Z/Z_⋆ versus M is shown. To note how the best-fit lines vary from the ones retrieved with the full dataset of 44 exoplanets, showing the variability and unreliability of such curves for model prediction. Howard also discussed in their paper the shift of the slope α between the dataset, considering only R_err<3%, or for the broader range considered of R_err<5%. In the bottom panels of Fig. B.1, it is shown how we assessed the difference between the data from this paper compared to the ones from Howard et al. (2025b). While a simple subtraction could provide information about the difference between the two mean values, we also wanted to investigate whether each ranges of values (so accounting for data uncertainty) lie within the other. To do so, a normalisation using the standard deviations of the data has been computed. The normalised difference can be expressed with the following formula: $$\sigma =\frac{CE/F{M}_{Howard}-CE/FM}{\sqrt{{\sigma }_{\text{Howard\hspace{0.17em}}}^2+{\sigma }_{\mathrm{C}\mathrm{E}/\mathrm{F}\mathrm{M}}^2}}.$$(B.1)

Fig. B.1 Output comparison between the data retrieved in this study, and the one found in Howard et al. (2025b). Upper panels: Bulk metallicity (normalised by the star’s metallicity) over mass for the initial 21 planets (as mentioned in Fig. A.1 and A.2) both for FM and CE structures. Bottom panels: Normalised difference (calculated with eq. B.1) between the two estimates.

Results regarding this comparison were considered positive. While there is a slight change in the predictive trends between the two structures CE and FM (slight shift in β for CE, and in α for FM compared to the ones found in Howard’s. The results for the normalised difference ultimately confirmed the similarities between the two retrieved results, and their respective methods; All values lie between a standard deviation of 1 and −1 (or |σ|< 1), showing statistical consistency for all 21 cases compared. Furthermore, the normalised difference does not show any systematic error related to the mass (since there is no trend in the normalised difference σ). We thus conclude by saying that there is no discouraging difference between the two results (as the difference between the two approaches can be explained to be within the data error range and thus not due to randomness). We then accepted the validity of our approach, explained in Section 2.3, and proceeded with the inspection of the additional structures.

Appendix C Model assumption-driven changes

To isolate the impact of structural assumptions on the inferred metallicity, we quantified model-dependent offsets in Z across the four models (FM, CE, DCA, DCA3) by forming differences relative to FM, namely ΔZ_FM-CE ≡ Z_FM − Z_CE, ΔZ_FM-DCA ≡ Z_FM − Z_DCA, and ΔZ_FM-DCA3 ≡ Z_FM − Z_DCA3. Figure C.1 aims at showing such changes in the assumed structure over the full mass regime of the inspected exoplanets. The three lines show the best-fit lines of the data that is computed by the subtraction between the two mean points according to the models compared. No σ in the models has been accounted for, as only the mean points have been inspected. For this reason, no error protraction has been shown on the best-fit lines, which have been kept only for visual reasons. Nonetheless, without counting model errors, the ΔZ for all three models being compared to FM show statistical validity as shown by the Kendall’s τ results (and respective p-values) shown on the Figure’s legend. We can see that, on average, all the models tend to reach 0 (No Z difference between two structurally different models) in the high-M regime. While the best-fit line for ΔZ_FM–CE approaches zero, the difference between the FM-DCA and FM-DCA3 goes to the negative (meaning that FM is no longer the structural endmember, but DCA and DCA3 overpredict Z more than FM). This is caused by the forcing of respectively stellar and three times stellar abundances in their structures. Low-Z scenarios (usually in the high-M region) will be ‘ruled out’ by forcing an enriched metallicity in the atmosphere. This could suggest that low-Z planets are unlikely to be super-enriched in heavy elements in their atmospheres. This agrees with the theory exposed in Fortney et al. (2013). Further studies must be made to investigate if that is the case, or if degeneracies between mass distribution in massive planets might overshadow the overall magnitude of the opacity effect proved to be dominant by this study.

Fig. C.1 Change in bulk metallicity according to the assumed structure and their difference relative to the FM model over mass.

Appendix D WLS χ² test

To assess the reliability of the power-law fits, we performed a WLS χ² test, which accounts for individual data uncertainties rather than relying solely on R². In WLS, each residual is scaled by its measurement error: $${\chi }_{\mathrm{o}\mathrm{b}\mathrm{s}}^2=\sum _{i=1}^N{\left(\frac{y_i-f(x_i;\theta )}{{\sigma }_i}\right)}^2,$$(D.1)

where y_i are the data, f(x_i; θ) the model, and σ_i the 1σ uncertainties. Under the null hypothesis that the model describes the data and errors are Gaussian, ${\chi }_{\text{obs}}^2$ follows a chi-square distribution with ν = N − p degrees of freedom.

For N = 43 and p = 2, ν = 41, giving ${\chi }_{(0.05;41)}^2\approx 56.9$. Fits are accepted at the 5% level if ${\chi }_{\mathrm{o}\mathrm{b}\mathrm{s}}^2\le {\chi }_{(0.05;\nu )}^2$. In all cases, ${\chi }_{\mathrm{o}\mathrm{b}\mathrm{s}}^2$ exceeds this threshold, indicating that the fits do not reproduce the data within the quoted errors. For the Z–M relation, the observed values are: $${\chi }_{\mathrm{o}\mathrm{b}\mathrm{s}}^2=580.4(\mathrm{C}\mathrm{E}),465.4(\mathrm{D}\mathrm{C}\mathrm{A}),311.7(\mathrm{D}\mathrm{C}\mathrm{A}3),243.7(\mathrm{F}\mathrm{M}),$$

corresponding to reduced ${\chi }_{\text{red}}^2$ values of 14.2, 11.4, 7.6, and 5.9, respectively. Since ${\chi }_{\text{red}}^2$ ≫ 1, we reject the null hypothesis and conclude that the observed scatter exceeds that expected from measurement errors alone. The diversity of planetary metallicities likely reflects intrinsic astrophysical variability rather than observational noise.

All Tables

All Figures

	Fig. 1 Three internal structures tested in this study, with their respective heavy element structure profiles: CE, all the heavy elements situated in a well-defined core; DC, a gradient applied to the planet’s core; and FM: a homogeneous envelope throughout the planet.
In the text

	Fig. 2 Retrieved curves for CE, FM, and DC0.25–0.5–0.75 structures (varying mdilute values) across three different planetary cases for different masses. The CE curves are shown in purple, FM in yellow, and DC in black, and can be distinguished through the different line style.
In the text

	Fig. 3 Retrieved curves for CE, FM, CEA, DCA, and DCA3. The dashed lines (in burgundy and orange) show the atmospheric abundances assumed for DCA and DCA3 models, respectively. Such lines act as a lower threshold for these models as the assumed Zatm was kept fixed for a given planet.
In the text

Fig. 4 Relation between the retrieved heavy-element mass MZ and planetary mass M, shown for the four main structures tested (CE, FM, DCA, and DCA3). Shown as dashed lines in a grey gradient are curves indicating the trends for fixed planetary compositions (0%, 50%, 90%, and 99% gas composition, respectively). Bayesian fits are shown using a power-law function MZ = β Mα. A summary of the fitted parameters and correlation coefficients is presented in Tables 1 and 2.

In the text

	Fig. 5 Relation between retrieved metallicity Z and planetary mass M, shown for the four main structures tested (CE, FM, DCA, and DCA3). Bayesian fits are shown using a power-law function, Z = β Mα. For additional quantitative results, we refer to Tables 1 and 2.
In the text

	Fig. 6 Top: metallicity estimates for adiabatic (Rρ = 0) and super-adiabatic (Rρ = 0.1, 0.4) models. This is shown for four different planetary cases. Bottom: corresponding temperature profiles over the planet’s normalised mass, using the structures from the median metallicity of the DCA model at a fixed age of 4 Gyr.
In the text

	Fig. 7 Kendall’s τ correlation p-value for MZ/M and Z/M with n-planets removed from the initial dataset. Planets are removed from the lowest-mass (left) and highest-mass (right) planets. The number of planets being iteratively removed can be seen on the x-axis. The four different models are respectively colour-coded. The horizontal line shows the p = 0.05 acceptance threshold.
In the text

	Fig. B.1 Output comparison between the data retrieved in this study, and the one found in Howard et al. (2025b). Upper panels: Bulk metallicity (normalised by the star’s metallicity) over mass for the initial 21 planets (as mentioned in Fig. A.1 and A.2) both for FM and CE structures. Bottom panels: Normalised difference (calculated with eq. B.1) between the two estimates.
In the text

	Fig. C.1 Change in bulk metallicity according to the assumed structure and their difference relative to the FM model over mass.
In the text