© The Authors 2024

1 Introduction

Uranus has only been visited by the Voyager 2 spacecraft. This Uranus flyby yielded valuable gravity data (for example French et al. 1988) that are used to constrain the planetary interior. However, since it was a single flyby, only the first two gravitational moments (J₂ and J₄) were measured with rather high uncertainties. These values have been later improved using observations of Uranus’ moons and rings (for example Jacobson 2014; French et al. 2024). Nonetheless, Uranus’ gravity data remain insufficient to constrain its internal structure and bulk composition.

At the moment, various key open questions regarding Uranus’ nature remain open. For example, the name ‘ice giant’ may not be appropriate for Uranus, since it remains unknown whether Uranus’ bulk composition is water- or rock-dominated (Helled et al. 2020). Another unsolved mystery is the origin of Uranus’ measured low heat flux. One possibility is that Uranus’ deep interior is still very hot and that the primordial heat cannot escape due to thermal boundary layers or non-convective regions caused by composition gradients (for example Nettelmann et al. 2016; Vazan & Helled 2020; Scheibe et al. 2019, 2021). Alternatively, it is possible that Uranus’ low luminosity corresponds to a cold interior, which could be a result of a giant impact that led to rapid cooling (Reinhardt et al. (2020)). Improving our understanding of Uranus’ interior and composition could aid in testing and constraining planet formation and evolution models.

Interior models are designed to fit the measured physical properties of the planets. These correspond to the planetary mass, radius, rotation rate, and gravitational field (gravitational moments). There are two widely used methods to model the planetary internal structure. The first corresponds to the standard approach, the so-called physical interior models, where the structure equations are solved together with an Equation of State (EoS) of the materials assumed by the modeller. This includes assuming a heat transport mechanism within the planet to get a temperature profile of the planetary interior. In the simplest case, the temperature profile is assumed to be adiabatic corresponding to a fully convective interior. Further details on interior modelling can also be found in Helled & Fortney (2020) and the references therein.

The second method, the so-called empirical interior models, infers density profiles compatible with the measurements without making an a piori assumption on the planetary composition and temperature profile. Density profiles are typically represented by polytropes (for example Neuenschwander et al. 2021; Neuenschwander & Helled 2022), polynomials (for example Helled et al. 2011; Movshovitz & Fortney 2022), or by an arbitrary function (for example Marley et al. 1995; Podolak et al. 2000, 2022). The pressure is inferred in a second step assuming hydrostatic equilibrium. Since no assumptions are made for the materials and heat transport mechanism, the composition and the temperature profile of the planet remain unknown. This approach finds significantly more density profiles than physical models, but may also include nonphysical ones. For example, empirical density profiles can posses a central density that is too high for any plausible material to sustain.

In general, determining the composition and internal structure of a planet is challenging due to the degenerate nature of this problem. Even for Jupiter, whose planetary mass, radius, rotation rate, gravitational field, and atmospheric composition have been accurately measured, the data are insufficient to uniquely determine Jupiter’s interior (for example Miguel et al. 2022; Howard et al. 2023; Militzer & Hubbard 2024), regardless of the method used.

For Uranus, the available data are limited in comparison to the gas giant planets. However, several studies have investigated Uranus’ internal structure in detail. For example, Nettelmann et al. (2013) presented physical interior models of Uranus (and Neptune) where the interior was divided into three layers: a rocky centre surrounded by two adiabatic and homogeneous envelopes of hydrogen, helium, and water. They found the outer envelope metallicity to be Z ≤ 8% for modified shape and rotation data and Z ≤ 17% for original values provided by Voyager 2. In each of these cases, solar metallicity models (Z = 0.015) were still found to be possible. Neuenschwander & Helled (2022) presented empirical interior models of Uranus and Neptune. It was shown that faster rotation and/or deep winds favour centrally concentrated density distributions. It was also argued that a more accurate determination of the gravity field can significantly reduce the possible range of internal structures. In addition, Neuenschwander & Helled (2022) predicted the values of J₆ and J8 for both planets based on their model. Since information about composition and temperature is crucial to get a complete understanding of the planet, several algorithms have been proposed to interpret empirical density profiles in terms of temperature and composition (for example Podolak et al. 2022; Neuenschwander et al. 2024).

Podolak et al. (2022) presented a method that combines a top-down with a down-up approach. For the top-down part, the algorithm starts at the planetary surface using the expected temperature and composition values of Uranus at the 1-bar pressure level. It then divides the planets into shells and moves inwards. At each consecutive shell, the algorithm attempts to minimise the mass fraction of high-Z material whenever possible (first water, then rock, then iron) while keeping the temperature constant. If no mixture matches the corresponding density and pressure, the temperature is increased while keeping the same composition. This procedure results in a composition that favours the material with the lowest mean molecular weight and keeps the temperature as low as possible. For the bottom-up part, the algorithm starts with a postulated temperature and composition of the centre. Continuing outwards, the algorithm checks whether the same composition matches the density and pressure while keeping the temperature monotonically decreasing. If the check fails, a composition is inferred that retains the previous temperature by decreasing the metal content (first iron, then rock, then water). The top-down and bottom-up algorithms were combined by comparing the inferred composition profiles, and choosing the solution which minimises the amount of water (layer by layer). Using this method, Podolak et al. (2022) concluded that the water-to-rock ratio in Uranus cannot be much smaller than 0.5.

The algorithm presented in Neuenschwander et al. (2024) is based on a top-down approach. It begins with the planetary surface values that include Uranus’ atmospheric composition and temperature. Then, the algorithm attempts to move inwards, assuming a constant composition and entropy. If this fails in some layer, the algorithm uses the Ledoux criterion (Ledoux 1947, also see Section 2.2) by assuming the fluid to be precisely at the Ledoux limit: the temperature gradient is compensated exactly by the composition gradient to yield neutral buoyancy. The algorithm therefore finds a unique temperature-composition profile for a given density profile. Neuenschwander et al. (2024) concluded that their empirical density profiles can be interpreted as having non-adiabatic interiors with composition gradients. Compared to fully adiabatic models, these solutions have higher central temperatures. The inferred water-to-rock ratio was found to be between 2.6 and 21.

In this paper, we present a new random algorithm to determine the composition-temperature profiles of improved empirical structure models. Introducing randomness yields a broader range of solutions which is more in line with the spirit of using empirical models. We improve our empirical structure models by calculating the Theory of Figures (Zharkov & Trubitsyn (1978)) up to tenth order. Our Theory of Figures (ToF) implementation includes the possibility of incorporating barotropic differential rotation. We also employ better approximations than previous ToF implementations which lead to a significantly enhanced accuracy. Finally, we include a treatment of the atmosphere which leads to more realistic outermost parts of Uranus in our empirical structure models.

Our paper is structured as follows: in Section 2 we explain the newly developed methods, such as the inclusion of the atmosphere, the development of a tenth-order ToF, and the new random algorithm that is used to infer the composition, temperature, and stability against convection throughout the planet. We present our results in Section 3. A comparison of our results with previous studies and a discussion on the limitations of this work are presented in Section 4. A summary and conclusions are presented in Section 5.

2 Methods

Our empirical models are based on polytropes which relate the density ρ with the pressure P inside the planet via: $$\rho (P)={\left(\frac{P}{K}\right)}^{\frac{n}{n+1}},$$(1)

where K and n are constants. Following Neuenschwander & Helled (2022), we begin with three consecutive polytropes defined by constants K_1,2,3 and n_1,2,3, where the polytropes are separated by two transition radii. Additionally, we added a fourth polytrope representing the outermost region of the planet which is consistent with atmospheric models of Uranus. Specifically, in the outermost region (atmosphere) we impose the following polytropic relation: $$\rho \left(r\ge r_{\ast }\right)={\left(\frac{P\left(r\ge r_{\ast }\right)}{K_{\text{atmos}}}\right)}^{\frac{n_{\text{atmos}}}{n_{\text{atmos}}+1}},$$(2)

where K_atmos and n_atmos are fitted to the atmosphere model taken from Hueso et al. (2020). P is the atmospheric pressure profile obtained by numerically integrating the condition for hydrostatic equilibrium, ρ(r)⁻¹ ∇P(r) = ∇ U(r), where U is the total potential (gravity and rotation). We let r∗ to be freely chosen by the algorithm as long as P(r_*) ≥ P_min and set P_min = 100 bar. In order to ensure mass conservation we applied: $$\rho (r)\to \frac{\rho (r)}{4\pi {\displaystyle {\int }_0^{r_{\max }}{r}^2}\rho (r)\text{d}r}\cdot M,$$(3)

after using Eq. (2) until convergence was accomplished.

The polytropic relation between the density ρ and the pressure P in Eq. (1) is only inspired by, not based on physical materials. The constants K_1,2,3 and n_1,2,3 are chosen to exclusively fit the gravity data regardless of any material properties, and therefore the density profiles are empirical. As we discuss in Section 2.2, the inferred P - ρ relations are then interpreted in terms of composition and temperature using physical Equations of State (EoS).

We employed the Markov-chain-Monte-Carlo (MCMC) emcee algorithm (Foreman-Mackey et al. (2013)) and used the ToF (Zharkov & Trubitsyn (1978)) with N = 4096 spheroids and no interpolation to find empirical profiles of Uranus. They assume uniform rotation and agree within 1σ with the measured gravitational moments listed in Table 1.

Figure 1 presents our 469 4-polytrope models for Uranus. The density ρ is given as a function of the average normalised radius calculated by the ToF (see discussion below). We also show the pressure as a function of density P(ρ) in the outermost part of the planet which fits an atmosphere model of Uranus (Hueso et al. 2020). The atmosphere model is shown by the white dots and is followed by our models until P(l*) ≥ P_min = 10⁷ Pa.

Overall, our model consists of four polytropes separated by three transition radii where only three polytropes are chosen freely to fit the gravitational moments J₂ and J₄ . The measured gravitational moments include the non-hydrostatic contributions from Uranus’ atmosphere dynamics (e.g Kaspi et al. 2013; Soyuer et al. 2020). Our models assume uniform rotation and correspond to a fully hydrostatic planet without differential rotation or winds.

Fig. 1 Uranus’ density ρ as a function of the average normalised radius r/R. Other density profiles that have been published by previous studies are also presented for comparison (Helled et al. 2011; Nettelmann et al. 2013; Vazan & Helled 2020). The panel on the right show the pressure as a function of density P(ρ) in the outermost region. The white dots show the atmosphere model by Hueso et al. (2020). The full P(ρ) distribution is presented in appendix C in Figure C.1.

2.1 Tenth order Theory of Figures

The models presented in Figure 1 are based on the Theory of Figures (ToF) by Zharkov & Trubitsyn (1978). The ToF inserts a spheroid Ansatz (surfaces parameterised by a function r_l(ϑ)) into the total potential U. ϑ and l are the polar angle and volumetric mean radius of the spheroid, respectively. U = V + Q is the sum of the Newtonian potential V = −G ʃ d³r′ρ/|r′ − r| and the centrifugal potential $Q=-\frac{1}{2}{\omega }^2{r}^2{\sin }^2(\vartheta )$. The details of inserting r_l(ϑ) into U = U(r) and expanding the resulting expression order by order are discussed in detail in Appendix A. Originally, ToF was developed up to third order (Zharkov & Trubitsyn 1978) and later was expanded to fourth (Nettelmann 2017) and seventh order (Nettelmann et al. 2021). For this work, we developed a tenth order ToF where the necessary expressions denoted as $f_n,f_n^{\prime }$, and A_n, which are themselves functions of so- called figure functions s_n , reach tenth order. Further details on our tenth order ToF implementation can be found in Appendix A. In Subsection A.1, we discuss potentials Q that correspond to barotropic differential rotation and demonstrate that our ToF implementation accurately incorporates barotropic differential rotation.

ToF calculates and converges on the figure and the gravitational moments J_2n of the planet with a density profile ρ(r). The gravitational moments appear in the expansion of the external part of the Newtonian potential into spherical harmonics: $$V_{\text{ext}}(r,\mu )=\frac{GM}{r}\left(1-{\displaystyle \sum _{n=1}^{\infty }{\left(\frac{R_{\text{eq}}}{r}\right)}^{2n}}{J}_{2n}{P}_{2n}(\mu )\right),$$(4)

where R_eq is the equatorial radius and µ = cos(ϑ).

The ToF assumes no contributions of odd J_2n+1 and no dependence on the azimuthal angle φ. In order to calculate the J_2n , the ToF divides the planet into a given number of spheroids N . A larger N leads to more accurate solutions for the figure and the planetary gravitational moments.

In order to test the accuracy of our numerical ToF implementation, we use data provided by Wisdom & Hubbard (2016). As a benchmark, they use a simple Jupiter model with a single polytropic Equation of State of p = Kρ² and assume a rotation period of 9^h55^m29.7^s as well as an equatorial radius of R_eq = 71 492 km. Thanks to the simplicity of this model, it is possible to independently calculate the gravitational moments J2_n using spherical Bessel functions. We use these Bessel values of J_2n as a reference and plot the relative difference to our solutions in Figure 2 for J₂ (top) and J₈ (bottom) as examples. We present different orders of the ToF as well as the implementation used for Movshovitz & Fortney (2022). We also investigate the influence of approximating relevant expressions with spline interpolation. The parameter n_x stands for the number of points that were calculated without interpolation. n_x = N hence means that no interpolation was used, while n_x = N/32 means that only every 32nd point was actually calculated in the corresponding arrays, with all the intermediate points being estimated by using a cubic spline interpolation.

Several conclusions can be made from this comparison. First, our implementation of the ToF is more accurate than the one presented in Movshovitz & Fortney (2022) and can reach relative precisions of the order of 10⁻⁷ instead of 10⁻³ –10⁻⁵ for the gravitational moments up to J₈. Second, using seventh and tenth order ToF yield better accuracies for all the gravitational moments when using a sufficient high number of layers. Tenth order ToF is superior to seventh order ToF for gravitational moments beyond and including J₆ . Third, the spline interpolation only influences accuracy when N is smaller than ∼2¹².

The implementation of the ToF presented here is more accurate than the one in Movshovitz & Fortney (2022) and the two implementations have differences. The first difference between the two concerns a small typo in equation (B.17) for f₄^′ in Nettelmann (2017). The correct equation is¹: $$\begin{array}{ll}{f}_4^{\prime }=\hfill & \frac{1}{3}{s}_4-\frac{9}{35}{s}_2^2-\frac{20}{77}{s}_2{s}_4-\frac{45}{143}{s}_2{s}_6-\frac{81}{1001}{s}_4^2\hfill \\ \hfill & +\frac{72}{385}{s}_2^3+\frac{4579}{5005}{s}_2^2{s}_4-\frac{12798}{25025}{s}_2^4.\hfill \end{array}$$(5)

Since Movshovitz & Fortney (2022) relies on Nettelmann (2017) instead of Eq. (5), they introduce a small error into their calculations for fourth order ToF. However, the resulting deviation is very marginal and only affects fourth order. This does not cause the differences visible in Figure 2. We believe that our improved accuracy is achieved by using the trapezoid rule to calculate the average density via M = 4π ƒ ρ(r)r²dr, instead of approximating $M\approx \sum _i\text{Δ}{\rho }_{i}4\pi r_i^3/3^{}$².

Figure 3 shows the runtime (in seconds) for the calculations depicted in Figure 2. The calculations were performed on a higher-end conventional laptop in Python. The ToF calculates all gravitational moments J_2n simultaneously, so the presented runtime does not depend on a specific gravitational moment. We find that for fourth order ToF, the implementation of Movshovitz & Fortney (2022) is faster compared to our code. A spline interpolation does not lead to significant time savings. We also find that for higher orders, our implementation is faster. This has been achieved by directly computing the ToF coefficients in a machine-readable format and hence avoiding external coefficient files that need to be loaded repeatedly. Furthermore, we find that a spline interpolation can save significant computing time.

Fig. 2 Relative difference of the Bessel gravitational moment solutions J2 (top) and J8 (bottom) compared to the ToF solutions. Results are shown for different numbers N of used spheroids. We compare our implementation (This work) to the one used in Movshovitz & Fortney (2022) (M&F 22) for different orders of the ToF. We also show the influence of using cubic spline interpolation with the parameter nx that stands for the number of points that were calculated without interpolation.

Fig. 3 Runtime in seconds of the calculations shown in Figure 2 as a function of the number N of used spheroids in the ToF. We refer to the legend and caption of Figure 2 for further explanations.

2.2 A new algorithm to infer the temperature profile and composition

In this subsection we describe the new algorithm to infer the composition and temperature T profiles for the empirical models. We improve over previous work (for example Podolak et al. 2022; Neuenschwander et al. 2024) by introducing randomness into the process. This allows us to find a significantly wider range of possible solutions for each density profile, which is more consistent with the spirit of empirical models to find ‘outside the box’ solutions. The main steps of the algorithm can be summarised as follows:

• We begin with an empirical structure model that provides the planetary density p = [ρ₀,…,ρ_N] and pressure P = [P₀,…, P_N]. Since our EoS for Hydrogen and Helium (H- He) are limited to temperatures above 100 K (Chabrier et al. 2019), the first few entries of ρ and P are skipped.

• The algorithm tries to find a starting tuple (ρ₀, P₀, T₀, X₀, Y₀, Z₀) that is consistent with the EoS and the empirical profile and chooses randomly if multiple exist.

• The algorithm continues to the next spheroid (henceforth called layers) with a stepsize Δn, which we took to be 50 or 100 for our N = 4096 models. Smaller step sizes only lead to larger computing times with no significant change to our results.

• The algorithm tries to find isentropic solutions first. With this choice, we bias our algorithm towards finding isentropic solutions (corresponding to convective regions).

• If there is no isentropic solution, we try changing the composition. Priority is given to changing the composition of the materials already present. If that fails, heavier components are added. In both cases, we check the Ledoux criterion introduced below and choose randomly if multiple options exist.

• This is repeated until the planetary centre is reached or it is impossible to generate new candidates (p_n, P_n, T_n, X_n, Y_n, Z_n) that fulfil all necessary conditions.

The details of our algorithm are summarised in the flow chart presented Figure 4. We allow the skipping and ignoring of up to three layers to ensure that changing the composition is truly necessary and not just a fluke due to the uncertainties of the EoS and the empirical model. Due to the uncertainties of the EoS, we set |S_n − S_n−∆n|/|S_n−∆n| ≤ 10⁻² for the isentropic condition.

Our implementation of the random algorithm in Figure 4 supports up to four materials at the same time, using the assumption of ideal mixing for the total density ρ. $$\frac{1}{\rho }=\frac{X+Y}{{\rho }_{X+Y}}+\frac{Z^0}{{\rho }_{Z^0}}+\frac{Z^1}{{\rho }_{Z^1}}+\frac{Z^2}{{\rho }_{Z^2}}.$$(6)

The ratio of the Hydrogen and Helium (H-He) mass fractions X and Y are fixed to be proto-solar. Therefore, they are treated as a single component. As a result, our algorithm can treat up to four materials simultaneously, distributed into three (composite) components. We use the Chabrier et al. (2019) EoSs for the H- He mixture. For the heavier materials Z we assume water (H₂O), rocks (SiO2), and iron (Fe) using the EoSs of More et al. (1988).

Often, the interiors of the giant planets in the Solar System are assumed to be isentropic due to large-scale convection. However, non-isentropic regions can exist and could be a result of boundary layers or composition gradients. Both need to be locally stable, which is tested by displacing a bubble of material. If the bubble oscillates around its original position (with a frequency N ∈ ℝ), the surrounding is locally stable. If the frequency N ∈ iℝ is imaginary, the bubble gets pushed further away from its original position. The surrounding is interpreted to be unstable, that is, subjected to convection. To check for local stability in our algorithm, we use the slightly modified Ledoux criterion from Paxton et al. (2013). Following Kippenhahn et al. (2013), the Brunt-Väisälä frequency N of a displaced bubble of material is (see also Unno et al. 1989, for details): $$N^2=\frac{g^2\rho }{P}\frac{{\chi }_T}{{\chi }_{\rho }}\left[{\nabla }_{ad}-{\nabla }_T+\underset{=:B}{\underbrace{\frac{{\chi }_{\rho }}{{\chi }_T}{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\mu }\right)}_{P,T}\frac{ \text{d\hspace{0.17em}ln\hspace{0.17em}}\mu }{ \text{d\hspace{0.17em}ln\hspace{0.17em}}\rho }\frac{ \text{d\hspace{0.17em}ln\hspace{0.17em}}\rho }{ \text{d\hspace{0.17em}ln\hspace{0.17em}}P}}}\right].$$(7)

Following similar steps as in Paxton et al. (2013), B can be calculated numerically as: $$B=\frac{{\chi }_{\rho }}{{\chi }_T}\frac{\ln \rho \left(P_n,T_n,A_{n+\text{Δ}n},\right)-\ln \rho \left(P_n,T_n,A_n\right)}{\ln P_{n+\text{Δ}n}-\ln P_n},$$(8)

where A_n := (X_n, Y_n, Z_n). We provide a detailed derivation in Appendix B. Requiring N² ≥ 0 leads to the Ledoux criterion stated in Figure 4: $${\nabla }_T\le \beta \left({\nabla }_{ad}+B\right),$$(9)

where $${\nabla }_{ad}:={\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}T}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}\right)}_S,{\nabla }_T:=\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}P}.$$(10)

The β pre-factor ∈ [0,1] assists the algorithm to find lower central temperature solutions more quickly. Excluded solutions with β(∇_ad + B) ≤ ∇_T ≤ ∇_ad + B would have higher temperature gradients. We routinely used β = 0.7 and still found plenty of high central temperature solutions (see Section 3).

The β pre-factor parameterises the uncertainty of the energy transport in inhomogeneous layers. If the composition gradient is stable there are two possibilities for the energy transport. The first is double-diffusive convection and the second is conduction. For planetary interior conditions, double-diffusive convection is still poorly understood. It is unclear whether layers can form and survive on long timescales (Wood et al. 2013; Fuentes et al. 2022). In addition, it is unclear what the appropriate temperature profile is when composition gradients exist (for example Garaud 2018). For thermal conduction, the temperature gradient is well-defined in principle. However, it depends on an unknown local luminosity in the deep interior and the thermal conductivity, which is highly uncertain for the relevant conditions. Given these issues, we remain agnostic about the temperature gradient and let the algorithm randomly pick from a temperature consistent with the empirical model and Ledoux stability. Our choice of β = 0.7 is motivated by the fact that it is unlikely that the temperature gradient is very similar to the Ledoux gradient. Such regions would only be marginally stable and could also lead to extreme temperatures.

Figure 5 presents the two different interior models we consider in this study to infer Uranus’ composition. The top panel corresponds to a model of Uranus admitting up to three materials per layer (U_3,comp). We begin with H-He and H₂O in Uranus’ atmosphere. As soon as the algorithm fails to find more solutions with these components, we allow the addition of SiO₂ and remove H-He. When no solutions can be found with H₂O and SiO₂ , then H₂O is removed and we allow for the addition of Fe. This model is more traditional in the sense that it follows a structure in which lighter elements do not exist in the deep interior.

The bottom panel on the other hand shows a model of Uranus admitting up to four materials per layer (U_4,comp). Just as before, we begin with H-He and H₂O in the outermost part of the planet. As soon as the algorithm fails to find more solutions with these components, we allow the addition of SiO₂, but still allow for the presence of H-He. If at some point no solutions with H-He, H₂O, and SiO₂ can be found, then H-He are removed and we allow for the addition of Fe. This configuration allows solutions with ‘fuzzy cores’ as we allow H-He to exist in Uranus’ deep interior (Helled & Stevenson 2017; Valletta & Helled 2022).

Fig. 4 Our algorithm for inferring the temperature profile and composition of the empirical density (and pressure) profiles.

Fig. 5 Schematic drawing of the two possible composition models considered in this work. Top: U3,comp model that allows the presence of up to three elements simultaneously in a given layer. Bottom: U4,comp model that allows the presence of up to four elements simultaneously in a given layer. The ratios between different components are chosen freely and composition gradients are allowed but not necessary. The transition radii (white gaps) vary and are placed by the algorithm.

3 Results

In total, the random algorithm derived composition-temperature profiles for 469 empirical density profiles (shown in Figure 1), which yielded 14659 U_4,comp and 13645 U_3,comp solutions. First, we present our predictions for the gravitational moments beyond J₆ up to J₁₄ using our updated empirical models (Figure 1) in Subsection 3.1. In Subsection 3.2 we present the inferred composition and temperature profiles using the new algorithm for three examples of our U_4,comp inferred models. We compare our complete results for the two cases U_3,comp and U_4,comp in Subsection 3.3 and discuss relevant correlations between parameters of interest in Subsection 3.4.

3.1 Gravitational moments predictions

Our tenth order ToF allows us to calculate Uranus’ gravitational moments accurately and to high order. We first use our models to predict Uranus’ gravitational coefficients J₆ up to J₁₄ presented in Table 2. We find that J₆ = (5.3078 ± 0.3312) • 10⁻⁷ and J₈ = (−1.1114 ± 0.1391) • 10⁻⁸. These values are in agreement with the values predicted by Neuenschwander & Helled (2022). For the higher order gravitational moments, we find J₁₀ = (2.8616 ± 0.5466) • 10⁻¹⁰, J₁₂ = (−8.4684 ± 2.0889) • 10⁻¹² and J₁₄ = (2.7508 ± 0.7944) • 10⁻¹³ . These values could be measured via a future Uranus mission and can constrain Uranus’ internal structure (for example Helled & Fortney 2020) and atmosphere dynamics (for example Kaspi et al. 2017; Soyuer et al. 2020).

3.2 Composition and temperature profiles

In this subsection we present the composition and temperature profiles and identify the stable regions against convection using the random algorithm. First, we show three different empirical density profiles (U_4,comp model) which are presented in Figure 6. These three models were chosen as examples to illustrate the diversity of possible solutions. The dotted lines correspond to individual solutions. The average of all individual solutions of a given profile is shown by a thick line and the 1σ range is shown by the coloured region.

The low central density profile presented in Figure 6a leads to solutions with a well-mixed and rock-rich inner region. We find that there are two large convection zones (with a corresponding Brunt-Väisälä frequency of N² < 0) separated by a broad stable layer at P ~ 20 GPa and ~0.75r_U, where r_U is the average surface radius of Uranus. H-He is present throughout the planet, even in the centre with small proportions. The inferred central temperatures are high, of the order of 50 000 K.

On the other hand, the model presented in Figure 6b which also has a low central density and a somewhat similar density distribution, is found to have a rather different solution. This model is incompatible with a well-mixed inner region and is found to have only one large convection zone in the outer region of the planet that stops at P ~ 20 GPa and ~0.75 r_U. The deep interior consists of composition gradients including H-He with typical central temperatures of ∼20 000 K. Although at first sight the two density profiles seem rather similar they are sufficiently different to result in significantly different solutions for the inferred temperature and composition profiles. This clearly illustrates the problem of degeneracy in interior modelling. A direct comparison of the two density profiles is presented in Appendix C in Figure C.2.

Finally, the high central density profile shown in Figure 6c is characterised by three distinct zones. First, there are two zones that are clearly convective and are separated by a sharp stable layer at ~0.8 r_U and P ~ 10 GPa. The water-to-rock ratio of the second middle zone varies significantly for different individual solutions. The final innermost zone is found to be iron dominated in composition, with only negligible amounts of rocks and water. H-He is present in the middle zone, but not in the planetary centre. The central temperatures are high and are found to be ∼40 000 K on average.

All of these profiles possess a convective ionic water region (according to the p-T line in Redmer et al. (2011)), consistent with the requirement to generate Uranus’ magnetic field. The three following criteria are fulfilled at the same time: N² < 0 holds, H₂O is present in the composition mix and the red temperature lines cross the ionic H₂O region. This could explain the origin of Uranus’ magnetic field that can be attributed to the high electrical conductivity profile in this region (for example Soyuer et al. 2020, and references therein). The majority (77% in the U_3,comp case and 83% in the U_4,comp case) of our profiles were inferred to have a convective ionic water region.

The inferred magnitudes of the Brunt-Väisälä frequencies in convective regions are unrealistically large. This is likely because the polytropic models can have density gradients different than expected from an adiabatic and homogeneous region, which would be too unstable (see also Neuenschwander et al. 2024). Here we only examine whether N² is negative or positive to interpret our profiles regarding convective stability of the regions with composition gradients. In Section 4.1 we discuss the temperature-sensitivity of the inferred composition. Future work should address this topic and ensure to provide realistic density gradients and Brunt-Väisälä frequencies, which could provide further constraints when modelling Uranus’ interior.

3.3 U_3,comp versus U_4,comp

Below, we present the results of all the empirical models combined. We begin with comparing these two cases focusing on the inferred composition and temperature profiles.

Uranus’ bulk composition

The assumed internal structure significantly affects the inferred bulk composition. Figure 7 shows the complete range of inferred compositions and central temperatures for cases U_3,comp (top) and U_4,comp (bottom). We find that most of the solutions for the U_3,comp model are water-rich and that rock-dominated interiors are rare. The few rock-dominated solutions we found require extreme central temperatures exceeding ∼50 000 K that are rather unrealistic. Figure 7 clearly shows for the U_3,comp case that the water mass fraction is inversely correlated with the central temperature, where lower temperatures correspond to more water-rich interiors. More specifically, the water mean mass fraction values possess a Spearman’s rank correlation coefficient of −0.57 at a p-value of <0.005 when compared to the mean central temperatures.

Unlike the U_3,comp case, the U_4,comp models allow H-He to exist in the deep interior. This significantly affects the inferred rock mass fraction and rock-dominated solutions are found to be much more common. Adding H-He to the mixture of materials compensates for higher abundances of heavier materials (such as SiO₂). This in turn influences the central temperatures, which are now always hotter than ∼15 000 K (compared to ≳8000 K for the U_3,comp models). There is also no correlation between the water mass fraction and the central temperature. The water mean mass fraction values now possess a Spearman’s rank correlation coefficient of −0.03 at a p-value of 0.55 when compared to the mean central temperatures. Furthermore, we find that typically water is less common than rock for all central temperatures.

We next summarise all composition results in Figure 8. It shows the overall mass fraction of the different materials for both scenarios U_3,comp (top) and U_4,comp (bottom). The dots displayed in Figure 7 are weighted averages to account for the fact that we did not necessarily find an equal amount of composition-temperature profiles for each empirical density profile. Figure 8 presents the mean and 1σ range of all these weighted averages. We find that for the U_3,comp models water is generally the most abundant material, followed by rocks, H-He, and iron. However, from the U_4,comp models, it is clear that if H-He can exist in the deeper interior we generally expect a water-to-rock ratio that tends to be significantly smaller than one (~0.6).

Fig. 6 Three example empirical density profiles of Uranus (depicted from left to right). Their density, composition, and temperature profile inferred by the new random algorithm is depicted from top to bottom, respectively. The extended phase diagram of water based on Redmer et al. (2011) is shown in addition to the temperature profiles. The stability against convection is presented using the Brunt-Väisälä frequency N2 in the bottom last plots.

Fig. 7 Inferred composition and temperature for our Uranus empirical models for U3,comp (top) and U4,comp (bottom). The x-axis corresponds to the mass fraction of a material and the y-axis shows the central temperature. The colours represent different materials (see legend of Figure 6) and the error bars indicate the 1σ range of solutions found by the random algorithm for each density profile. To improve readability, error bars for the central temperature Tc are not shown. They are included in Figure 11.

Convectivity of Uranus’ interior

Figure 9 summarises the convective regions inferred by our algorithm for the two cases presented in Figure 5. We begin by considering the histogram for the U_3,comp case (top) and note that we can distinguish between two broad categories: Roughly half of our solutions are mostly non-convective (~20% of the planet w.r.t. its radius is convective) while the other half is mostly convective (~90% of the planet is convective).

In the U_4,comp case (bottom), we can identify the same two categories of solutions in the histogram. However, there are more convective solutions overall. This is not surprising given the fact that there is more freedom to mix different elements, that is, more options to choose from for the U_4,comp case. Some of these additional options are isentropic and compatible with convection, and since our algorithm is biased towards such solutions, we tend to find more convective solutions in the U_4,comp case. The peaks of solutions that are convective for only ~20% of the planet correspond to profiles possessing only a convective outermost part. The profile depicted in Figure 6b serves as an excellent example.

Considering now the other plots, we see good agreement with the above discussion: almost all profiles posses a convective outermost part. And almost all profiles possess a transition region around ~0.75r_U. This seems to be a robust prediction given our 4-polytrope empirical models. The profiles then start to disagree below ~0.75r_U, splitting up into the two aforementioned categories in roughly equally sized groups.

Fig. 8 Uranus’ inferred composition using the new random algorithm. We depict results for the U3,comp (top) and U4,comp (bottom) models. The pie charts summarise the composition histograms from Figure 7, and the numbers are the means and 1 σ ranges of all empirical models shown as dots in Figure 7.

3.4 Correlations

The random algorithm yielded many solutions for each empirical density profile, enabling the statistical analysis of the solutions. Here, we show the correlations of a few selected planetary properties for the U_4,comp model: The mass fractions of different materials, central temperature and density, gravitational moments, and the mean location of the ionic H₂O regions that are convective. Unless specifically stated, the following correlations are also valid for the U_3,comp model. The complete data for the U_3,comp and U_4,comp models are provided for a qualitative analysis in Appendix C (Figures C.3 and C.4). The stated p-values, which are all below 0.005, stay below 0.005 even when combining them to account for the Bonferroni correction³.

Fig. 9 First and third plot: histograms showing the percentage of Uranus’ radius that is convective versus the number of solutions. In the U4,comp case, more than 100 empirical profiles were inferred to be convective for ~90% of Uranus’ radius. Second and last plot: radial distributions of the convectivity for all solutions. In both cases, only ~20% of the profiles were found to be convective at ~0.75rU.

Uranus’ magnetic field generation

Uranus has been observed to possess a magnetic field (Ness et al. 1989). Hence, interior modelling needs to explain this observation. According to Maxwell’s equations, a magnetic field can be generated by moving charged particles. Given the right temperature and pressure conditions, water can become ionic (charged). Most of our models (77% in the U_3,comp case and 83% in the U_4,comp case) possess a convective region with ionic water. The convection would determine the velocity of the charged water particles generating the magnetic field. We now want to further investigate the properties of these regions. First, we explore whether a higher abundance of a given element influences the location of the convective ionic water region. Figure 10 shows the mean location of the convective ionic water region as a function of the total H-He mass fraction (M_X+Y/M).

We find that a higher H-He mass fraction correlates with a deeper location of the convective ionic H₂O region. This can be explained by the larger extent of the H-He dominated outermost part of the planet, which pushes the H₂O rich regions further into Uranus’ interior. Typically, the water mass fraction in our Uranus models is ~0.2 in the convective ionic H₂O regions. Some U_4,comp models consist of a water mass fraction of ~0.5 and nearly 1 for U_3,comp models in these convective regions. The full data for the water mass fraction in convective ionic H₂O regions is presented in Appendix C. At the moment, the minimum ionic water mass fraction that is required to ensure a high enough electrical conductivity to generate Uranus’ magnetic field remains unknown (for example Soyuer et al. 2020; Soderlund & Stanley 2020). Our work suggests that the inferred water mass fraction of ~0.2 is sufficient to explain Uranus’ dynamo. Clearly, this topic should be investigated in future research in detail.

It should be noted that not all of our U_4,comp solutions possess a convective ionic H₂O region in the first place. We therefore only show the 83% of solutions that possess one in Figure 10. For the extension of the convective ionic H₂O region, we found an upper limit of ~0.15r_U (see Appendix C).

Fig. 10 H-He mass in Uranus versus mean location of the convective ionic H2O region. The error bars indicate the 1σ range of solutions found by the algorithm for density profiles possessing an ionic-water region. Profiles without a convective ionic H2O region are not presented.

H-He in the deep interior of Uranus

Here we investigate the correlation between the centre H-He mass fraction M_X+Y,c/M_c as a function of the central temperature T_c. The results are shown in Figure 11.

We find that the central temperature is inversely correlated with H-He abundance for lower central density empirical U_4,comp models. Very high central temperatures are incompatible with a large H-He centre mass fraction. Overall, we find an upper limit of ~0.1 for the H-He mass fraction in Uranus’ centre. Higher central density models do not include any H-He in the centre, since including H-He and Fe simultaneously as H₂O and SiO₂ is not allowed by construction (see Figure 5). U_3,comp solutions could only find H-He in the central region if the central density is sufficiently low so it can be reproduced using a pure H-He-H₂O mixture.

Fig. 11 H-He mass fraction in Uranus’ centre versus the central temperature Tc. The error bars indicate the 1σ range of solutions found by our algorithm.

The influence of future missions on Uranus’ interior

Figure 12 shows the total H-He mass fraction M_X+Y/M as a function of the higher order gravitational moments J₆ and J₈. We find clear correlations for both cases. We can therefore conclude that an accurate measurement of J₆ and J₈ can constrain M_X+Y/M in Uranus: J₆ and J₈ are most sensitive in the outermost region of the planet, where H-He is most abundant. Since the mean location of the convective ionic H₂O region is also correlated with M_X+Y/M (see Figure 10), accurate measurements of J₆ and J₈ can further constrain the location of the convective ionic water region. The fact that the spearman rank-order correlation coefficients are almost opposite to each other is no coincidence since the gravitational moments themselves are strongly correlated.

It is worth noting that a future mission will also more tightly constrain J₂ and J₄. Our work did not find significant correlations between the above quantities of interest and the gravitational moments J₂ and J₄.

4 Discussion

4.1 Temperature sensitivity of the composition

In Section 3.2, we presented composition-temperature profiles and discussed that the homogeneous regions may exhibit a negative Brunt-Väisälä frequency too far away from zero. Another way to interpret this is that the inferred temperature profile is not adiabatic even though the region is homogeneous in composition and the algorithm labels it as a convective region. This occurs because empirical models can sometimes have a density gradient that is flatter than expected for an isentrope at a constant composition. This is true for empirical models in general and not just the ones presented here. This behaviour is noticeable when a physical EoS is used to interpret the composition and temperature of these models. Otherwise, it may remain hidden.

In a region with a small density gradient, the two ways to match the low density of a deeper layer are to either increase the entropy (and leave the isentrope) or decrease the mean molecular weight of the mixture. Since the latter is not allowed in our algorithm, this can cause the entropy to increase inwards. As noted previously, this is not a stable configuration, because convection should quickly drive the temperature gradient towards an adiabatic one.

In order to investigate how this issue affects our solutions, we forced an adiabatic temperature on the pressure-density profiles in regions labelled as convective by the algorithm. Forcing an adiabatic temperature gradient in convective regions leads to colder temperatures, the planet could have more lighter materials. While this approach is not fully self-consistent, it allows us to estimate the expected change in the inferred composition. We removed the constraint of an inwards monotonically increasing mean molecular weight to ensure that solutions are found despite the flat density gradients, for example. We then compared these new composition profiles with our previous results.

Figure 13 shows the differences in composition for U_3,comp (top) and U_4,comp (bottom) models.

For the U_3,comp models, the differences between the models are similar to Normal distributions peaked around zero. Therefore, despite the non-adiabatic temperature gradients in convective regions, we can conclude that our findings for the U_3,comp are rather robust. On the other hand, for U_4,comp models, the inferred composition differs more and follows a clear trend. While the changes in H-He and iron are small, the new inferred composition is biased towards more water and less rock. The magnitude of the shift is similar to the uncertainties shown in Figure 8.

The changes in composition directly depend on the temperature. Figure 14 shows the temperature-density relation for a protosolar H-He mixture, pure water, rock, and iron. We used our EoS to calculate the density as a function of pressure at three different constant pressures of 100, 500, and 1000 GPa. These pressures were chosen to represent the region inside the planet where most of the mass is.

As can be seen from the figure, at lower temperatures between about 1000 to 10 000 K, the density does not strongly depend on the temperature, particularly for the more refractory rocks and iron. However, for temperatures beyond ~ 10 000 K the density strongly decreases with temperature. As discussed in Section 3.2, our models tend to be very hot, because of the composition gradients and the deviations from the isentropes in the convective layers. The combination of the hot interiors and the temperature-density relation causes the inferred composition to be sensitive to the temperature.

Our test reveals that the inferred composition for the U_3,comp should be robust, while this is less true for the U_4,comp models, where we find a significant change of the water and rock mass fractions. To improve future constraints on the composition and temperature of Uranus, empirical models could restrict the poly-tropic indices so that density gradients do not become too flat, for example. We hope to address this in future studies.

Fig. 12 Total H-He mass fraction MX+Y/M versus the higher order gravitational moments J6 (top) and J8 (bottom). The errorbars indicate the 1σ range of solutions found by the random algorithm for each density profile.

Fig. 13 Absolute change in composition (H-He, water, rock, and iron) when enforcing a perfect adiabatic temperature profile. Positive numbers indicate that the forced adiabatic temperature profiles possess more of a given component. The results are shown for both U3,comp (top) and U4,comp (bottom) models.

Fig. 14 Density as a function of the temperature for three isobars at 100, 500, and 1000 GPa, with the pressure increasing from the bottom to the top of the figures. The four panels show the density-temperature relation of a protosolar H-He mixture, pure water, rock, and iron.

4.2 Comparison with previous work: methods

Our random algorithm differs in various ways compared to previous studies which focused on interpreting empirical models in terms of temperature and composition. For example, the algorithm of Neuenschwander et al. (2024) starts in the outermost region of the planet, where the composition and temperature are reasonably well known. Then it moves towards the planetary centre layer by layer searching for isentropic, constant composition solution first. The algorithms are rather different when no isentropic and constant composition solutions are found: Neuenschwander et al. (2024) numerically solve Eq. (7) for N² = 0 which leads to a unique solution. This approach has two limitations: First, it can lead to very large temperature gradients if the composition changes rapidly. Second, it is unclear whether these solutions should be interpreted as marginally stable against large-scale convection or are essentially already unstable. On the other hand, our algorithm proposes numerous candidates that satisfies the inequality N² ≥ 0 and chooses randomly amongst them. The results of this study are therefore more agnostic in terms of the temperature gradient in non-convective layers. As a result, we find a large range of possible temperature and composition profiles for a single empirical density profile.

The algorithm proposed by Podolak et al. (2022) is rather different from the one we use here. As we discuss in detail in Section 1 it combines top-down and bottom-up approaches and minimises the amount of water layer by layer. Furthermore, their algorithm focuses on finding constant temperature rather than constant entropy solutions for both approaches. This leads to mostly isothermal solutions compared to our mostly isentropic solutions. As a result, they do not have information on entropy and can not identify the stability against convection in their models.

4.3 Comparison with previous work: results

Despite the above mentioned differences, our results for the internal structure and composition are consistent with these previous studies. Neuenschwander et al. (2024) concluded that Uranus has a convective atmosphere and mantle on top of a non-convective inner region which is stable against convection. The transition between the convective and non-convective region was highly model-dependent and varied between ~40 and ~85% of Uranus’ radius. Neuenschwander et al. (2024) found much higher central temperatures in comparison to standard adiabatic models. Our results can be viewed as an extension of those in Neuenschwander et al. (2024). We find numerous models that agree with the aforementioned description but also infer other type of solutions, such as almost purely convective interiors. This is due to the broader range of solutions found with our random algorithm. We conclude that the work presented in this paper is compatible with the results of Neuenschwander et al. (2024) and shows even more possibilities for Uranus’ internal structure.

Podolak et al. (2022) found that the water-to-rock ratio in Uranus cannot be lower than 0.5 for their models that were analogously built to our U_3,comp case. We confirm this trend as can be seen from Figure 7: Aa low water-to-rock ratio is rare for U_3,comp models and tends to require high central temperatures, as found by Podolak et al. (2022). We note that the models of Podolak et al. (2022) were based on random density profiles which do not necessarily fit the available gravitational moments J₂ and J₄ (as only the moment of inertia was fitted). Despite the significantly different methods, the results are in agreement and should be considered to be robust.

Finally, Malamud et al. (2024) recently proposed an alternative to H-He being present in the core as a mechanism to explain a low water-to-rock ratio of Uranus and Neptune that is consistent with the expected composition of outer solar system plan-etesimals. They showed that chemical reactions between plan-etesimals dominated by organic-rich refractory materials and the hydrogen in gaseous atmospheres of protoplanets can form large amounts of methane. It was concluded that Uranus and Neptune could be compatible with having accreted refractory-dominated planetesimals, while still remaining ‘icy’ thanks to significant methane enrichment instead of water ice enrichment (see Malamud et al. 2024, for details).

4.4 Suggested future investigations

While our study presents an important advancement in empirical models, our work includes limitations and further improvements can be made in the future as listed below: • Our models are based on polytropic relations between the pressure and the density within the planet. A more agnostic approach would forego this for completely random density profiles that are consistent with the gravity data.

• Our models assume uniform rotation and no winds. In a refined study, the dynamical wind contributions to the gravity data should be accounted for. However, the difference due to winds on the overall composition is expected to be small (Neuenschwander et al. 2024).

• The atmospheric model (Hueso et al. 2020) and the EoS (More et al. 1988; Redmer et al. 2011; Chabrier et al. 2019) used for this work come with uncertainties themselves. These propagate directly to the results presented in this report.

• The β pre-factor used for our algorithm was chosen ad hoc. It remains unclear what β value would correspond to a physically realistic temperature gradient in non-convective regions. This problem directly propagates to the inferred central temperatures, which can be rather high.

• The order of magnitude of the Brunt-Väisälä frequency N² is unrealistically high. Future work should aim to decrease the typical order of magnitude further.

We suggest that a combination of both empirical (using randomness and relying primarily on gravity data) and physical (realistic EoS and density gradients) methods is a promising direction for future research. Such a semi-empirical approach would likely overcome several limitations of this work and provide both more robust and yet less biased predictions for planetary interior modelling.

5 Summary and conclusions

Empirical structure models have the advantage of probing solutions that are not covered by traditional interior models. First, we present new and improved empirical models of Uranus using four polytropes and tenth order ToF. These models are improved compared to previous work due to the more physical representation of Uranus’ atmosphere. The addition of the atmospheric polytrope increases the success rate of the random algorithm presented in this work by a factor of 10. The new random algorithm to infer the planetary temperature and composition profiles can account for two different structure models for Uranus.

We find that our empirical models would only match an adi-abatic density gradient by an extraordinary coincidence. This is because the density profiles of empirical models are, by design, only made to fit the gravitational data. Consequently, in several cases we infer super-adiabatic temperatures that result in unrealistic Brunt-Väisälä frequencies. As we discussed earlier, in several cases differences in temperature can also affect the inferred planetary composition. We therefore suggest that the Brunt-Väisälä frequency could be used as an additional constraint when inferring the planetary composition from empirical models and we hope to address this topic in future research.

Overall, several key conclusions from our study can be made:

• The tenth order ToF provides accurate predictions of the gravitational moments up to J₁₄ assuming uniform rotation: J₆ = (5.3078 ± 0.3312)⋅10⁻⁷, J₈ = (−1.1114 ± 0.1391)⋅10⁻⁸, J₁₀ = (2.8616 ± 0.5466)⋅10⁻¹⁰, J₁₂ = (−8.4684 ± 2.0889)⋅ 10⁻¹², and J₁₄ = (2.7508 ± 0.7944)⋅10⁻¹³.

• Allowing H-He to be mixed into the deeper interior (as is expected from formation models) causes the inferred composition of Uranus to be rock-dominated: SiO₂ ~ 56%, H₂O ~ 33%, H-He ~ 9%, and Fe ~ 2% of the total mass.

• We find an upper limit of ~0.1 for the H-He mass fraction in the centre.

• No H-He in deeper parts of the planet leads to lower central temperatures and necessitates a higher H₂O abundance: H₂O ~ 70%, SiO2 ~ 24%, H-He ~ 6%, and Fe ~ 1% of the total mass.

• Almost all of our Uranus models are either mostly convective except for a few sharp transition regions for ~90% of Uranus’ radius r_U or only convective in the outermost region above ~0.8 r_U.

• Almost all of our 4-polytrope based models possess a non-convective region around ~0.75 r_U.

• The majority of our models possess a region, ranging between ~(0.75–0.9)r_U, that is both convective and possess ionic H₂O. This could explain the observed magnetic activity of Uranus.

• The mean location of the convective ionic H₂O range is strongly correlated with the H-He mass fraction. It’s length can range anywhere from barely above 0, up to a maximum of ~0.l5r_U.

• An accurate measurement of J₆ and J₈ would significantly constrain the composition solution space, most notably the H-He mass fractions and consequently the mean location of the convective ionic H₂O range.

Our work suggests that a rock-dominated Uranus is a viable scenario when allowing the presence of H-He in its centre. We clearly show that Uranus’ interior is complex and is likely to harbour composition gradients and non-adiabatic regions. This affects the inferred bulk composition and importantly the inferred water-to-rock ratio. These results should be linked to formation and evolution models of Uranus.

Finally, we suggest that a dedicated space mission to Uranus is required to improve our understanding of Uranus’ internal structure and bulk composition. Ideally, such a mission should measure Uranus’ gravitational and magnetic field with high accuracy as well as the temperature and composition of its atmosphere. Such measurements will provide much-needed constraints to interior models that will eventually reveal key information on Uranus’ interior, origin, and evolution.

Data availability

A Mathematica file is available at https://zenodo.org/doi/10.5281/zenodo.13340864

Acknowledgements

The random algorithm presented in this work is an adapted version of the Master Thesis done by Sara Engeli at the Department of Astrophysics, University of Zurich. 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. We are thankful to the anonymous referee for valuable improvement suggestions and to Naor Movshovitz, his publicly available Theory of Figures code served as a foundation for our work.

Appendix A The Theory of Figures (ToF)

Consider a planet with a total potential U. ToF (Zharkov & Trubitsyn (1978)) aims at characterising surfaces r_l(ϑ) of constant U, that is, finding solutions to the equation ${U|}_{r=r_l(\vartheta )}$ = constant. r_l(ϑ) parametrises the surface and is given by a spheroid Ansatz $$r_l(\vartheta )=l\left(1+{\displaystyle \sum _{n=0}^{\infty }{s}_{2n}}(l)P_{2n}(\cos \vartheta )\right),$$(A.1)

where ϑ is the polar angle. l denotes the average radius of the corresponding spheroid, that is $$(4\pi /3)l^3=2\pi {\displaystyle {\int }_{-1}^{1}d}\cos \vartheta {\displaystyle {\int }_0^{r_l(\vartheta )}d}{r}^{\prime }{r}^{\prime 2}.$$(A.2)

Equation A.2 is called the condition of equal volume and can be used to calculate s₀ as a function of s₂, s4,.… The Legendre polynomials are denoted by P₂n(cos ϑ). s_2n(l) are called figure functions and are to be determined. U = V + Q is the sum of the Newtonian potential V = −G ∫ d³ r′ρ/|r′ − r| and the centrifugal potential Q, which is given by $Q=-\frac{1}{2}{\omega }^2{r}^2{\sin }^2(\vartheta )$ for the case of solid body rotation with a frequency ω. By inserting r_l (ϑ) into U(r) = V(r) + Q(r) and expanding every expression order by order one arrives at a general formula of the form $$U(l,\vartheta )=-\frac{4\pi }{3}G\overline{\rho }{l}^2{\displaystyle \sum _{k=0}^{\infty }{A}_{2k}}(l)P_{2k}(\mu ),$$(A.3)

where µ : = cos ϑ. The functions A_2k (l) are determined uniquely by the figure functions s_2n (l), the frequency ω, and the functions S _2n (l) defined below. By assumption, it holds U(l, ϑ) = constant as a function of A on surfaces of constant total potential. Therefore, A_2k(l) = 0 for k > 0 and U(l, A) = A₀(l) must hold on these surfaces. This conclusion will be crucial for the numerical algorithm described further down. The aforementioned expansion to arrive at equation A.3 is carried out by expanding the Newtonian potential $$V(r,\vartheta )=-\frac{G}{r}{\displaystyle \sum _{n=0}^{\infty }\left(r^{-2n}{D}_{2n}(r)+r^{2n+1}{D}_{2n}^{\prime }(r)\right)}{P}_{2n}(\mu ),$$(A.4)

where we have introduced $$\begin{array}{l}{D}_n(l)=\frac{2\pi }{n+3}{\displaystyle {\int }_0^{l}d}{l}^{\prime }\rho \left(l^{\prime }\right){\displaystyle {\int }_{-1}^{1}d}{\mu }^{\prime }{P}_n\left({\mu }^{\prime }\right)\frac{d{r}^{n+3}}{dl},\\ D_n^{\prime }(l)=\frac{2\pi }{2-n}{\displaystyle {\int }_l^{R_{\text{m}}}d}{l}^{\prime }\rho \left(l^{\prime }\right){\displaystyle {\int }_{-1}^{1}d}{\mu }^{\prime }{P}_n\left({\mu }^{\prime }\right)\frac{d{r}^{(2-n)}}{d{l}^{\prime }},(n\ne 2)\\ D_2^{\prime }(l)=2\pi {\displaystyle {\int }_l^{R_{\text{m}}}d}{l}^{\prime }\rho \left(l^{\prime }\right){\displaystyle {\int }_{-1}^{1}d}{\mu }^{\prime }{P}_2\left({\mu }^{\prime }\right)\frac{d\ln r}{d{l}^{\prime }}.\end{array}$$(A.5)

R_m denotes the mean radius of the planet, that is R_m = l_max. A dimensionless form can be obtained through the definitions $$S_n(z)=\frac{3}{4\pi \overline{\rho }{l}^{n+3}}{D}_n(l),S_n^{\prime }(z)=\frac{3}{4\pi \overline{\rho }{l}^{2-n}}{D}_n^{\prime }(l),$$(A.6)

At this point we have everything to get the functions A_2k (l) defined in equation A.3 as a function of the figure functions s_2n(l), the frequency ω, and the functions S_2n(l). However, to implement and understand the numerical algorithm, we need to do additional work: By defining r_l(ϑ) = l(1 + Σ(l, ϑ)) and $$\begin{array}{l}{S}_n(z)=\frac{1}{z^{n+3}}{\displaystyle {\int }_0^{z}d}{z}^{\prime }\frac{\rho \left(z^{\prime }\right)}{\overline{\rho }}\frac{d}{d{z}^{\prime }}\left[z^{\prime n+3}{f}_n\left(z^{\prime }\right)\right],\\ S_n^{\prime }(z)=\frac{1}{z^{(2-n)}}{\displaystyle {\int }_z^{1}d}{z}^{\prime }\frac{\rho \left(z^{\prime }\right)}{\overline{\rho }}\frac{d}{d{z}^{\prime }}\left[z^{\prime 2-n}{f}_n^{\prime }\left(z^{\prime }\right)\right],\\ S_0(z)=\frac{m(z)}{M{z}^3},\end{array}$$(A.7)

where $$\begin{array}{l}{f}_n(z)=\frac{3}{2(n+3)}{\displaystyle {\int }_{-1}^{1}d}\mu P_n(\mu ){(}^1,\\ f_n^{\prime }(z)=\frac{3}{2(2-n)}{\displaystyle {\int }_{-1}^{1}d}\mu P_n(\mu ){(}^1(n\ne 2),\\ f_2^{\prime }(z)=\frac{3}{2}{\displaystyle {\int }_{-1}^{1}d}\mu P_n(\mu )\ln (1+\text{Σ}),\end{array}$$(A.8)

we can rewrite, after performing a partial integration and assuming dρ/dz to be finite, $$\begin{array}{l}{S}_n(z)=\frac{\rho (z)}{\overline{\rho }}{f}_n(z)-\frac{1}{z^{n+3}}{\displaystyle {\int }_0^z\frac{d\rho }{\overline{\rho }}}{z}^{\prime n+3}{f}_n\left(z^{\prime }\right),\\ S_n^{\prime }(z)=-\frac{\rho (z)}{\overline{\rho }}{f}_n^{\prime }(z)+\frac{1}{z^{2-n}}\left(\frac{\rho (1)}{\overline{\rho }}{f}_n^{\prime }(1)-{\displaystyle {\int }_z^1\frac{d\rho }{\overline{\rho }}}{z}^{\prime 2-n}{f}_n^{\prime }\left(z^{\prime }\right)\right).\end{array}$$(A.9)

The integration over µ can be performed analytically for the expressions f_n(z) and f′_n(z). We now posses convenient expressions for numerical evaluation: Divide the planet into a given number of spheroids N with volumetric mean radii l₁,…, l_N. Start with a (bad) guess for the value of the figure functions s2n at positions l_i for i ∈ {1,…, N}. Equation A.8 then allows us to calculate the coefficients f_n and f′_n which in turn gives us a guess for the functions S_n and S′_n thanks to equation A.9. Now one can check the conclusion from equation A.3: If A_2k(l_i) = 0 for k > 0 and all i ∈ {1,…, N} then the (bad) guess was actually good and the figure function values for surfaces of constant potential have been found. Otherwise, one can solve A_2k (l) = 0 for k > 0 explicitly to linear order⁴ to get a better guess for the figure functions s_2n(l). This can be repeated until convergence is reached. The gravitational moments J_2n can then be obtained by comparing equation A.4 with equation 4 which yields $$J_{2n}=-{\left(\frac{R_m}{R_{\text{eq}}}\right)}^{2n}{S}_{2n}(1).$$(A.10)

A.1 Barotropic differential rotation

We need to take a few steps back in order to understand what forms of differential rotation are allowed in the ToF discussion, that is, what kind of rotational potentials Q. We start with the assumptions used in Zharkov & Trubitsyn (1978) and consider a planet consisting of a fluid with a density ρ, pressure p and velocity field υ. We use spherical coordinates r (radius), ϑ (polar angle), and φ (azimuthal angle). We assume that electromagnetic forces can be neglected and postulate

• Conservation of momentum for a static fluid with no viscosity: $$(\mathit{v}\cdot \mathbf{\nabla })\mathit{v}=-\frac{1}{\rho }\mathbf{\nabla }p+\mathbf{\nabla }V,$$(A.11)

• Conservation of mass for a static fluid: $$\mathbf{\nabla }\cdot (\rho \mathit{v})=0,$$(A.12)

• Newtonian gravity: $$V(\mathit{r})=G{\displaystyle \int \frac{\rho \left({\mathit{r}}^{\prime }\right)}{\left|\mathit{r}-{\mathit{r}}^{\prime }\right|}}\text{d}{\mathit{r}}^{\prime },$$(A.13)

• p only depends on ρ, not the temperature T (barotropic): $$p=p(\rho ),$$(A.14)

• υ only has an azimuthal component υ_φ that is independent of φ: $$\mathit{v}=v_{\phi }(r,\vartheta ){\widehat{e}}_{\phi }.$$(A.15)

Thanks to equation A.14 we can define $$P(p)={\displaystyle {\int }_{p_0}^p\frac{ \text{d}{p}^{\prime }}{\rho \left(p^{\prime }\right)}}\Rightarrow \mathbf{\nabla }P=\frac{1}{\rho }\mathbf{\nabla }p$$(A.16)

due to the Fundamental Theorem of Calculus. Consequently, we are able to rewrite Equation A.11 as $$(\mathit{v}\cdot \mathbf{\nabla })\mathit{v}=\mathbf{\nabla }(-P+V),$$(A.17)

which implies $$\mathbf{\nabla }\times ((\mathit{v}\cdot \mathbf{\nabla })\mathit{v})=0.$$(A.18)

Equation A.18 is a testable condition for any velocity field υ that is consistent with our assumptions. We can simplify equation A.18 by using equation A.15 for (υ ⋅ ∇)υ: $$\begin{array}{l}(\mathit{v}\cdot \mathbf{\nabla })\mathit{v}=\left(v_r\frac{\partial v_r}{\partial r}+\frac{v_{\vartheta }}{r}\frac{\partial v_r}{\partial \vartheta }+\frac{v_{\phi }}{r\sin \vartheta }\frac{\partial v_r}{\partial \phi }-\frac{v_{\vartheta }{v}_{\vartheta }+v_{\phi }{v}_{\phi }}{r}\right){\widehat{e}}_r\\ +\left(v_r\frac{\partial v_{\vartheta }}{\partial r}+\frac{v_{\vartheta }}{r}\frac{\partial v_{\vartheta }}{\partial \vartheta }+\frac{v_{\phi }}{r\sin \vartheta }\frac{\partial v_{\vartheta }}{\partial \phi }+\frac{v_{\vartheta }{v}_r}{r}-\frac{v_{\phi }{v}_{\phi }\cot \vartheta }{r}\right){\widehat{e}}_{\vartheta }\\ +\left(v_r\frac{\partial v_{\phi }}{\partial r}+\frac{v_{\vartheta }}{r}\frac{\partial v_{\phi }}{\partial \vartheta }+\frac{v_{\phi }}{r\sin \vartheta }\frac{\partial v_{\phi }}{\partial \phi }+\frac{v_{\phi }{v}_r}{r}+\frac{v_{\phi }{v}_{\vartheta }\cot \vartheta }{r}\right){\widehat{e}}_{\phi }\\ =\left(\frac{-v_{\phi }^2}{r}\right){\widehat{e}}_r+\left(\frac{-v_{\phi }^2\cot \vartheta }{r}\right){\widehat{e}}_{\vartheta }.\end{array}$$(A.19)

Since the curl of a vector field A in spherical coordinates is given by $$\begin{array}{cc}& \mathbf{\nabla }\times A=\frac{1}{r\sin \vartheta }\left(\frac{\partial }{\partial \vartheta }\left(A_{\phi }\sin \vartheta \right)-\frac{\partial A_{\vartheta }}{\partial \phi }\right){\widehat{e}}_r\\ & +\frac{1}{r}\left(\frac{1}{\sin \vartheta }\frac{\partial A_r}{\partial \phi }-\frac{\partial }{\partial r}\left(r{A}_{\phi }\right)\right){\widehat{e}}_{\vartheta }\\ & +\frac{1}{r}\left(\frac{\partial }{\partial r}\left(r{A}_{\vartheta }\right)-\frac{\partial A_r}{\partial \vartheta }\right){\widehat{e}}_{\phi },\end{array}$$(A.20)

equation A.18 states $$\begin{array}{cc}& \mathbf{\nabla }\times ((\mathit{v}\cdot \mathbf{\nabla })\mathit{v})=\frac{1}{r\sin \vartheta }\left(0-\frac{\partial }{\partial \phi }\left(\frac{-v_{\phi }^2\cot \vartheta }{r}\right)\right){\widehat{e}}_r\\ & +\frac{1}{r}\left(\frac{1}{\sin \vartheta }\frac{\partial }{\partial \phi }\left(\frac{-v_{\phi }^2}{r}\right)-0\right){\widehat{e}}_{\theta }\\ & +\frac{1}{r}\left(\frac{\partial }{\partial r}\left(r\left(\frac{-v_{\phi }^2\cot \vartheta }{r}\right)\right)-\frac{\partial }{\partial \vartheta }\left(\frac{-v_{\phi }^2}{r}\right)\right){\widehat{e}}_{\phi }=0.\end{array}$$(A.21)

The ê_r and ê_ϑ components are already zero thanks to the assumption that υ_φ (r, ϑ) does not depend on φ. The ê_φ component is not necessarily zero. Rather, it yields the condition $$\begin{array}{ll}\frac{\partial }{\partial r}\left(r\left(\frac{-v_{\phi }^2\cot \vartheta }{r}\right)\right)\hfill & =\frac{\partial }{\partial \vartheta }\left(\frac{-v_{\phi }^2}{r}\right)\hfill \\ \iff \frac{\partial v_{\phi }}{\partial r}\cot \vartheta \hfill & =\frac{\partial v_{\phi }}{\partial \vartheta }\frac{1}{r}\hfill \end{array}$$(A.22)

that quantifies velocity fields υ that are acceptable given our assumptions. We now aim to define a total potential U that facilitates the calculation of the pressure p given a density ρ. Plugging in the result from equation A.19 into equation A.11 component-wise, we arrive at $$\begin{array}{ll}\frac{\partial p}{\partial r}\hfill & =\rho \left(\frac{v_{\phi }^2}{r}+\frac{\partial V}{\partial r}\right)=\rho \frac{\partial }{\partial r}\left({{\displaystyle \int }}^{\text{}}\frac{v_{\phi }^2}{r}dr+V\right)\hfill \\ \frac{1}{r}\frac{\partial p}{\partial \vartheta }\hfill & =\rho \frac{1}{r}\left(v_{\phi }^2\cot \vartheta +\frac{\partial V}{\partial \vartheta }\right)=\rho \frac{1}{r}\frac{\partial }{\partial \vartheta }\left({{\displaystyle \int }}^{\text{}}{v}_{\phi }^2\cot \vartheta d\vartheta +V\right).\hfill \end{array}$$(A.23)

Thanks to the condition derived in equation A.22, these two equations are actually one: ∇p = ρ∇(V + Q), where we defined $$Q={\displaystyle {\int }_0^r\frac{v_{\phi }^2\left(r^{\prime },\vartheta \right)}{r^{\prime }}}d{r}^{\prime }={\displaystyle {\int }_0^{\vartheta }{v}_{\phi }^2}\left(r,{\vartheta }^{\prime }\right)\cot {\vartheta }^{\prime }d{\vartheta }^{\prime }.$$(A.24)

U = V + Q is called the total potential. Surfaces of constant density, pressure and potential coincide, since ∇p and ∇U are pointing in the same direction (remember ∇p = ρ∇U). Furthermore, ∇ρ = ρ∇ U ⇒ ∇ρ × V U = 0. So the normal vectors of U and ρ coincide at every point in space. One can solve $$\frac{\partial p}{\partial r}=\rho \frac{\partial U}{\partial r}$$(A.25)

to infer the pressure p as a function of the radius in the planet as soon as the total potential U and density ρ are known. To summarise, we can indeed define U = V + Q as a total potential and recover the condition for hydrostatic equilibrium, p⁻¹ ∇p = ∇U, as long as equation A.22 holds. Uniform rotation is characterised by υ_φ = ωr sin ϑ, satisfying equation A.22 and yielding $Q=\frac{1}{2}{\omega }^2{r}^2{\sin }^2\vartheta$. A more general form is given by $$Q={\omega }^2{r}^2{\sin }^2\vartheta \left(\frac{1}{2}+{\displaystyle \sum _{i=0}^{\infty }\frac{{\alpha }_{2i}}{2(i+1)}}\frac{r^{2i}}{R_m^{2i}}{\sin }^{2i}\vartheta \right),$$(A.26)

and allows for higher order terms with free and unitless parameters α_2i deviating from uniform rotation. In the literature (for example Wisdom & Hubbard (2016)), Q is often given as Q = Q₀ + ∆Q, where Q₀ = ω²r² sm²(ϑ)/2 and $$\text{Δ}Q={\displaystyle \sum _{i=1}^n{C}_{2i}}\frac{r^{2i}}{R_{\text{eq}}^{2i}}{\sin }^{2i}\vartheta .$$(A.27)

The coefficients C_2i have units of (length/time)² . These two formulations are equivalent via the relation $${\alpha }_{2i}=\frac{2(i+1)R_{\text{m}}^{2i}{C}_{2(i+1)}}{{\omega }^2{R}_{\text{eq}}^{2(i+1)}}.$$(A.28)

We now showcase the validity of our ToF implementation with barotropic differential rotation by comparing it to Wisdom & Hubbard (2016). In addition to their simple Jupiter model with solid body rotation, Wisdom & Hubbard (2016) also presented results for a differentially rotating Jupiter. They called these models DR2 and DR3. We depict the relative difference for J₂ and J₈ between our implementation and Wisdom & Hubbard (2016) for the DR2 case in Figure A.1. The DR3 case produces analogous results.

We find that, unlike in Figure 2, convergence is now stopped at a relative precision of 10⁻⁴-10⁻⁵. This is not indicative of a precision loss, but rather due to the fact that the CMS DR2 solutions are not perfectly accurate themselves. The CMS DR2 solutions were calculated using 2⁹ gridpoints and do not represent an analytical result. Our accuracy should hence exceed the accuracy of Wisdom & Hubbard (2016) for large enough N, which is consistent with a convergence halt that is visible in Figure A.1. For future reference, we therefore state our DR2 and DR3 results (using tenth order ToF with N = 2¹⁵ spheroids and n_x = N) explicitly. The calculated gravitational moments beyond J₁₄ are increasingly inaccurate, a higher order ToF implementation would be required.

Finally, we note that cubic spline interpolation is now significantly less accurate compared to the solid body rotation case depicted in Figure 2, even for high N.

A.2 Mathematica file

With this work, we provide a Mathematica file⁵ that is able to calculate all ToF coefficients up to arbitrary order. However, the runtime increases exponentially with the order of the ToF. The Mathematica file is structured as follows:

• The order up to which ToF should be calculated is defined with the parameter ToF0rder.

• The order up to which differential rotation should be accounted for is defined with the parameter DiffRotOrder, that is $$\begin{array}{l}Q=-{\omega }^2{r}^2{\sin }^2\vartheta (\frac{1}{2}\\ +{\displaystyle \sum _{i=1}^{\text{DiffRot0rder }}\frac{{\alpha }_{2i}}{2(i+1)}}\frac{r^{2i}}{R_m^{2i}}{\sin }^{2i}\vartheta ).\end{array}$$(A.29)

• The expressions $${(1+x)}^{n+3},{(1+x)}^{2-n},\log (1+x),{(\alpha +x)}^{-n-1},{(\alpha +x)}^n$$(A.30)

  found in equations A.4 and A.8 are Taylor expanded up to order ToF0rder. x and α are placeholders to be replaced later.

• The first figure function s₀ is calculated as a function of s₂, s₄, … based on the condition of equal volume (equation A.2), we replace s₀ → s₀(s₂, s₄, … ).

• The surfaces of constant total potential $$\begin{array}{l}{r}_l(\vartheta )=l(1+s_0\left(s_2,s_4,\dots \right)\\ +{\displaystyle \sum _{n=1}^{\text{ToForder }}{(\text{OR})}^n}{s}_{2n}(l)P_{2n}(\cos \vartheta ))=:l\left(1+{\text{Σ}}_l(\mu )\right)\end{array}$$(A.31)

  are defined. OR is a helper parameter to keep track of the order of an expression. Subsequently, powers of Σ_l(µ)ⁱ are calculated iteratively, only keeping terms of order (OR)ⁱ with i ≤ ToF0rder.

• The now calculated powers of Σ_l(µ)ⁱ can be inserted into the Taylor expanded expressions from equation A.30. More concretely, we replace X → (Σ_l(µ)ⁱ or lⁱ Σ_l(µ)ⁱ) and αⁱ → lⁱ according to equations A.4 and A.8.

• Finally, the necessary integrations over the angle µ are performed to recover all coefficients for the numerical algorithm.

Fig. A.1 Relative differences to the CMS DR2 solutions by Wisdom & Hubbard (2016) of the gravitational moments J2 (top) and J8 (bottom) compared to our ToF solutions. We depict results for different numbers N of used spheroids. We compare our implementation for different orders of the ToF and show the influence of using spline interpolation with the parameter nx that stands for the number of points that were calculated without interpolation.

Appendix B The Ledoux criterion

In this section, we derive the form of the Ledoux criterion used in this work. In the following derivation, g is the gravitational acceleration, and ρ, P, T, S denote density, pressure, temperature, and entropy. We use the mean molecular weight µ to parameterise the composition. The gradients are defined as $${\nabla }_{ad}:={\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}T}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}\right)}_S,{\nabla }_T:=\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}P},{\nabla }_{\mu }:=\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mu }{d\text{\hspace{0.17em}ln\hspace{0.17em}}P}.$$(B.1)

We start with the Brunt-Väisälä frequency as in Kippenhahn et al. (2013) $$\begin{array}{l}{N}^2=g\frac{-{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}T}\right)}_{P,\mu }}{-\frac{dr}{d\text{\hspace{0.17em}ln\hspace{0.17em}}P}}\left({\nabla }_{ad}-{\nabla }_T\right)+g\frac{{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\mu }\right)}_{P,T}}{-\frac{dr}{d\text{\hspace{0.17em}ln\hspace{0.17em}}P}}{\nabla }_{\mu }\\ =\frac{g}{r}\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}P}{d\text{\hspace{0.17em}ln\hspace{0.17em}}r}\left(\frac{-1}{{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}T}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}\right)}_{\rho ,\mu }{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\right)}_{T,\mu }}\left({\nabla }_{ad}-{\nabla }_T\right)-{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\mu }\right)}_{P,T}{\nabla }_{\mu }\right),\end{array}$$(B.2)

where we used the triple product rule $$-1={\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}T}\right)}_{P,\mu }{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}T}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}\right)}_{\rho ,\mu }{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\right)}_{T,\mu }.$$(B.3)

By defining $$V:=-\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}P}{d\text{\hspace{0.17em}ln\hspace{0.17em}}r},{\chi }_{\rho }:={\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\right)}_{T,\mu },{\chi }_T:={\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}T}\right)}_{\rho ,\mu },$$(B.4)

we can write $$N^2=\frac{gV}{r}\frac{{\chi }_T}{{\chi }_{\rho }}\left({\nabla }_{ad}-{\nabla }_T+\frac{{\chi }_{\rho }}{{\chi }_T}{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\mu }\right)}_{P,T}{\nabla }_{\mu }\right).$$(B.5)

In hydrostatic equilibrium ρ⁻¹∇P = ∇U and therefore V is equal to GM(r)ρ/(rP). With this we arrive at the following expression (see also Unno et al. (1989)): $$N^2=\frac{g^2\rho }{P}\frac{{\chi }_T}{{\chi }_{\rho }}\left({\nabla }_{ad}-{\nabla }_T+B\right),$$(B.6)

where the composition term B is $$\begin{array}{l}\begin{array}{ll}B\hfill & =\frac{{\chi }_{\rho }}{{\chi }_T}{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\mu }\right)}_{P,T}\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mu }{d\text{\hspace{0.17em}ln\hspace{0.17em}}p}\hfill \\ \hfill & =\frac{{\chi }_{\rho }}{{\chi }_T}{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\mu }\right)}_{P,T}\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mu }{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{d\text{\hspace{0.17em}ln\hspace{0.17em}}p}\hfill \end{array}\\ =\frac{{\chi }_{\rho }}{{\chi }_T}\underset{\delta \text{\hspace{0.17em}ln\hspace{0.17em}}\rho \to 0}{\lim }\frac{\text{\hspace{0.17em}ln\hspace{0.17em}}\rho \left(P,T,\mu +\frac{d\mu }{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\delta \text{\hspace{0.17em}ln\hspace{0.17em}}\rho \right)-\text{\hspace{0.17em}ln\hspace{0.17em}}\rho (P,T,\mu )}{\delta \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{d\text{\hspace{0.17em}ln\hspace{0.17em}}p}.\end{array}$$(B.7)

The last line equality follows from $$\begin{array}{l}\text{\hspace{0.17em}ln\hspace{0.17em}}\rho \left(P,T,\mu +\frac{d\mu }{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\delta \text{\hspace{0.17em}ln\hspace{0.17em}}\rho \right)=\text{\hspace{0.17em}ln\hspace{0.17em}}\rho (P,T,\mu )+\frac{d\mu }{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\partial \mu }\right)}_{P,T}\delta \text{\hspace{0.17em}ln\hspace{0.17em}}\rho \\ +𝒪\left(\delta \text{\hspace{0.17em}ln\hspace{0.17em}}{\rho }^2\right)\\ =\text{\hspace{0.17em}ln\hspace{0.17em}}\rho (P,T,\mu )+\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mu }{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\left(\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}\mu }\right)}_{P,T}\delta \text{\hspace{0.17em}ln\hspace{0.17em}}\rho \\ +𝒪\left(\delta \text{\hspace{0.17em}ln\hspace{0.17em}}{\rho }^2\right).\end{array}$$(B.8)

Therefore, B can be evaluated numerically as $$\begin{array}{l}B=\frac{{\chi }_{\rho }}{{\chi }_T}\underset{\delta \text{\hspace{0.17em}ln\hspace{0.17em}}\rho \to 0}{\lim }\frac{\text{\hspace{0.17em}ln\hspace{0.17em}}\rho \left(P,T,\mu +\frac{d\mu }{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\delta \text{\hspace{0.17em}ln\hspace{0.17em}}\rho \right)-\text{\hspace{0.17em}ln\hspace{0.17em}}\rho (P,T,\mu )}{\delta \text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }{d\text{\hspace{0.17em}ln\hspace{0.17em}}p}\\ =\frac{{\chi }_{\rho }}{{\chi }_T}\frac{\text{\hspace{0.17em}ln\hspace{0.17em}}\rho \left(P_n,T_n,{\mu }_{n+\text{Δ}n}\right)-\text{\hspace{0.17em}ln\hspace{0.17em}}\rho \left(P_n,T_n,{\mu }_n\right)}{\text{\hspace{0.17em}ln\hspace{0.17em}}{\rho }_{n+\text{Δ}n}-\text{\hspace{0.17em}ln\hspace{0.17em}}{\rho }_n}\frac{\text{\hspace{0.17em}ln\hspace{0.17em}}{\rho }_{n+\text{Δ}n}-\text{\hspace{0.17em}ln\hspace{0.17em}}{\rho }_n}{\text{\hspace{0.17em}ln\hspace{0.17em}}{P}_{n+\text{Δ}n}-\text{\hspace{0.17em}ln\hspace{0.17em}}{P}_n}\\ =\frac{{\chi }_{\rho }}{{\chi }_T}\frac{\text{\hspace{0.17em}ln\hspace{0.17em}}\rho \left(P_n,T_n,{\mu }_{n+\text{Δ}n}\right)-\text{\hspace{0.17em}ln\hspace{0.17em}}\rho \left(P_n,T_n,{\mu }_n\right)}{\text{\hspace{0.17em}ln\hspace{0.17em}}{P}_{n+\text{Δ}n}-\text{\hspace{0.17em}ln\hspace{0.17em}}{P}_n}.\end{array}$$(B.9)

Appendix C Further data and results

Figure C.1 shows the pressure P of our empirical interior models as a function of their density ρ. We depict the same zoomedin subplot that is used in Figure 1 and now include the whole solution range.

Figure C.2 displays a direct comparison of the density profiles depicted in Figures 3.2 (blue) and 3.2 (orange).

We find that although the density profiles are indeed somewhat similar, they are distinct enough to justify the fact that the random algorithm infers significantly different temperature and composition profiles for the two of them. The blue density profile (from Figure 3.2) has a core density that is ∼ 15% lower and a less steep density gradient compared to the orange density profile (from Figure 3.2).

Figures C.3 and C.4 present all correlation plots of our results for the U_3,comp and U_4,comp case, respectively. The plots involving the convective ionic H₂ O range only include models that posses a convective ionic H₂O region in the first place. These are 77% and 83% of solutions, respectively.

Fig. C.1 Uranus’ pressure P as a function of the density ρ. Other pressure-density profiles that have been published by previous studies are also presented for comparison (Helled et al. 2011; Nettelmann et al. 2013; Vazan & Helled 2020). The panel on the right shows the outermost region. The white dots show the atmosphere model by Hueso et al. (2020).

Fig. C.2 Direct comparison of the density profiles depicted in Figures 3.2 (blue) and 3.2 (orange).

Fig. C.3 Summarised results of the random algorithm for 4-polytrope empirical models of Uranus (U3,comp case). Each individual point represents one empirical profile. There are no errorbars indicating the range for a single profile to enhance readability.

Fig. C.4 Summarised results of the random algorithm for 4-polytrope empirical models of Uranus (U4,comp case). Each individual point represents one empirical profile. There are no errorbars indicating the range for a single profile to enhance readability.

References

All Tables

All Figures

Fig. 1 Uranus’ density ρ as a function of the average normalised radius r/R. Other density profiles that have been published by previous studies are also presented for comparison (Helled et al. 2011; Nettelmann et al. 2013; Vazan & Helled 2020). The panel on the right show the pressure as a function of density P(ρ) in the outermost region. The white dots show the atmosphere model by Hueso et al. (2020). The full P(ρ) distribution is presented in appendix C in Figure C.1.

In the text

Fig. 2 Relative difference of the Bessel gravitational moment solutions J2 (top) and J8 (bottom) compared to the ToF solutions. Results are shown for different numbers N of used spheroids. We compare our implementation (This work) to the one used in Movshovitz & Fortney (2022) (M&F 22) for different orders of the ToF. We also show the influence of using cubic spline interpolation with the parameter nx that stands for the number of points that were calculated without interpolation.

In the text

	Fig. 3 Runtime in seconds of the calculations shown in Figure 2 as a function of the number N of used spheroids in the ToF. We refer to the legend and caption of Figure 2 for further explanations.
In the text

	Fig. 4 Our algorithm for inferring the temperature profile and composition of the empirical density (and pressure) profiles.
In the text

Fig. 5 Schematic drawing of the two possible composition models considered in this work. Top: U3,comp model that allows the presence of up to three elements simultaneously in a given layer. Bottom: U4,comp model that allows the presence of up to four elements simultaneously in a given layer. The ratios between different components are chosen freely and composition gradients are allowed but not necessary. The transition radii (white gaps) vary and are placed by the algorithm.

In the text

Fig. 6 Three example empirical density profiles of Uranus (depicted from left to right). Their density, composition, and temperature profile inferred by the new random algorithm is depicted from top to bottom, respectively. The extended phase diagram of water based on Redmer et al. (2011) is shown in addition to the temperature profiles. The stability against convection is presented using the Brunt-Väisälä frequency N2 in the bottom last plots.

In the text

Fig. 7 Inferred composition and temperature for our Uranus empirical models for U3,comp (top) and U4,comp (bottom). The x-axis corresponds to the mass fraction of a material and the y-axis shows the central temperature. The colours represent different materials (see legend of Figure 6) and the error bars indicate the 1σ range of solutions found by the random algorithm for each density profile. To improve readability, error bars for the central temperature Tc are not shown. They are included in Figure 11.

In the text

	Fig. 8 Uranus’ inferred composition using the new random algorithm. We depict results for the U3,comp (top) and U4,comp (bottom) models. The pie charts summarise the composition histograms from Figure 7, and the numbers are the means and 1 σ ranges of all empirical models shown as dots in Figure 7.
In the text

Fig. 9 First and third plot: histograms showing the percentage of Uranus’ radius that is convective versus the number of solutions. In the U4,comp case, more than 100 empirical profiles were inferred to be convective for ~90% of Uranus’ radius. Second and last plot: radial distributions of the convectivity for all solutions. In both cases, only ~20% of the profiles were found to be convective at ~0.75rU.

In the text

	Fig. 10 H-He mass in Uranus versus mean location of the convective ionic H2O region. The error bars indicate the 1σ range of solutions found by the algorithm for density profiles possessing an ionic-water region. Profiles without a convective ionic H2O region are not presented.
In the text

	Fig. 11 H-He mass fraction in Uranus’ centre versus the central temperature Tc. The error bars indicate the 1σ range of solutions found by our algorithm.
In the text

	Fig. 12 Total H-He mass fraction MX+Y/M versus the higher order gravitational moments J6 (top) and J8 (bottom). The errorbars indicate the 1σ range of solutions found by the random algorithm for each density profile.
In the text

	Fig. 13 Absolute change in composition (H-He, water, rock, and iron) when enforcing a perfect adiabatic temperature profile. Positive numbers indicate that the forced adiabatic temperature profiles possess more of a given component. The results are shown for both U3,comp (top) and U4,comp (bottom) models.
In the text

	Fig. 14 Density as a function of the temperature for three isobars at 100, 500, and 1000 GPa, with the pressure increasing from the bottom to the top of the figures. The four panels show the density-temperature relation of a protosolar H-He mixture, pure water, rock, and iron.
In the text

Fig. A.1 Relative differences to the CMS DR2 solutions by Wisdom & Hubbard (2016) of the gravitational moments J2 (top) and J8 (bottom) compared to our ToF solutions. We depict results for different numbers N of used spheroids. We compare our implementation for different orders of the ToF and show the influence of using spline interpolation with the parameter nx that stands for the number of points that were calculated without interpolation.

In the text

	Fig. C.1 Uranus’ pressure P as a function of the density ρ. Other pressure-density profiles that have been published by previous studies are also presented for comparison (Helled et al. 2011; Nettelmann et al. 2013; Vazan & Helled 2020). The panel on the right shows the outermost region. The white dots show the atmosphere model by Hueso et al. (2020).
In the text

	Fig. C.2 Direct comparison of the density profiles depicted in Figures 3.2 (blue) and 3.2 (orange).
In the text

	Fig. C.3 Summarised results of the random algorithm for 4-polytrope empirical models of Uranus (U3,comp case). Each individual point represents one empirical profile. There are no errorbars indicating the range for a single profile to enhance readability.
In the text

	Fig. C.4 Summarised results of the random algorithm for 4-polytrope empirical models of Uranus (U4,comp case). Each individual point represents one empirical profile. There are no errorbars indicating the range for a single profile to enhance readability.
In the text