ABSTRACT (original article): We present improved empirical density profiles of Uranus and interpret them in terms of their temperature and composition using a new random algorithm. The algorithm to determine the temperature and composition is agnostic with respect to the temperature gradient in non-isentropic regions and chooses randomly amongst all possible gradients that are stable against convection and correspond to an Equation of State compatible composition. Our empirical models are based on an efficient implementation of the Theory of Figures up to 10th order including a proper treatment of the atmosphere. The accuracy of 10th order ToF enables us to present accurate calculations of the gravitational moments of Uranus up to
J
14
:
J
6
=(5.3078±0.3312)·
-7
10
,
J
8
=(-1.1114±0.1391)·
-8
10
,
J
10
=(2.8616±0.5466)·
-10
10
,
J
12
=(-8.4684±2.0889)·
-12
10
and
J
14
=(2.7508±0.7944)·
-13
10
. We consider two interior models of Uranus that differ with respect to the maximal number of materials allowed per layer of Uranus (three vs. four composition components). The case with three materials does not allow Hydrogen and Helium in deeper parts of Uranus and results in a higher water abundance which leads to lower central temperatures. On the other hand, the models with four materials allow H-He to be mixed into the deeper interior and lead to rock-dominated solutions. We find that these four composition components models are less reliable due to the underlying empirical models incompatibility with realistic Brunt frequencies. Most of our models are found to be either purely convective with the exception of boundary layers, or only convective in the outermost region. Almost all of our models possess a region that is convective and consists of ionic HO which could explain the generation of Uranus’ magnetic field. CITATION (original article): Luca Morf, Simon Müller, Ravit Helled, The interior of Uranus - Thermal profile, bulk composition and the distribution of rock, water and hydrogen and helium, A&A, 690 (2024) A105. https://doi.org/10.1051/0004-6361/202450698

Introduction

Consider a planet with a total potential
U
that has a north-south and an azimuthal symmetry. The Theory of Figures (Zharkov and Trubitsyn, 1978) aims at characterising surfaces
r
l
(ϑ)
of constant
U
, that is, finding solutions to the equation
U

r=
r
l
(ϑ)
=constant
. We make an Ansatz
r
l
(ϑ)
that parametrises these surfaces and is given by
r
l
(ϑ)=l1+
∞
∑
n=0
s
2n
(l)
P
2n
(cosϑ)
,
where
ϑ∈[0π]
is the polar angle.
l
denotes the average radius of a corresponding spheroid, that is
4π
3
3
l
=2π
1
∫
-1
dcosϑ
r
l
(ϑ)
∫
0
d
′
r
2
′
r
.
Each
l
therefore uniquely identifies a spheroid. The above equation is called the condition of equal volume and can in particular be used to calculate
s
0
as a function of
s
2
,
s
4
,….
The Legendre polynomials are denoted by
P
2n
(cosϑ)
.
s
2n
(l)
are called figure functions and are to be determined in order to ensure
U

r=
r
l
(ϑ)
=constant
.
U=V+Q
is the sum of the Newtonian potential
V=-G∫
3
d
′
r
ρ
′
r
-r
and the centrifugal potential
Q
, which is given by
Q=-
1
2
2
ω
2
r
2
sin
(ϑ)
for the case of solid body rotation with a frequency
ω
.
By inserting
r
l
(ϑ)
into
U=V+Q
and expanding every expression order by order one arrives at a general formula of the form
U(l,ϑ)-
4π
3
G
ρ
2
l
∞
∑
n0
A
2n
(l)
P
2n
(cosϑ)
.
We will calculate the functions
A
2n
(l)
below. Remember that by assumption,
U(l,ϑ)constant
as a function of
ϑ
must hold on surfaces of constant total potential. Therefore,
A
2n
(l)=0
for
n>0
and
U(l,ϑ)=
A
0
(l)
must be true on these surfaces. This conclusion is crucial: We can test if we have found a solution for the figure functions
s
2n
(l)
by checking if the resulting
A
2n
(l)
are indeed 0 for
n>0
.
Finally, note that instead of solid body rotation, one can also consider a more general form given by
Q
2
ω
2
r
2
sin
ϑ
1
2
+
∞
∑
i0
α
2i
2(i+1)
2i
r
2i
R
m
2i
sin
ϑ
and allow for higher order terms with free and unitless parameters
α
2i
deviating from uniform rotation. This is called (barotropic) differential rotation. Note that
R
m
=
l
max
=max{l:ρ(
r
l
(ϑ))>0∀ϑ}
.

Definitions and preliminary calculations

Set the following variable to True to generate .txt files with all results in a notation compatible with numerical implementation:
In[]:=
SaveFiles=False;
If you work with a .m file, replace NotebookDirectory[] with DirectoryName[$InputFileName], then type echo “wolframscript -file ToF.m” | at now to run it on Linux detached from the terminal.
We cannot expand sums to infinity, so we need to set limits for the above introduced sums:
In[]:=
ToFOrder=4;​​DiffRotOrder=0;
Now we introduce variables for the average density
ρ
av
and a dimensionless alternative
z
to the average spheroid radii
l
:
In[]:=
ρ
av
=
M
4π
3
3
Rm
;​​l=z*Rm;
Finally, we define lists for the figure functions, differential rotation parameters, etc. to display our results. These lists are not relevant for the calculations:
In[]:=
ss=List[];SS=List[];SSdash=List[];alphas=List[];​​For[i=0,i<2*ToFOrder+1,i++,ss=Append[ss,ToExpression["s"<>ToString[i]]]];​​For[i=0,i<2*ToFOrder+1,i++,SS=Append[SS,ToExpression["S"<>ToString[i]]]];​​For[i=0,i<2*ToFOrder+1,i++,SSdash=Append[SSdash,ToExpression["S"<>ToString[i]<>"p"]]];​​For[i=0,i<DiffRotOrder+1,i++,alphas=Append[alphas,ToExpression["alpha"<>ToString[i]]]];
For a north-south and azimuthally symmetric planet, the Newtonian potential can be expanded as (see for example Zharkov and Trubitsyn, 1978)​
,where we have introduced
and​
,​
,​
.​
R
m
denotes the mean radius of the planet, that is
R
m
=
l
max
=max{l:ρ(
r
l
(ϑ))>0∀ϑ}
. Note that we use a one-dimensional density function
instead of
). This is justified by assuming hydrostatic equilibrium
(
denotes the pressure). Hydrostatic equilibrium implies that surfaces of constant density and potential coincide:
, so the normal vectors of
U
and
coincide at every point in space. Therefore, we can assume that just like surfaces of constant potential are characterized by the variable
only
(U=U(l))
, the same holds true for
ρ
(and even
P
).
A dimensionless form of the above expressions can be obtained through the definitions​
, ​
,where
and
.
At this point we have everything to get the functions
A
2n
(l)
introduced above as a function of the figure functions
s
2n
(l)
, the frequency ω, and the functions
S
2n
(l)
. However, to implement a numerical algorithm, we need to do additional work: By defining
and​
,​
,​
,where​
,​
,​
,we can rewrite, after performing a partial integration and assuming
dρ
dz
to be finite,​
,​
.This form is more convenient for numerical evaluation. Note that the integration over
can be performed analytically for the expressions
and
.
We define the expressions
and
,
x
will later be replaced with
:
In[]:=
PowerLawfn[n_,x_]=Normal[Series[
n+3
(1+x)
,{x,0,ToFOrder}]]/.
β_
x
->f[x,β];​​PowerLawfdashn[n_,x_]=Normal[Series[
2-n
(1+x)
,{x,0,ToFOrder}]]/.
β_
x
->f[x,β];​​PowerLawfdash2[x_]=Normal[Series[Log[1+x],{x,0,ToFOrder}]]/.
β_
x
->f[x,β];
We also define the expressions in front of
and
in
.
α
will later be replaced with
.
x
will later be replaced with
:
In[]:=
PowerLawDn[n_,x_,α_]=Normal[Series[
-n-1
(α+x)
,{x,0,ToFOrder}]]/.
β_
x
->f[x,β];​​PowerLawDdashn[n_,x_,α_]=Normal[Series[
n
(α+x)
,{x,0,ToFOrder}]]/.
β_
x
->f[x,β];​​
We define
for the first time.
OR
is a helper parameter that we use to only calculate expressions up to the relevant order in order to save time.
s
0
will get replaced soon, therefore no
OR
parameter for
s
0
:
In[]:=
rpre[l_,μ_]=l(1+s[0,l]+Sum[
n
OR
*s[2n,l]*LegendreP[2n,μ],{n,1,ToFOrder}]);
We can already calculate the ratio between the equatorial radius (
) and the mean radius
for each spheroid:
In[]:=
Print"
R
equatorial

R
mean
=",rpre[l,0]l/.OR->1/.s[α_,zRm]->Indexed[ss,α+1];
R
equatorial

R
mean
=1+s0-
s2
2
+
3s4
8
-
5s6
16
+
35s8
128
Notation compatible with numerical implementation:
In[]:=
Print"
R
equatorial

R
mean
=",FortranForm[rpre[l,0]/l/.OR->1/.s[α_,zRm]->Indexed[ss,α+1]];​​If[SaveFiles,Export[FileNameJoin[{NotebookDirectory[],"R_ratio.txt"}],FortranForm[rpre[l,0]/l/.OR->1/.s[α_,zRm]->Indexed[ss,α+1]],"Text"]];
Notation compatible with numerical implementation:
We use our results until this point and insert them into earlier definitions:
We define storage lists for the correct labels:
Notation compatible with numerical implementation:
Notation compatible with numerical implementation:
We define storage lists for the correct labels:
Notation compatible with numerical implementation:
Notation compatible with numerical implementation:

Numerical algorithm

CITE THIS NOTEBOOK

The interior of Uranus: Thermal profile, bulk composition and the distribution of rock, water and hydrogen and helium​
by Luca Morf, Simon Müller, Ravit Helled
Wolfram Community, STAFF PICKS, August 21, 2025
​https://community.wolfram.com/groups/-/m/t/3534064