© The Authors 2025

1 Introduction

The outer solar system harbours intriguing, and yet the least understood planets in the Solar System. Uranus and Neptune, often referred to as ‘ice giants’, remain enigmatic in terms of their composition and formation and evolution history (for example Helled 2025, and references therein). These two planets have always been key to understanding the solar system, but their importance has become even more apparent in recent years. Intermediate-mass exoplanets, with masses between Earth and Neptune, are found to be the most common planets in the galaxy (for example Zhu & Dong 2021). A better understanding of these giants could provide answers to several fundamental open problems in planetary physics. For example, it could reveal the mechanisms most relevant for the formation of a planet at a given distance from a specific star. It might help quantify the relationship between the observable atmospheres and the hidden interiors of a planet. And it could shed light on how a planet’s evolution after formation influences what we are able to observe today. To find answers to such problems, advanced planet formation and evolution models are needed. Such models can be improved by obtaining constraints from planetary interior models. In particular, it is important to ensure that such interior models are built with a minimal set of assumptions. They should provide an unbiased view regarding the possible internal structures of the planets. In this work, we present a novel and agnostic method to generate self-consistent planetary interior models by improving over previous work (Morf et al. 2024).

Interior models are designed to fit the measured physical properties of a planet. These correspond to the planetary mass M, equatorial radius R_eq, rotation period P_rot, and external (Newtonian) gravitational field V_ext: $$V_{\mathrm{e}\mathrm{x}\mathrm{t}}(\mathit{r})=\frac{GM}{r}(1-\sum _{n=1}^{\mathrm{\infty }}{\left(\frac{R_{\mathrm{r}\mathrm{e}\mathrm{f}}}{r}\right)}^{2n}{J}_{2n}{P}_{2n}(\cos \vartheta )),$$(1)

expressed by gravitational moments J_2n. We assume a northsouth and azimuthally symmetric planet described by spherical coordinates r = (r, ϑ, φ). The Legendre polynomials are denoted as P_2n and appear during the expansion of V using spherical harmonics (for example Zharkov & Trubitsyn 1978). The reference radius R_ref is an arbitrary normalisation constant that can be equal to R_eq, but does not have to be. The gravitational moments J_2n have only been measured for planets within our solar system, where Voyager 2 (for example French et al. 1988) provided the first data for Uranus and Neptune. These values have later been improved using moon and ring observations (Jacobson 2009, 2014; Wang et al. 2023; French et al. 2024; Jacobson & Park 2025). In a corotating reference frame, the total exterior potential U_ext is then given by: $$U_{\mathrm{e}\mathrm{x}\mathrm{t}}(\mathit{r})=V_{\mathrm{e}\mathrm{x}\mathrm{t}}(\mathit{r})+Q(\mathit{r})=V_{\mathrm{e}\mathrm{x}\mathrm{t}}(\mathit{r})+\frac{1}{2}{\left(\frac{2\pi }{P_{\mathrm{r}\mathrm{o}\mathrm{t}}}\right)}^2{r}^2{\sin }^2\vartheta ,$$(2)

assuming uniform rotation.

We briefly summarise recent progress in interior modelling of Uranus or Neptune. This includes work with so-called physical and empirical models. Physical models are self-consistent and provide a complete picture of a planet’s interior. They include a density, pressure, temperature, and composition profile compatible with a given Equation of State (EoS). However, physical models strongly depend on the assumptions used by the modeller. They only provide a single and rather biased possibility among many into a given planet’s interior. Empirical models a priori provide only the pressure and density profile within a planet. The interior temperature and composition can be inferred in a second step with an algorithm. Since density and pressure were first inferred without any EoS knowledge, such algorithms typically cannot adhere to the same standards as used for physical modelling. Empirical results for the temperature and composition hence suffer from inconsistencies or unphysical features (for example Neuenschwander et al. 2024; Morf et al. 2024). However, empirical models employ a minimal set of assumptions, enabling the modellers to achieve a more general view of the possible solutions for the planetary interior.

Recent work on physical interior models of Uranus or Neptune includes Militzer (2024); Cano Amoros et al. (2024); Tejada Arevalo (2025). Militzer (2024) base their models on ab-initio simulations for a mixture of water, methane, and ammonia. They show that these components phase separate under the pressure–temperature conditions in the interiors of Uranus and Neptune. The inferred interiors that match the observed gravity field and are compatible with magnetic field measurements are rock-poor: 2% and 7% of the total mass for Uranus and Neptune, respectively. Cano Amoros et al. (2024) provide hydrogen-water phase diagrams from available experimental and computational data. They find interior adiabatic structure models for Uranus and Neptune and infer a strong water depletion in the top layers of both planets. This is caused by the immiscibility of hydrogen and water at certain conditions. Tejada Arevalo (2025) present non-adiabatic and inhomogeneous evolution models for Uranus and Neptune assuming an interior composition of hydrogen, helium, methane, ammonia, water, and rocks. By setting the relevant initial parameters, they produce a Uranus model that preserves much of its primordial internal heat. In contrast, their Neptune model undergoes adiabatic cooling of the outer envelope, in agreement with observations.

Recent work on empirical interior models of Uranus or Neptune includes Movshovitz & Fortney (2022); Neuenschwander et al. (2024); Malamud et al. (2024); Morf et al. (2024); Lin et al. (2025); Mankovich et al. (2025). Movshovitz & Fortney (2022) investigate the constraining power of a high-precision orbiter measuring the gravitational fields of Uranus and Neptune. They find that while the ambiguity of Uranus’ and Neptune’s compositional profiles will still persist with high-precision data, determining the locations of any large composition gradient regions could become feasible. Neuenschwander et al. (2024) model the density of Uranus with up to three polytropes. Their models for Uranus’ interior include non-adiabatic regions. This leads to significantly hotter internal temperatures and higher rock-to-water ratios compared to physical models. Malamud et al. (2024) show that Uranus and Neptune could have accreted refractory-dominated planetesimals, while still remaining icy. They investigate chemical reactions of organic-rich refractory materials and the hydrogen in gaseous atmospheres of protoplanets. Uranus and Neptune could hence consist of large amounts of methane ‘ice’. Morf et al. (2024) provide a random algorithm to interpret empirical models of Uranus in terms of their temperature and composition. Additionally, they supply the Theory of Figures (Zharkov & Trubitsyn 1978) to tenth order, allowing for accurate calculations of higher-order gravitational moments. Their models in particular constrain the properties of the ionic water convection zone, likely responsible for Uranus’ magnetic field. Lin et al. (2025) argue that high degrees of mixing are required for Uranus interior models to be consistent with the measured gravitational moments. Furthermore, they show that the gravity measurements of a future orbiter could distinguish between high and low atmospheric metallicity scenarios. Mankovich et al. (2025) provide the spectra that could manifest in resonances with ring orbits or in Doppler imaging of Uranus’s visible surface. Their approach, originally developed for stellar oscillations, has the potential to significantly reduce the degeneracy in the set of otherwise allowed interior profiles.

The aim of this study is to develop methods that generate interior models that are both agnostic and physical, to bridge the gap between physical and empirical modelling. Our paper is organised as follows: In Section 2, we introduce our methods in more detail and give a general overview of our algorithm. In Sections 3 and 4, we show the inferred interior models of Uranus and Neptune, respectively. In Section 5, we elaborate and compare the inferred properties for Uranus and Neptune such as their hydrogen-helium abundances and rock-to-water ratios (Section 5.1). We investigate the models’ consistency with the observed magnetic fields (Section 5.2) and immiscibility constraints from hydrogen-helium-water mixtures (Section 5.3). A direct comparison of the Uranus versus Neptune models can be found in Section 5.4. We discuss limitations and the required follow-up work in Section 6 and present conclusions in Section 7.

2 Methods

Interior models are designed to fit the measured planetary data. In this study we focused on Uranus and Neptune and used the data listed in Table 1.

We present a global algorithm to generate agnostic and physical interior models. A schematic overview of the algorithm is shown in Figure 1. A more detailed description is presented in Figure B.1. The global algorithm starts by generating an initial random density profile ρ⁽⁰⁾(r): We assume a perfectly spherical planet with radius R and mass M. Spherical symmetry allows to describe the interior by a one-dimensional density function ρ = ρ(r), where r denotes the distance to the centre. Regarding the form of said function, every monotonically decreasing function f : r ∈ [0, R] → f(r) ∈ [1, 0] can be rescaled to a one-dimensional density profile ρ via $$\rho (r):=\frac{f(r)\cdot M}{4\pi \int r^{2}f(r)\mathrm{d}r}.$$(3)

We generated a random initial function f analogous to Podolak et al. (2022): Start defining f with the two tuples (r = 0, f(r) = 1) and (r = R, f(r) = 0). These two tuples are interpreted as corners of the rectangle [0, R] × [1, 0]. Now choose a random point (r*, f(r*)) in the rectangle, adding a third tuple to the definition of f. The point (r*, f(r*)) defines two new and smaller rectangles: [0, r*] × [1, f(r*)] and [r*, R] × [f(r*), 0]. Within these smaller rectangles, we can again choose new tuples, increasing the number of points where f is defined from three to five. Now repeat this procedure by choosing new tuples for each new rectangle as often as desired. This procedure ensures that f is indeed monotonically decreasing and generated in a random manner. For this work, we chose 1 + 2 + 4 + 8 + 16 = 31 tuples in addition to the two initial tuples and used linear interpolation between the tuples to get to N = 512 equidistantly spaced radial points with a corresponding value of f. Applying a Gaussian filter to f (as in Podolak et al. 2022) and using Equation (3) yielded the initial random density ρ⁽⁰⁾(r).

As a next step, ρ⁽⁰⁾(r) was given to the Theory of Figures¹ (Zharkov & Trubitsyn 1978; Morf et al. 2024). The Theory of Figures (ToF) takes the one-dimensional density function ρ(r), the mass M, the equatorial radius R_eq, and the rotation period P_rot as inputs and calculates the implied two-dimensional spheroidal shape of the planet. In particular, the gravitational moments J_n,ToF and the internal pressure P_ToF(r) are calculated as an output. The pressure P_ToF(r) is calculated by integrating the hydrostatic equilibrium equation: $$\mathbf{\nabla }P=\rho \mathbf{\nabla }U,$$(4)

where U is the total interior potential consisting of gravitational and centrifugal contributions (Appendix A.1 in Morf et al. 2024, for details).

Having obtained initial density ρ⁽⁰⁾(r) and pressure $P_{\text{ToF}}^{(0)}(r)$ profiles, the compositional algorithm is called next. The details of the compositional algorithm are discussed in Appendix B. In a nutshell, the compositional algorithm uses Equations of State (EoS) for different materials to infer a temperature $T_{\text{EoS}}^{(0)}(r)$ and composition profile ${\mathit{X}}_{\text{EoS}}^{(0)}(r)$ based on the given density-pressure profile. However, the densities implied by the EoS ${\rho }_{\text{EoS}}^{(0)}(r)$ can deviate from the originally provided densities ρ⁽⁰⁾(r). Although the compositional algorithm tries to keep the deviations as small as possible, they are necessary because a random density profile will almost never represent a physically realistic interior.

Due to these deviations, the density profile has to be renormalised (like in Equation (3)) to ensure that the planet still has the correct mass M. The renormalised density profile is denoted as ρ⁽¹⁾(r). Using the ToF a second time provides new values for the pressure $P_{\text{ToF}}^{(1)}(r)$ and gravitational moments $J_{n,\text{ToF}}^{(1)}$ that are consistent with ρ⁽¹⁾(r). At this point, one hence has a self-consistent tuple: $$({\rho }^{(1)}(r),P_{\mathrm{T}\mathrm{o}\mathrm{F}}^{(1)}(r),J_{n,\mathrm{T}\mathrm{o}\mathrm{F}}^{(1)}),$$(5)

and an inferred temperature and composition tuple: $$(T_{\mathrm{E}\mathrm{o}\mathrm{S}}^{(0)}(r),{\mathit{X}}_{\mathrm{E}\mathrm{o}\mathrm{S}}^{(0)}(r)).$$(6)

Combining Equations (5) and (6) to a full interior model is not self-consistent: The tuple in Equation (6) is consistent with ${\rho }_{\mathrm{E}\mathrm{o}\mathrm{S}}^{(0)}(r)\ne {\rho }^{(1)}(r)$ and $P_{\mathrm{T}\mathrm{o}\mathrm{F}}^{(0)}(r)\ne P_{\mathrm{T}\mathrm{o}\mathrm{F}}^{(1)}(r)$.

However, ρ⁽¹⁾(r) should be much more compatible with the EoS than ρ⁽⁰⁾(r) was initially. By calling the compositional algorithm again, one can expect that it can more easily infer temperature $T_{\mathrm{E}\mathrm{o}\mathrm{S}}^{(1)}(r)$ and composition ${\mathit{X}}_{\mathrm{E}\mathrm{o}\mathrm{S}}^{(1)}(r)$ profiles, whose implied densities ${\rho }_{\text{EoS}}^{(1)}(r)$ are close to the provided densities ρ⁽¹⁾(r). Consequently, ρ⁽²⁾(r) and ρ⁽¹⁾(r) should be much closer to each other compared to how close ρ⁽¹⁾(r) is to ρ⁽⁰⁾(r). Therefore, by repeating the above outlined procedure, eventually there should be a k ∈ ℕ, where the difference between ρ^(k+1)(r) and ρ^(k)(r) is so small for all r, that the tuple: $$({\rho }^{(k+1)}(r),P_{\mathrm{T}\mathrm{o}\mathrm{F}}^{(k+1)}(r),J_{n,\mathrm{T}\mathrm{o}\mathrm{F}}^{(k+1)},T_{\mathrm{E}\mathrm{o}\mathrm{S}}^{(k)}(r),{\mathit{X}}_{\mathrm{E}\mathrm{o}\mathrm{S}}^{(k)}(r))$$(7)

is practically self-consistent. Then, Equation (7) is a complete physical model of a planetary interior which has been generated in a random manner. We employed the criterion: $$\text{\hspace{0.17em}for all\hspace{0.17em}}r:\frac{|P_{\mathrm{T}\mathrm{o}\mathrm{F}}^{(k+1)}(r)-P_{\mathrm{T}\mathrm{o}\mathrm{F}}^{(k)}(r)|}{P_{\mathrm{T}\mathrm{o}\mathrm{F}}^{(k+1)}(r)}<ϵ$$(8)

to check whether a solution is self-consistent. Since the density and pressure are directly linked to each other via Equation (4), there is no need to define a second criterion for the density. In this work, we chose ϵ = 0.02. This is a conservative number given the uncertainties associated with the EoS of the different materials and their mixtures. For example, the assumption of ideal mixing (Equation (B.1)) alone can introduce deviations of up to 5% in density (Darafeyeu et al. 2024). A detailed discussion about algorithm convergence is given in Appendix C.

Finally, in general $J_{n,\text{ToF}}^{(k+1)}\ne J_{n,\text{ToF}}^{(0)}$ holds true at the end of the global algorithm. So there is no point in ensuring a precise match between the initial gravitational moments $J_{n,\text{ToF}}^{(0)}$ and the available data in Table 1. Only $J_{n,\text{ToF}}^{(k+1)}$ has to be data consistent, promoting Equation (7) from a physical interior model to a model also compatible with the data of the considered planet.

To summarise, we iterated until the density, pressure, temperature, and composition profile have converged and all agreed with the data in Table 1 and the used EoS. Our method is capable of generating interior models that are both agnostic and self-consistent. Being agnostic refers to using a minimal set of assumptions like for empirical models. However, we still maintain a level of self-consistency typical for physical models. As shown below, our results highlight that no single interior model can be regarded as the standard for a given planet. A wide range of self-consistent models must be considered to get a full picture that can reliably inform formation and evolution models.

Fig. 1 Overview of our global algorithm to infer agnostic and self-consistent planetary interior models. The difference between $\tilde{\rho }$ used by the compositional algorithm and ρ compatible with the interior pressure and planetary mass decreases until it vanishes using our iterative process. Traditionally, empirical models only use a single iteration, leading to inconsistencies. Figure B 1 contains more details.

Fig. 2 Inferred composition and convective and radiative structures for the four Uranus interior models. The legends summarise the total mass fractions for each component. The uppermost small arc of each slice is empty because we do not infer a composition for pressures below 100 bars.

3 Uranus

We present four Uranus models: U1, U2, U3, and U4. The inferred internal structures and compositions of the Uranus models are summarised in Figure 2. Models U1 and U3 correspond to water-rich interior solutions. Their respective rock-to-water ratios are low, just 0.04 and 0.08, consistent with Uranus being an ‘ice giant’. These models also possess similar total H-He abundances (14% and 11% of Uranus’ mass, respectively). However, their convection and temperature profiles vary significantly. Model U1 is mostly convective, has four different convection zones and reaches a central temperature below 13000 K. On the other hand, model U3 has two large stable regions. The innermost stable region covers more than half (54%) of Uranus’s radius and the central temperature is above 20 000 K.

On the other hand, models U2 and U4 are very similar in terms of their convection and thermal profiles. Both models reach a central temperature of ~6000 K. They are almost exclusively convective with two large convection zones separated by a small stable layer at ~75% of Uranus’ radius. However, the inferred compositions in models U2 and U4 are very different. In model U2, Uranus’ interior is rock-dominated with a rock-to-water ratio of 3.92, while the rock-to-water ratio in model U4 is 0.64. The H-He abundances are more comparable: 23% and 20% of Uranus’ mass for models U2 and U4, respectively. In addition, H-He is present also in the planetary deep interior in both cases. The H–He mass fractions for models U 2 and U 4 in the centre are 16% and 12%, respectively.

The density profiles and temperature-pressure profiles are presented in Figure 3. We find that the density profiles for models U2 and U4 are very similar. They both have a density discontinuity (Δρ > 1 g cm⁻³) at ~75% of Uranus’ radius and reach central densities of ~4.5 g cm⁻³. This density discontinuity around 3 · 10⁵ bar leads to a temperature jump that is 500 K higher for U2 compared to U4. Models U1 and U3 reach higher central densities (6.6 and 7.8 g cm⁻³, respectively) than U 2 and U 4. There are also no clear jumps in density in models U1 and U3. We further observe that all models leave the solution space discovered by Neuenschwander & Helled (2022) around a normalised radius of 0.8. This finding demonstrates how assumptions made by the modeller (such as the use of three consecutive polytropes in Neuenschwander & Helled 2022) can result in biased outcomes and exclude parts of the solution space. In contrast, our randomised approach represents a notable improvement, expanding the scope of the solution space.

For comparison, we show the measured gravitational moments from Jacobson & Park (2025) (in grey) alongside those from French et al. (2024) (in black) in Figure 3. Jacobson & Park (2025) only report formal (not realistic) uncertainties for their two results based on precession rates of either the ring centres of opacity (COO) or the geometric ring midlines (COR). They conclude that their knowledge of the gravitational moments is not better compared to that postulated by French et al. (2024).

Overall, our results show that both rock-dominated and water-dominated solutions are possible for Uranus. The same holds true for further properties: A wide variety of density, temperature, and convection profiles are feasible. The temperature and entropy profiles of our Uranus models can be found in Appendix A in Figure A.1. Relative pressure deviations according to Equation (8) are shown in Figure A.2.

Fig. 3 Density, pressure, and temperature profiles for the four Uranus models. Top: density versus normalised radius. For comparison, also shown are empirical solutions from Neuenschwander & Helled (2022) (grey area). The panel in the top right shows Uranus’ measured (French et al. 2024, black) and calculated (coloured) gravitational moments. We also show Uranus’ gravitational moments according to Jacobson & Park (2025) (grey) who find different results (COR vs COO) corresponding to different ring centre definitions. The coloured uncertainties depict estimated errors from the Theory of Figures. Bottom: temperature versus pressure. Also shown is the predicted region of ionic water inferred by Redmer et al. (2011). The dashed lines show the demixing line, below which hydrogen, helium, and water are expected to become immiscible (Howard et al. 2025).

4 Neptune

We present four Neptune models N1, N2, N3 and N4. Figure 4 shows sketches of the inferred structure and bulk composition of our Neptune models. Models N1 and N4 are consistent with the notion of Neptune being an ‘ice giant’. They have water-dominated interiors with rock-to-water ratios of 0.20 and 0.26, respectively, and H-He mass fractions of 12% and 10%, respectively. However, there are also significant differences. Model N1 has a central density higher than 18 g cm⁻³, while the central density of model N4 is only 6.2 g cm⁻³. Although both models include stable regions, they are more prevalent for model N1. For example, model N4 has a fully convective deep interior while the innermost region of model N1 is stable against convection and has a composition gradient. Model N4 initially has a larger temperature gradient compared to model N1, leading to a first appearance of ionic water at relatively larger normalised radii (0.89 vs 0.83) with higher temperatures (3100 K vs. 2500 K).

Models N2 and N3 represent Neptune models that are rock-dominated, supporting the idea of Neptune being a ‘rock giant’ (for example Helled & Fortney 2020; Teanby et al. 2020). Their rock-to-water ratios are found to be 1.78 and 1.72, with more than 50% of the planetary mass consisting of rocks. Although both models have large convective zones, the interior temperature profiles are rather different. Model N2 has a central temperature of 28 000 K, while the central temperature of model N3 is 12 000 K. Both of these inferred temperatures are significantly higher than the temperature inferred by standard adiabatic models which are of the order of 6000 K (for example Nettelmann et al. 2013). In addition, in model N3 H–He is found also in the planetary centre with 10% of the composition by mass. The existence of H-He in the deep interior is consistent with formation models of Neptune with composition gradients (for example Valletta & Helled 2022). The central region in model N2, on the other hand, is enriched with iron (18% by mass).

The inferred density profiles and pressure-temperature profiles are shown in Figure 5. For Neptune, we find four very unique density profiles. Model N1 has the highest central density of 18.4 g cm⁻³. The density in model N2 is below 3 g cm⁻³ in the outer half of the planet until it reaches a large density discontinuity (Δρ > 3 g cm⁻³). Models N3 and N4 have similar central densities (5.8 and 6.2 g cm⁻³, respectively), but differ with respect to the location of their density discontinuity. The discontinuities are located at a normalised radius of ~0.7 for model N3 and ~0.6 for model N2. We further find that models N1, N2, and N4 have a temperature jump larger than 500 K between pressures of 200 and 2000 bar. Also in the case of Neptune, there are significant differences in the density profiles in comparison to the solution space found by Neuenschwander & Helled (2022). This is observed at multiple locations (normalised radii of around 0.05, 0.4, and 0.75) with the differences being more notable compared to Uranus. Our randomised approach indeed covers a wider space of solutions compared to parametrised approaches.

We conclude that, like for Uranus, both rock-dominated and water-dominated solutions are valid, and that a wide range of compositions and internal structures is consistent with Neptune’s observed data. Our models consist of large convective zones or several stable regions, and can differ by more than 10 000 K with respect to their core temperatures. The temperature and entropy profiles of the four Neptune models can be found in Appendix A in Figure A.3. Relative pressure deviations according to Equation (8) are shown in Figure A.4.

Fig. 4 Same as Figure 2, but for Neptune.

5 Implications

5.1 Bulk composition

The inferred bulk compositions of our Uranus and Neptune models are summarised in Figure 6. We also list the total H-He mass fractions and compare them to previous work in Table 2. We find that the total H-He mass fractions (0.11–0.23 for Uranus and 0.10–0.19 for Neptune) are consistent with previous studies with the exception of Militzer (2024), who obtains lower values. The derived rock-to-water ratios of our models (and selected previous studies) are listed in Table 3. We find that for both Uranus and Neptune, the rock-to-water ratio is poorly constrained and ranges between 0.04 and 3.92 for Uranus, and between 0.20 and 1.78 for Neptune. Overall, our models include solutions with rock-to-water ratios that are significantly higher compared to previous studies (up to 3.92 for Uranus and 1.78 for Neptune). At the same time, we show that solutions with low rock-to-water ratios (as low as 0.04 for Uranus and 0.20 for Neptune) are also valid. We can therefore conclude that given the current data, there is no clear reason to prefer water-rich or rock-rich interiors for Uranus and Neptune. Therefore, referring to Uranus and Neptune as ‘ice giants’ may be inappropriate since the bulk composition of the planets remain unclear (for example Helled & Fortney 2020; Hofstadter et al. 2024). Constraining the rock-to-water ratio would require additional (and more accurate) measurements (for example Nimmo et al. 2024) or additional constraints from planetary formation models (for example Helled 2025).

Fig. 5 Same as Figure 3, but for Neptune.

Fig. 6 Inferred bulk composition for Uranus and Neptune.

5.2 Magnetic field generation

Thanks to Voyager 2, we know that both Uranus and Neptune possess considerable magnetic fields (Connerney et al. 1987; Ness et al. 1989). The fields are multipolar and tilted significantly relative to the rotation axes of the planets (for example Jones 2011). These observations are relevant when modelling the planetary interiors, as solutions have to be consistent with dynamo generation. A common explanation for the magnetic fields of the planets is the presence of ionic water in a convective region (Redmer et al. 2011). For example, Stanley & Bloxham (2004, 2006) found that a thin convective and conductive shell in the outer planetary region can suppress large-scale dipolar modes and hence explain the observed multipolar fields. This occurs for both liquid (+stably stratified) and solid (+insulating) cores. And Gastine et al. (2012) showed that multipolar fields can be reproduced for an even wider range of scenarios.

Figure 7 shows the extent of the convective and stable regions, as well as the pressure-temperature region where (pure) water is expected to be ionic. We find that all the models indeed possess a thin layer of electrically conducting ionic water that is convective. We list the upper (r₂) and lower (r₁) limits of these dynamo regions in Table 4 in units of normalised radii and compare our results to previous work. For Uranus, we find that our dynamo regions are located deeper in the interior compared to Morf et al. (2024). This can be explained by the less steep temperature gradients in this work, leading to a later onset of the ionic phase transition. Our results for r₂ (0.69–0.74 for Uranus and 0.78–0.92 for Neptune) are in agreement with Militzer (2024). However, r₁ is found to be larger: 0.62–0.67 for Uranus and 0.67–0.86 for Neptune. This in turn leads a higher r₁/r₂ ratio: 0.86–0.97 for Uranus and 0.82–0.93 for Neptune, compared to Militzer (2024) and the preferred scenarios in Stanley & Bloxham (2004, 2006). However, our findings for r₁ should be interpreted as upper bound values. This is because the region where the magnetic field is generated can change when considering different composition components. Here we only consider pure water and follow the phase diagram of Redmer et al. (2011). In principle, other mixtures (such as a mixture of water and rock) could enable dynamo generation (for example Huang et al. 2020; Gao et al. 2022). In that case our magnetic dynamo regions could extend deeper into the interior than indicated in Table 4. This would bring our ratio of r₁/r₂ closer to Stanley & Bloxham (2004, 2006); Militzer (2024). We note that a mixture of rocks with hydrogen or mixtures of other volatiles such as ammonia and methane could also be electronically conducting. Furthermore, other studies (for example Soderlund et al. 2013) found multipolar dynamos while considering different geometries and despite the presence of large, solid, electrically conducting inner cores. A better understanding of planetary dynamos as well as material properties at planetary conditions for different elements and mixtures is clearly desirable.

Fig. 7 Convective and stable regions in our models for Uranus and Neptune. Also shown are regions where ionic water is expected.

5.3 Immiscibility of hydrogen, helium and water mixtures

Our models are based on the ideal mixing approximation (Appendix B and Equation (B.1) contain more details). This assumption can lead to differences in density of the order of a few percent (for example Darafeyeu et al. 2024). In particular, under this assumption immiscibility between different elements is not considered. Howard et al. (2025) provide miscibility results based on ab initio calculations of the hydrogen-water phase diagram. However, since our models include a mixture of hydrogen, helium, and water, a correction is needed. Using their Equation (A.1), we can calculate the pressure-temperature line, below which hydrogen, helium, and water become immiscible. We employ: $$x_{{\mathrm{H}}_2\mathrm{O}}=\frac{Z_{{\mathrm{H}}_2\mathrm{O}}(1+\frac{\alpha }{2})}{9-Z_{{\mathrm{H}}_2\mathrm{O}}(9-(1+\frac{\alpha }{2}))},$$(9)

when calculating the number fraction of water x_H2O from the mass fraction of water Z_H2O. The variable α = 0.705/0.275 denotes the proto-solar mass ratio between hydrogen and helium. Equation (9) ensures that x_H2O + x_H2 = 1 − x_He, m_H2O/m_H2 = 9, and m_He/m_H2 = 2. The inferred demixing lines are presented with dashed lines in Figures 3 and 5. We find that in all models, the P − T lines remain above the demixing curves. We can therefore conclude that our models are valid, and no demixing is expected.

5.4 Uranus versus Neptune

Uranus and Neptune are often thought of as being similar, but it remains unclear how alike these planets truly are (for example Helled & Fortney 2020). There are obvious differences between the planets (equatorial radii, masses, rotation periods, gravity fields), which are listed in Table 1.

In the outermost convective layers, we tend to find higher H-He mass fractions for Uranus (0.62–0.73) compared to Neptune (0.25–0.49). This is consistent with the findings of Nettelmann et al. (2013); Cano Amoros et al. (2024); Tejada Arevalo (2025) and could indicate that Uranus indeed possesses more H-He in the outermost regions compared to Neptune. Furthermore, the outermost edge r₂ (Table 4) of the magnetic dynamo regions tend to be slightly higher up for Neptune (0.78–0.92) compared to Uranus (0.69–0.74). This is consistent with the findings of Militzer (2024) and therefore also likely indicative of a dichotomy between Uranus and Neptune. Atmospheric probes sent to the atmospheres of both planets could test these predictions. Regarding the general bulk compositions (Figure 6) or the inferred rock-to-water ratios (Table 3), we remain agnostic about inferring concrete differences between Uranus and Neptune. We can only confirm that a wide range of options are possible for both planets.

We further note that the two planets have different measured heat fluxes (Pearl & Conrath 1991) which is often explained with different structures such as gradients or boundary layers (for example Vazan & Helled 2020; Scheibe et al. 2019, 2021) or giant impacts (for example Reinhardt et al. 2020). Our models do not include energy balance or flux calculations. If differences in structures are indeed the reason for the different heat fluxes, that would favour models U1 and U3 for Uranus and models N2 and N3 for Neptune. Models U1 and U3 contain extensive composition gradients where heat transport cannot be facilitated efficiently via convection. As a consequence, comparatively little heat would be able to escape Uranus, while the large convection zones in models N2 and N3 would enable a more efficient heat transport and therefore a larger heat flux in comparison to Uranus. Interestingly, models U1 and U3 are rock-poor, while models N2 and N3 are rock-rich. While this could be indicative of a bulk composition dichotomy between Uranus and Neptune, we refrain from making any definite statements. For a robust conclusion, both more models and heat flux calculations (and measurements) would be required. We hence remain agnostic about this issue and hope to investigate this topic in detail in future work.

6 Discussion

This work introduces a novel method for modelling the internal structures of planetary interiors. By uniting the flexibility of empirical models with the physical validity of traditional approaches, we produce agnostic and self-consistent interior solutions. Our iterative algorithm allows for the convergence of randomly initialised density profiles. They converge towards solutions that simultaneously satisfy hydrostatic equilibrium, fit the gravity data, and respect constraints imposed by the employed EoS. We note that interior models can be further constrained by seismology Mankovich et al. (2025). We provide realistic Brunt–Väisälä frequencies necessary for interpreting oscillation modes in planetary interiors in Figures A.1 and A.3. Our models are well suited for comparison with asteroseismic observations by future missions. While this study represents a major step forward in planetary modelling, some limitations remain and future research is clearly desirable.

First, the EoS of materials at planetary conditions themselves may contain errors that can affect the results (More et al. 1988; Chabrier & Debras 2021; Haldemann et al. 2020; Cozza et al. 2025). For example, the calculations in Haldemann et al. (2020) contain an error that affects the entropy, but not its derivatives (Aguichine et al. 2025). Tejada Arevalo (2025) find that correcting for this error produces qualitatively similar results. All EoSs of pure materials possess uncertainties and it is hence clear that the choice of the EoS affects the results. In Appendix D, we compare two different EoS for SiO₂ to quantify how different their results are as an additional example. We find that the differences in density can be of the order of a few percent. Furthermore, the assumption of ideal mixing (Equation (B.1)) can yet again lead to differences in density of the order of a few percent (for example Darafeyeu et al. 2024) and therefore affect the inferred composition. So even if the EoSs of the pure materials would serve as perfect constraints, a complexity arises from the fact that the mixtures are treated in an over-simplified manner. Finally, we only considered immiscibility for the hydrogen-helium-water mixture, but did not include other immiscibilities (like a mixture of water and rock). Overall, we can therefore expect that the uncertainties associated with all of the aforementioned issues would also affect the inferred bulk compositions.

Second, for the elements included in Uranus and Neptune, we only considered hydrogen, helium, water, rocks, and iron. Other elements such as methane or ammonia are also expected to exist in the planetary interiors and could also influence the results (for example Malamud et al. 2024; Tejada Arevalo 2025). For example, methane and ammonia, being less dense than water and rock, can contribute to the reduced densities required by gravitational harmonic constraints. Consequently, our models may over-predict the hydrogen and helium abundance. Including methane and ammonia could also elevate the temperature profiles due to these lower densities, assuming a constant specific entropy. In this work we do not include methane and ammonia since the exact abundance ratios between different elements remain uncertain. In other words, one would either have to rely on an educated guess (Tejada Arevalo 2025), or add even more free parameters to the problem.

Third, our models rely on the pressure-density relation provided in Hueso et al. (2020) for pressures below 100 bars. We rely on an atmospheric model because our approach is not well suited to handle the complex chemical processes occurring in planetary atmospheres such as condensation. Atmospheric models are generally more detailed and therefore more accurate in representing the conditions in the planetary atmosphere. Therefore, we rather not assume an atmospheric composition that could bias the results of the interior model. However, we do ensure that our models have density-pressure relations that are compatible with these more accurate atmospheric models (Morf et al. 2024, for details). Although the used atmosphere model by Hueso et al. (2020) should be more accurate than our approach below 100 bar, it also includes uncertainties. To account for this, we allowed the initial temperature $T_{\text{Eos}}^0$ to deviate from the atmosphere model. Nonetheless, a different atmosphere model could lead to different results.

Fourth, we did not consider dynamical wind corrections to the measured gravitational moments (Kaspi et al. 2013; Soyuer et al. 2023). This can affect the density profile solution space and therefore our findings. However, Neuenschwander et al. (2024) investigate this effect for Uranus and suggest that wind corrections should not have a large influence on the inferred temperature and composition profiles. Since Uranus’ gravity data are significantly more constrained than Neptune’s, we expect that the same holds true for Neptune.

Finally, the rotation periods of Uranus and Neptune might deviate from the rotation periods listed in Table 1 (Helled et al. 2010). Changing the rotation periods has a large impact on the density profile solutions space and the inferred composition (for example Nettelmann et al. 2013; Neuenschwander et al. 2024). For this work, we rely on purely observational constraints (Desch et al. 1986; Karkoschka 2011) as there is no broader consensus or further validation of the alternative estimates by Helled et al. (2010). We acknowledge the importance of this issue and plan to explore it in future work.

7 Conclusions

The key conclusions of our study can be summarised as follows:

• Bridging empirical and physical interior modelling:

  By using both random density profiles and a random compositional algorithm, we minimise bias for our planetary interior models. At the same time, we ensure a level of self-consistency typically reserved for physical interior models by using a novel iterative procedure. Our results thus cover a wide range of options for interior density, temperature, and composition profiles that are physical and robust.

• Uranus and Neptune as ‘rock giants’:

  We find that different bulk compositions are consistent with current data for both Uranus and Neptune. They encompass rock-dominated and water-rich interiors, with rock-to-water ratios varying by more than an order of magnitude (0.04–3.92 for Uranus, 0.20–1.78 for Neptune; Table 3). We hence confirm that the label ‘ice giants’ for Uranus and Neptune may be more of a historical artifact rather than a robust physical classification.

• Differences between Uranus and Neptune:

  In the outermost convective zones, we find a higher local H-He mass fraction for Uranus (0.62–0.73) than for Neptune (0.25–0.49). Furthermore, it seems that Uranus’ magnetic field is generated deeper in the planetary interior compared to Neptune. We find upper limits of 0.69–0.74 (Uranus) vs. 0.78–0.92 (Neptune) for the outer edges of the dynamo regions in units of normalised radii (Table 4).

• Dynamo and immiscibility constraints matched:

  All our models contain convective layers where (pure) water exists in its ionic phase. Such regions could serve as the dynamo source of Uranus’ and Neptune’s multipolar magnetic fields. Furthermore, none of the modelled tempera-ture-pressure profiles cross the hydrogen–water immiscibility boundaries derived from recent ab-initio calculations.

With the potential for future dedicated missions to Uranus and Neptune, our method also provides a flexible and unbiased tool for interpreting forthcoming data. Ultimately, the interiors of Uranus and Neptune remain enigmatic, not because they are beyond reach, but because the data required to resolve their secrets are still out of grasp. Until then, only a plurality of models, not a singular one, can capture the full extent of possibilities for their hidden depths.

Acknowledgements

We thank the anonymous referee for his valuable suggestions that helped strengthen this work. We also would like to thank Simon Müller and Christian Reinhardt for providing us with numerical implementations of the employed EoS. This work was supported by the Swiss National Science Foundation (SNSF) through a grant provided as a part of project number 215634: https://data.snf.ch/grants/grant/215634

References

Appendix A Additional interior data

We show interior temperature, entropy, and frequency profiles for the models of Uranus and Neptune in Figures A.1 and A.3. We find that, as expected, the Brunt-Väisälä frequency N is nearly zero (|N²| < 10⁻⁸ h⁻²) in adiabatic regions. This is orders of magnitudes smaller than analogous results in Morf et al. (2024). The frequency N² should be smaller than zero throughout all adiabatic regions subjected to convection. This is violated in the outermost parts (normalised radii above 0.9) in Figures A.1 and A.3. However, these are caused by numerical floating point precision limitations with near-zero numbers, and they do not contain any physical meaning. The wiggles visible in Ledouxstable regions (N² > 10⁻² h⁻²) originate from the random nature of the algorithm.

We show the relative pressure deviation criterion introduced in Equation 8 in Figure A.2 for Uranus and Figure A.4 for Neptune. We find that all the models are consistent with the ϵ = 0.02 limit introduced in Section 2. In this work, the number k in Equation 8 ranges from 3 to 8.

Fig. A.1 Brunt-Väisälä frequencies (top), entropies (middle), and temperature (bottom) profiles of the four Uranus models. For comparison, we also show purely adiabatic temperature profiles from Nettelmann et al. (2013) (dashed grey).

Fig. A.2 Relative pressure deviations (Equation 8) of the four Uranus models.

Fig. A.3 Same as Figure A.1, but for Neptune.

Fig. A.4 Same as Figure A.2, but for Neptune.

Appendix B Compositional algorithm

We now discuss the compositional algorithm in detail, a part of the global algorithm depicted in Figure B.1. The compositional algorithm is depicted in Figure B.2 and represents an improved version of the work presented in Morf et al. (2024). As input, it requires both density values [ρ⁰,. . ., ρ^N] and pressure values [P⁰,. . ., P^N] within the planet. The index 0 denotes the outermost layer, where N corresponds to the planetary centre. In this work, we defined the first layer that has a pressure above 100 bars as the outermost layer. For layers below a pressure of 100 bar, we did not infer a composition and just fixed the profile to follow the density-pressure atmosphere models proposed by Hueso et al. (2020). This is motivated by more complex chemical processes in the atmosphere (such as in Hueso et al. 2020), which our approach is not well suited for.

We considered hydrogen (X) and helium (Y), water (Z₁), rocks (Z₂), and iron (Z₃) for a composition mass fraction X_Eos = (X, Y, Z₁, Z₂, Z₃). To that end, we relied on the EoS provided by More et al. (1988); Chabrier & Debras (2021); Haldemann et al. (2020). The abundance of hydrogen and helium relative to each other was always fixed at a proto-solar mass ratio of 0.705/0.275. We assumed the ideal mixing equation $$\frac{1}{\rho }=\frac{X}{{\rho }_X}+\frac{Y}{{\rho }_Y}+\frac{Z_1}{{\rho }_{Z_1}}+\frac{Z_2}{{\rho }_{Z_2}}+\frac{Z_3}{{\rho }_{Z_3}},$$(B.1)

to combine multiple different materials and their respective EoS.

In the beginning, the compositional algorithm chooses a random temperature $T_{\text{Eos}}^0$. For this work, we randomly chose $T_{\text{Eos}}^0$ ∈ [250, 450] K for the first layer that has a pressure above 100 bars. This temperature range is based on the results of Hueso et al. (2020). Given $T_{\text{Eos}}^0$, the compositional algorithm tries to find a composition ${\mathit{X}}_{\text{EoS}}^0$ that results in a density ${\rho }_{\text{EoS}}^0$ that is close enough to the provided density ρ⁰, meaning $$\frac{|{\rho }^0-{\rho }_{\mathrm{E}\mathrm{o}\mathrm{S}}^0|}{{\rho }^0}<\delta .$$(B.2)

Here, we used δ = 0.01. If multiple such ${\rho }_{\text{EoS}}^0$ exist, the algorithm chooses one of them randomly.

Moving towards the planetary interior, the compositional algorithm prioritises finding regions of constant composition. These regions are adiabatic due to large scale convection mixing the materials homogeneously. To find them, the compositional algorithm checks whether it can keep the composition ${\mathit{X}}_{\text{EoS}}^{\ast }$ from the previous layer and follow an adiabatic gradient for the next Δ layers. The check is successful if the densities inferred by the EoS obey Equation B.2 for all Δ layers. We only accepted such regions if Δ > N/8 = 64. This limit was introduced to prevent small convection cells. While small convection cells are possible, sustaining them over long enough time scales is tricky and highly dependent on dynamical processes such as oscillatory double-diffusive convection (Leconte & Chabrier 2012; Tulekeyev et al. 2024). We do not investigate such dynamical processes better treated by evolution models (such as Tejada Arevalo 2025; Eberlein & Helled 2025), which motivated the aforementioned limit.

Fig. B.1 Global algorithm with more details than Figure 1. As input, it requires a density profile ρ(0)(r) of a planetary interior. In this work, ρ(0)(r) is generated in a random manner. The global algorithm then repeatedly calls the Theory of Figures (ToF) and a compositional algorithm to infer the pressure, temperature, and composition profile of the planet. The algorithm is applied until everything is self-consistent up to a certain tolerance. The algorithm is then deemed successful if the gravitational moments calculated by the ToF are compatible with the provided data.

Fig. B.2 Compositional algorithm. As input, it requires density values [ρ0,. . .,ρN] and pressure values [P0,. . .,PN] for the planetary interior. If successful, the algorithm infers a temperature and composition profile. The densities implied by the Equations of State (EoS) can deviate from the input densities [ρ0,. . .,ρN].

As a second priority, the compositional algorithm checks if it can find a homogeneous and adiabatic region by ignoring Equation B.2 for Δ − 1 layers. Only satisfying Equation B.2 again for the Δ-th layer violates self-consistency to a high degree. However, due to the iterative nature of the global algorithm, the necessity for such violations should become less and less frequent.

Finally, as a third priority, the algorithm allows the composition to change: In addition to the inferred ρ_EoS satisfying Equation B.2, one now needs to check for stability against convection. This is done by calculating the Brunt-Väisälä frequency N. For an arbitrary density ρ, pressure P, and temperature T, N is given by (for example Unno et al. 1989; Kippenhahn et al. 2013) $$N^2=\frac{g^2\rho }{P}\frac{{\chi }_T}{{\chi }_{\rho }}({\mathrm{\nabla }}_{\mathrm{a}\mathrm{d}}-{\mathrm{\nabla }}_T+B),$$(B.3)

where g is the gravity acceleration. The thermodynamic quantities χ_ρ and χ_T are defined as $${\chi }_{\rho }={\left(\frac{\mathrm{\partial }\ln P}{\mathrm{\partial }\ln \rho }\right)}_{T,\mu },{\chi }_T={\left(\frac{\mathrm{\partial }\ln P}{\mathrm{\partial }\ln T}\right)}_{\rho ,\mu },$$(B.4)

where the molecular weight μ parametrises the composition. In particular, N depends on the difference between the adiabatic gradient ∇_ad and temperature gradient ∇_T $${\mathrm{\nabla }}_{\mathrm{a}\mathrm{d}}:={\left(\frac{\mathrm{\partial }\ln T}{\mathrm{\partial }\ln P}\right)}_S,{\mathrm{\nabla }}_T:=\frac{d\ln T}{d\ln P},$$(B.5)

as well as on the mixing term B $$B=\frac{{\chi }_{\rho }}{{\chi }_T}\frac{\ln \rho (P,T,\mu +\frac{d\mu }{d\ln \rho }\delta \ln \rho )-\ln \rho (P,T,\mu )}{\delta \ln \rho }\frac{d\ln \rho }{d\ln P},$$(B.6)

taken in the limit δ ln ρ → 0. The mixing term B can be approximated numerically (Appendix B of Morf et al. 2024). The Brunt-Väisälä frequency N denotes the oscillation frequency of a parcel of material that is slightly displaced compared to its surroundings while not exchanging any heat with it. For stability, one requires N² ≥ 0. The criterion N² ≥ 0 is called the Ledoux criterion. If multiple tuples (${\rho }_{\text{EoS}}^i,T_{\text{EoS}}^i,{\mathit{X}}_{\text{EoS}}^i$) satisfy both Equation B.2 and N² ≥ 0, the algorithm chooses one of them randomly.

The compositional algorithm succeeds by repeatedly applying the above steps in accordance with their priorities for the whole planet. The compositional algorithm fails if at some point none of the above steps are possible. If successful, it provides a temperature and composition profile for the planet. We stress again that the densities implied by the EoS can deviate from the input densities, sometimes even more than allowed by Equation B.2. The inferred temperature and composition profiles are hence only physical, if the entire model (Equation 7) satisfies the relative pressure deviation criterion stated in Equation 8.

Appendix C Algorithm convergence

In this section we investigate the convergence properties of the global algorithm. Figure C.1 shows the maximal deviation in the gravitational moments $$\underset{\text{\hspace{0.17em}moment deviation\hspace{0.17em}}}{\mathrm{maximal}\text{\hspace{0.17em}gravitational\hspace{0.17em}}}=\underset{i\in \{2,4\}}{max}\frac{{(J_{i,\text{ToF}}^{(k+1)}-J_{i,\text{data}})}^2}{{\sigma }_i^2},$$(C.1)

Fig. C.1 Maximal gravitational moment deviation versus maximal relative pressure deviation for different generations k ∈ {1, 3, 5} and in total over ten thousand Neptune models. Histograms of the depicted points are included for both axes. Points falling within the hatched area represent physical models. The arrows highlight one single path between generations.

versus the maximal relative pressure deviation $$\underset{\text{\hspace{0.17em}pressure deviation\hspace{0.17em}}}{\mathrm{maximal}\text{\hspace{0.17em}relative\hspace{0.17em}}}=\underset{r}{max}\frac{|P_{\mathrm{T}\mathrm{o}\mathrm{F}}^{(k+1)}(r)-P_{\mathrm{T}\mathrm{o}\mathrm{F}}^{(k)}(r)|}{P_{\mathrm{T}\mathrm{o}\mathrm{F}}^{(k+1)}(r)},$$(C.2)

for different generations k ∈ {1, 3, 5}. Each data point represents one model, a complete tuple according to Equation 7. Points within the hatched area are both consistent with the ϵ = 0.02 criterion (Equation 8) and with the measured gravitational moments within their respective uncertainties (Table 1). Such points therefore correspond to physical models.

Figure C.1 highlights two properties of our algorithm: First, the histograms show that as we progress to higher generations, our models become more physical. They align more with the measured gravitational moments and also become more self-consistent according to Equation 8. Second, considering the single highlighted path (black arrows) between the generations, the random nature of our compositional algorithm becomes apparent. While the models become more physical on average, that must not be true for any given model. The first model (blue point with black circle, k = 1) has pressure deviations of almost 100% and gravitational moments more than ten standard deviations away from the measured values. After just two generations (orange point with black circle, k = 3), it leads to a physical model thanks to random choices made by the compositional algorithm. However, two additional generations later (green point with black circle, k = 5), the result is again far away from the measured gravitational moments.

We can theretofore conclude that while the non-random parts of the algorithm lead towards physical models on average, the random nature of other algorithm parts ensure that a wide variety of options are explored. The aforementioned parameter ϵ should be thought of as a selection criterion, not a convergence criterion. There is no reason to prefer, for example, a generation k = 5 models over a generation k = 3 model. Higher generations are not necessarily refined versions of what came before. Finally, Figure C.1 demonstrates that one could easily choose lower values for ϵ. However, as discussed earlier, lower values for ϵ would not translate to a closer approach to physical reality. Each limitation listed in Section 6 is capable of introducing deviations comparable to or easily surpassing our ϵ = 0.02 threshold.

Appendix D Comparing different EoS for SiO₂

Throughout this work, we have emphasised that the EoS employed carry inherent uncertainties, which can lead to discrepancies when compared with alternative EoS formulations. To illustrate this, Figure D.1 presents a comparison focused solely on the SiO₂ EoS, where we contrast the EoS used in this study (More et al. 1988) with the ANEOS formulation for forsterite (Stewart et al. 2019), constructed using M-ANEOS which is based on fitting analytic expressions of the Helmholtz free energy (Thompson et al. 1974; Melosh 2007; Thompson et al. 2019).

Fig. D.1 Pressure versus density of SiO2 for our Uranus and Neptune models. The solid lines correspond to the EoS of More et al. (1988), while the dashed lines to the ANEOS EoS. The density of SiO2 (ρZ2) is not equal to the model density (ρ) since there are no regions with pure rock in the models (Equation B.1). We only show the regions where SiO2 is present in our models.

We find that the differences in density between the two SiO₂ EoS are generally below 10%, with many models showing even smaller deviations. These differences only affect the SiO₂ component in our ideal mixing formulation (Equation B.1). When comparing pressures at constant density instead by inverting the relation, discrepancies in pressure can reach up to 50%. Our adopted consistency threshold of ϵ = 0.02 in Equation 8 is clearly conservative in comparison. Finally, while both EoSs are developed for pure SiO₂, the behaviour of materials in mixtures can differ significantly (immiscibility, phase transitions, etc.). It is therefore desirable to perform detailed simulations of different mixtures at planetary conditions. Currently, there is no clear justification for preferring one EoS over the other, especially for materials ’heavier’ than hydrogen and helium.

All Tables

All Figures

Fig. 1 Overview of our global algorithm to infer agnostic and self-consistent planetary interior models. The difference between $\tilde{\rho }$ used by the compositional algorithm and ρ compatible with the interior pressure and planetary mass decreases until it vanishes using our iterative process. Traditionally, empirical models only use a single iteration, leading to inconsistencies. Figure B 1 contains more details.

In the text

	Fig. 2 Inferred composition and convective and radiative structures for the four Uranus interior models. The legends summarise the total mass fractions for each component. The uppermost small arc of each slice is empty because we do not infer a composition for pressures below 100 bars.
In the text

Fig. 3 Density, pressure, and temperature profiles for the four Uranus models. Top: density versus normalised radius. For comparison, also shown are empirical solutions from Neuenschwander & Helled (2022) (grey area). The panel in the top right shows Uranus’ measured (French et al. 2024, black) and calculated (coloured) gravitational moments. We also show Uranus’ gravitational moments according to Jacobson & Park (2025) (grey) who find different results (COR vs COO) corresponding to different ring centre definitions. The coloured uncertainties depict estimated errors from the Theory of Figures. Bottom: temperature versus pressure. Also shown is the predicted region of ionic water inferred by Redmer et al. (2011). The dashed lines show the demixing line, below which hydrogen, helium, and water are expected to become immiscible (Howard et al. 2025).

In the text

	Fig. 4 Same as Figure 2, but for Neptune.
In the text

	Fig. 5 Same as Figure 3, but for Neptune.
In the text

	Fig. 6 Inferred bulk composition for Uranus and Neptune.
In the text

	Fig. 7 Convective and stable regions in our models for Uranus and Neptune. Also shown are regions where ionic water is expected.
In the text

	Fig. A.1 Brunt-Väisälä frequencies (top), entropies (middle), and temperature (bottom) profiles of the four Uranus models. For comparison, we also show purely adiabatic temperature profiles from Nettelmann et al. (2013) (dashed grey).
In the text

	Fig. A.2 Relative pressure deviations (Equation 8) of the four Uranus models.
In the text

	Fig. A.3 Same as Figure A.1, but for Neptune.
In the text

	Fig. A.4 Same as Figure A.2, but for Neptune.
In the text

Fig. B.1 Global algorithm with more details than Figure 1. As input, it requires a density profile ρ(0)(r) of a planetary interior. In this work, ρ(0)(r) is generated in a random manner. The global algorithm then repeatedly calls the Theory of Figures (ToF) and a compositional algorithm to infer the pressure, temperature, and composition profile of the planet. The algorithm is applied until everything is self-consistent up to a certain tolerance. The algorithm is then deemed successful if the gravitational moments calculated by the ToF are compatible with the provided data.

In the text

	Fig. B.2 Compositional algorithm. As input, it requires density values [ρ0,. . .,ρN] and pressure values [P0,. . .,PN] for the planetary interior. If successful, the algorithm infers a temperature and composition profile. The densities implied by the Equations of State (EoS) can deviate from the input densities [ρ0,. . .,ρN].
In the text

	Fig. C.1 Maximal gravitational moment deviation versus maximal relative pressure deviation for different generations k ∈ {1, 3, 5} and in total over ten thousand Neptune models. Histograms of the depicted points are included for both axes. Points falling within the hatched area represent physical models. The arrows highlight one single path between generations.
In the text

	Fig. D.1 Pressure versus density of SiO2 for our Uranus and Neptune models. The solid lines correspond to the EoS of More et al. (1988), while the dashed lines to the ANEOS EoS. The density of SiO2 (ρZ2) is not equal to the model density (ρ) since there are no regions with pure rock in the models (Equation B.1). We only show the regions where SiO2 is present in our models.
In the text