# Mass–Radius Relationships and Contraction of Condensed Planets by Cooling or Despinning

Yanick Ricard, Frédéric Chambat

► To cite this version:

Yanick Ricard, Frédéric Chambat. Mass–Radius Relationships and Contraction of Condensed Planets by Cooling or Despinning. The Astrophysical Journal, 2024, 967 (2), pp.163. 10.3847/1538-4357/ad4113. hal-04811320

HAL Id: hal-04811320

https://hal.science/hal-04811320v1

Submitted on 30 Nov 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Mass–Radius Relationships and Contraction of Condensed Planets by Cooling or Despinning

Yanick Ricard and Frédéric Chambat

Ecole normale supérieure de Lyon, CNRS, 15 parvis René Descartes, Lyon, F-69007, France; Ricard@ens-lyon.fr, Frederic.Chambat@ens-lyon.fr

Received 2023 November 25; revised 2024 April 14; accepted 2024 April 18; published 2024 May 30

## Abstract

Condensed planets contract or expand as their temperature changes. With the exception of the effect of phase changes, this phenomenon is generally interpreted as being solely related to the thermal expansivity of the planet's components. However, changes in density affect pressure and gravity and, consequently, the planet's compressibility. A planet's radius is also linked to its rate of rotation. Here again, changes in pressure, gravity, and compressibility are coupled. In this article we clarify how the radius of a condensed planet changes with temperature and rotation, using a simple and rigorous thermodynamic model. We consider condensed materials to obey a simple equation of state which generalizes a polytopic EoS as temperature varies. Using this equation, we build simple models of condensed planet's interiors including exoplanets, derive their mass–radius relationships, and study the dependence of their radius on temperature and rotation rate. We show that it depends crucially on the value of $\rho_s\delta R/K_s$ ($\rho_s$ being surface density, $g$ gravity, $R$ radius, $K_s$ surface incompressibility). This nondimensional number is also the ratio of the dissipation number which appears in compressible convection and the Gruneisen mineralogic parameter. While the radius of small planets depends on temperature, this is not the case for large planets with large dissipation numbers; Earth and a super-Earth like CoRoT-7b are in something of an intermediate state, with a moderately temperature-dependent radius. Similarly, while the radius of these two planets is a function of their rotation rates, this is not the case for smaller or larger planets.

Unified Astronomy Thesaurus concepts: Solar system terrestrial planets (797); Planetary thermal histories (2290); Extrasolar rocky planets (511)

## 1. Introduction

In the 19th century, the contraction of our planet during its secular cooling was sometimes invoked to explain the Earth's topography (Dana 1847). This interpretation was abandoned after the advent of plate tectonic theory, and mountain ranges have since been explained by plate collisions and interactions. However, for other planets and satellites in the solar system that do not appear to have plate tectonics and have not undergone a resurfacing event like appears to be the case for Venus, many features such as compression scarps are attributed to a reduction in the planet's radius. Conversely, the absence of obvious signs of planet compression or extension has been used to put constraints on the planet's thermal evolution. The amount of planet radius variation has been discussed in the case of Mercury (3–7 km of contraction, e.g., Byrne et al. 2014), Mars (0–4 km of contraction since the Early Noachian, e.g., Nahm & Schultz 2011) or the Moon (negligible radius variation, e.g., Solomon & Chaiken 1976).

The change in the radius of a planet can have various origins. One obvious cause of contraction is thermal cooling, i.e., the increase in density of the mineralogical phases when the temperature decreases. Another possible contraction linked to thermal cooling is the phase changes that can occur. For example, the crystallization of a liquid metallic core to form a denser solid inner core is necessarily accompanied by a decrease in planetary radius. Another cause may be a change in angular velocity. The despinning of a planet has two effects.

The first reduces the hydrostatic flattening (Chambat et al. 2010) and induces a stress pattern that favors lithospheric cracks oriented parallel to the rotation axis. Such a preferential orientation of faults has been invoked on Mercury (Melosh & McKinnon 1988). This early weakening of the lithosphere may later be reactivated by thermal cooling (Watters et al. 1998; Byrne et al. 2014). The second effect contracts the planet by decreasing the spherically averaged centrifugal force (Saito 1974).

In a convecting planet, the temperature is controlled by the adiabatic gradient, and cooling or heating occurs at all depths. A simple relation between a planet's cooling rate and the variation $\delta R$ of its radius $R$ is often used, for example

$$ \frac{\delta R}{R} = \frac{1}{R^3} \int_{0}^{R} \alpha(r) \delta T(r) r^2 dr, \qquad (1) $$

which takes into account the variations of thermal expansivity $\alpha$ and temperature $T$, with depth (Solomon & Chaiken 1976), or even

$$ \frac{\delta R}{R} = \frac{1}{3} \bar{\alpha} \overline{\delta T}, \qquad (2) $$

where $\bar{\alpha}$ and $\overline{\delta T}$ are uniform (Hauck et al. 2004). For the Earth, an estimate of the cooling can be obtained from petrologic observation of ancient rocks. For example, the composition of non-arc basaltic rocks as a function of age, suggests a cooling rate of 50–100 K Gyr$^{-1}$ (Herzberg et al. 2010). The thermal expansivity is decreasing with depth in the mantle but a typical value of $\alpha = 3 \times 10^{-5}$ K$^{-1}$ would imply a radius reduction of 10–20 km in 3 Gyr.

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
These relations (1)–(2) are in fact problematic. Choosing an appropriate thermal expansivity for a quantity that varies with time and depth is the first difficulty. But the real difficulty lies elsewhere: changing the temperature affects the density, the gravity, and therefore the pressure. It is well known that for large objects, density is mainly related to pressure, rather than temperature, so it is not obvious that the change in radius of a cooling planet is so directly related to temperature. Jaupart et al. (2015) give an approximated estimate of the pressure change influence by supposing a small and uniform compressibility. Of course, a model that gives us the radial variations of thermal expansion, incompressibility, density, and temperature as functions of depth, makes it possible to accurately compute the radius change as a function of temperature by perturbing the planet’s elastogravitational equations but such a model is only known for the Earth.

To calculate the variations of the Earth’s radius with the variations of its angular velocity $\Omega$, Saito (1974) has applied this method of perturbing the elastogravitational equations to obtain

$$ \frac{\delta R}{R} = \frac{2}{3} h_0^* \frac{\Omega^2 R}{g} \frac{\delta \Omega}{\Omega}, \quad (3) $$

where $g$ being Earth’s surface gravity and $h_0^*$ is a degree 0 rotational Love number estimated to be $h_0^* = 0.09835$ using a radial model very close to PREM.

Is seems therefore that precise calculation of the contraction of a planet by cooling or despinning is only possible for the Earth. In this paper, we aim to show that a realistic estimate can, however, be obtained using a simple model for a generic condensed planet.

## 2. A Lane–Emden Planet

### 2.1. Equation of State

To get an estimate of the relation between the radius and temperature of a cooling planet, we must first choose an equation of state (EoS) to describe the thermodynamic relation between a planet’s pressure $P$, temperature $T$, and density $\rho$. For a condensed planet (solid or liquid, silicate of metal), Murnaghan’s EoS (Murnaghan 1937) gives a fairly simple and versatile expression that fits well with various high-pressure, high-temperature experiments on silicates and metals, as well as with the radius properties of the Earth (Ricard et al. 2022; Ricard & Alboussière 2023):

$$ P = \frac{K_0}{p} \left[ \left( \frac{\rho}{\rho_0} \right)^p - 1 \right] + \alpha_0 K_0 (T - T_0). \quad (4) $$

In this expression $\alpha_0$, $K_0$, $\rho_0$ and $T_0$ are the thermal expansivity, isothermal incompressibility, density and temperature under reference conditions (e.g., zero pressure and 25°C), and $p$ a nondimensional empirical exponent. This EoS leads to simple expressions for thermal expansivity and isothermal incompressibility, which are related only to density

$$ \alpha = \alpha_0 \left( \frac{\rho_0}{\rho} \right)^p, \quad K_T = K_0 \left( \frac{\rho}{\rho_0} \right)^p, \quad (5) $$

and these relations are reasonably verified experimentally (Anderson 1979).

When $T$ is uniform, this Equation (4) belongs to the family of polytropic EoS, on the form

$$ P = a \rho^{1+\frac{1}{n}} + P_b, \quad (6) $$

where $n$ is called the polytropic index (the exponent $p$ in (4) is therefore $1+1/n$). Polytropic equations have been extensively used in astrophysical literature since Chandrasekhar (1958) where the focus is often on stars and gaseous planets, and the index $n$ is large ($n = +\infty$ for a perfect gas). Neutron stars are well modeled with a polytropic index close to but lower than 1. In a solid planet $p \approx 3-4$ or $n \approx 1/3-1/2$ (Stixrude & Lithgow-Bertelloni 2005). The exponent $p$ is also $\partial K_T / \partial P$ and therefore it quantifies the increase of incompressibility with pressure.

### 2.2. Density Equation

If we consider a nonrotating planet with spherical symmetry (we account for the planet rotation in Section 5), its density verifies an equation that simply reformulates the gravity equation (Poisson’s equation)

$$ 4\pi G\rho = -\nabla \cdot \mathbf{g} = -\nabla \cdot \left(\frac{\nabla P}{\rho}\right) = -\nabla \cdot \left(K \frac{\nabla \rho}{\rho^2}\right). \quad (7) $$

In this expression, the quantity $K$ is defined by

$$ K = \rho \frac{dP}{d\rho}. \quad (8) $$

If the planet is convecting, its state is close to an adiabatic state with an isentropic incompressibility $K_S$, and $K$ can therefore be identified with $K_S$ (Bullen 1975). However, for a condensed planet, $K_S$ is very close to $K_T$, and similarly, the two heat capacities $C_P$ and $C_V$ are comparable. This is a consequence of the thermodynamic rules (Mayer’s relation) that imply $K_S/K_T = C_P/C_V = 1 + \Gamma\alpha T$ where the Grüneisen parameter $\Gamma = \alpha K_T / \rho C_V$ is of order 1 for solid silicates, liquid metal (and even for ideal gases), and decreases slightly with depth, while $\alpha T$ is a small quantity of order $10^{-2}$ (in planets, $T$ increases with depth but $\alpha$ decreases more, see Equation (5)). Whether the planet is convecting or not, the influence of temperature on density remains negligible compared to that of pressure and we can therefore identify $K$ with $K_T$ or $K_S$. Similarly, we confuse the heat capacities $C_P$ and $C_V$ we later note $C$.

### 2.3. Nondimensional Variables

We now define a dimensionless density $\tilde{\rho} = \rho / \rho_c$ and a dimensionless radius $\tilde{r} = r / r_c$ where $\rho_c$ is the density at the center of the planet (in the paper, the subscript $c$ will always refer to the properties at the center, the subscript $s$ at the surface, the subscript 0 to the reference values of the EoS (4), and the tilde symbol will refer to variables without dimensions) and $r_c$ is defined by:

$$ r_c^2 = r_0^2 \left( \frac{\rho_c}{\rho_0} \right)^{p-2} \quad (9) $$

where the length $r_0$ is

$$ r_0^2 = \frac{K_0}{4\pi G \rho_0^2}. \quad (10) $$
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The nondimensional density profiles $\tilde{\rho}(r/r_c)$, for $p=2$ (solid black), $p=3$ (solid red), and $p=4$ (solid green).</figcaption></figure>

Using these variables, Equation (7) becomes

$$ \frac{1}{\tilde{r}^2} \frac{d}{d\tilde{r}} \left( \tilde{r}^2 \tilde{\rho}^{p-2} \frac{d\tilde{\rho}}{d\tilde{r}} \right) + \tilde{\rho} = 0. \quad (11) $$

This equation must be solved with $\tilde{\rho}(0) = 1$ and $d\tilde{\rho}/d\tilde{r}(0) = 0$, until the outer nondimensional radius $\tilde{r} = \tilde{R}$, where

$$ \tilde{R} = \frac{R}{r_c}. \quad (12) $$

The resolution of (11) is analytical for $p=2$, in which case the solution is the first arc of a sinc function (sinc(x) is (sin x)/x) and for $p=\infty$ and $p=6/5$ (i.e., for a polytropic index 0 or 5) which are not in the appropriate range of parameters for a condensed planet. For the other exponents $p$, the dimensionless Equation (11) can be easily solved numerically using a Runge-Kutta method. The solutions of (11) as a function of $\tilde{r}$ are shown in Figure 1 for $p=2, 3$, and $4$. These curves are computed until $\tilde{\rho} = 0$ but the density with dimensions are only related to these curves for $\tilde{r} \le \tilde{R}$, i.e., for $\tilde{\rho}(\tilde{r}) \ge \rho_s/\rho_c$ where $\rho_s$ is the surface density.

Notice that, following Chandrasekhar (1958) the density of stellar objects in astrophysical literature is often sought in the form $\tilde{\rho} = \theta^n = \theta^{\frac{1}{p-1}}$ (Horedt 2004). In which case (11) becomes the well-known Lane–Emden equation

$$ \nabla^2 \theta + \theta^n = 0 \quad (13) $$

when the radius in the Laplace operator $\nabla^2$ is normalized by $r_c/\sqrt{p-1}$. However, solving numerically the Lane–Emden Equation (13) or the Equation (11) presents the same degree of difficulty and we thought that the planetary physicists might prefer thinking in terms of $\rho$ than $\theta$.

# 3. Density with Physical Dimension

## 3.1. Input Parameters

To obtain a solution with physical dimensions, various quantities need to be specified. We need to know the central density $\rho_c$, the planet radius $R$, and the nondimensionalizing length $r_c$, but the choice of $r_c$ is related to that of two EoS parameters $K_0$ and $\rho_0$. Therefore, four quantities need to be specified, ($\rho_c$, $R$, $K_0$, $\rho_0$). Alternatively, we can choose another

set of four independent parameters from which ($\rho_c$, $R$, $K_0$, $\rho_0$) can be derived. For example, we will calculate the mass of planets when ($\rho_s$, $R$, $K_0$, $\rho_0$) are chosen, i.e., when the planet's surface density and radius $R$ are specified. To apply our model to the Earth or to all planets for which mass and radius are known, we will instead choose ($M$, $R$, $K_0$, $\rho_0$) and calculate the appropriate values of the central and surface densities. Finally, in the paragraph where we discuss the effects of rotation, ($\rho_s$, $M$, $K_0$, $\rho_0$) will be fixed values, and $R$ will vary with the rate of rotation.

When ($\rho_s$, $R$, $K_0$, $\rho_0$) are known, the definition of $r_c = R/\tilde{R}$, (see (9)) using $\rho_c = \rho_s/\tilde{\rho}(\tilde{R})$ implies an implicit equation for $\tilde{R}$

$$ \left(\frac{R}{\tilde{R}}\right)^2 = \frac{K_0}{4\pi G \rho_0^p} \left(\frac{\rho_s}{\tilde{\rho}(\tilde{R})}\right)^{p-2}. \quad (14) $$

Similarly when ($M$, $R$, $K_0$, $\rho_0$) are known, we start from

$$ M = 4\pi \int_{0}^{R} \rho r^2 dr = 4\pi \rho_c r_c^3 \int_{0}^{\tilde{R}} \tilde{\rho} \tilde{r}^2 d\tilde{r}, \quad (15) $$

which can also be written, using (11), as

$$ \begin{aligned} M &= -4\pi\rho_c r_c^3 \int_0^{\tilde{R}} \frac{d}{d\tilde{r}} \left( \tilde{r}^2 \tilde{\rho}^{p-2} \frac{d\tilde{\rho}}{d\tilde{r}} \right) d\tilde{r} \\ &= -\frac{4\pi\rho_c R^3}{\tilde{R}} \tilde{\rho}^{p-2}(\tilde{R}) \frac{d\tilde{\rho}}{d\tilde{r}}(\tilde{R}). \end{aligned} \quad (16) $$

To get again an implicit equation for $\tilde{R}$, we can eliminate $\rho_c$ from this equation using (9) and (12), to get

$$ M = - \left( \frac{(4\pi)^{p-1} G \rho_0^p R^{3p-4}}{K_0 \tilde{R}^p} \right)^{\frac{1}{p-2}} \tilde{\rho}^{p-2} (\tilde{R}) \frac{d\tilde{\rho}}{d\tilde{r}} (\tilde{R}). \quad (17) $$

Finally when ($\rho_s$, $M$, $K_0$, $\rho_0$) are known, we eliminate $R$ from (17) using (14) which leads to a third implicit equation of $\tilde{R}$. In each of the three cases, when $\tilde{R}$ is determined, the central and surface densities or the mass are readily obtained.

In the remainder of this article, we take $K_0 = 130$ GPa which is typical incompressibility near the surface of the Earth, and $\rho_0 = 4000$ kg m$^{-3}$ (for planets made of silicates and metal, we thought it reasonable to choose a reference density somewhat larger than that of crustal rocks or very shallow mantle—see also Ricard et al. 2022—but this choice is not crucial). The length $r_0$ is therefore $r_0 = 3113$ km. At depth, this incompressibility increases continuously with the density according to (5). For example, a seismologic model of the Earth like PREM indicates that mantle density increases with depth by a factor of 1.7 (from $\approx 3200$ kg m$^{-3}$ to $\approx 5400$ kg m$^{-3}$) while incompressibility increases by a factor of 5.0 (to $\approx 650$ GPa). Using (5), this implies $p = \ln 5.0/\ln 1.7 = 3.03$ which again confirms that $p \approx 3$ is an appropriate choice. The Earth's density has also several large discontinuities with depth (in the transition zone, the core–mantle and inner–outer core boundaries) while the incompressibility only exhibits minor discontinuities with depth. This is another argument that led us to prefer, for our continuous model, a choice of numerical values that corresponds to the observation of the relatively continuous behavior of incompressibility.

Notice that $K_0$ and $\rho_0$ appear in the equations only by the ratio $K_0/\rho_0^p$ in the definition (9) of $r_c$. Compositionally denser (resp. lighter) materials have often a larger (resp. smaller) incompressibility, for example, $K_0/\rho_0^3$ for silicates and water
<table><caption>**Table 1**
Mass, Radius, and Average Densities of Various Planets</caption><thead><tr><th></th><th><i>M</i><br>(10<sup>22</sup> kg)</th><th><i>p̄</i><br>(kg m<sup>-3</sup>)</th><th><i>R</i><br>(km)</th></tr></thead><tbody><tr><td>CoRoT-7b</td><td>3620</td><td>9360</td><td>9735</td></tr><tr><td>Earth</td><td>597.2</td><td>5515</td><td>6371</td></tr><tr><td>Mars</td><td>64.17</td><td>3933</td><td>3389</td></tr><tr><td>Mercury</td><td>33.01</td><td>5427</td><td>2439</td></tr><tr><td>Ganymede</td><td>14.82</td><td>1940</td><td>2631</td></tr><tr><td>Moon</td><td>7.342</td><td>3344</td><td>1737</td></tr></tbody></table>

are similar (2.03 and 2.1 Pa m kg³, using $K_0 = 2.1$ GPa and $\rho_0 = 1000$ kg m⁻³ for water). The exponent *p* appropriate for water or ices is also in the range of those appropriate for silicates, close to 4 (Fei et al. 1993). Our model can therefore be extended to a larger variety of compositions than silicate planets.

## 3.2. Density Profile

We first apply our model to various planets for which *M* and *R* are known, see Table 1. We consider five planets or satellites of the solar system (Moon, Ganymede, Mars, Mercury, and Earth) and also examine the case for the exoplanet CoRot-7b using the mass and radius determinations of John et al. (2022) and Barros et al. (2014). The uncertainties suggested in these two papers may be underestimated as other articles have proposed values outside the corresponding confidence intervals. Among the planets selected, we include Ganymede whose composition probably has equal parts of rocky material and water, liquid, or ices.

Figure 2 shows the density profiles calculated for the Earth for *p* = 2, 3, and 4 (black, red, blue, solid lines). The length $r_c$ and density values obtained for various *p* are listed in Table 2. The PREM density profile (black, long-dashed) is shown for comparison. Of course, the large discontinuity between the core and mantle is not reproduced by our simple model. The prediction with *p* = 2 gives an overall better fit to PREM density, although the gradients of curves with *p* = 3 or 4 give a better fit to the gradient of PREM, at least in the mantle, which is a better fit to incompressibility.

The choice of *p* = 2 is convenient for checking the accuracy of numerical solutions because the solution is analytical with $\rho = \rho_c \text{sinc}(r/r_c)$, and $r_c$ is independent of $\rho_c$ ($r_c = r_0 = 3113$ km, see (9)). However, in addition to the fact that *p* = 2 is smaller than experimentally observed, density positivity requires $R \le \pi r_c = 10,408$ km. Clearly, *p* = 2 is not appropriate for large exoplanets like CoRot-7b. For *p* > 2, the dependence of $r_c$ on $\rho_c$ maintains the positivity of the surface density even for very large planets.

Figure 2 clearly shows that the question of planet differentiation requires a more complex model. However, the overall behavior of a compressible planet depends on the volumes subjected to compression, and the volume of Earth's core represents only 16% of the Earth's volume, so the discrepancy between the actual density and the modeled density near the center is not very significant for the purposes of this paper. Furthermore, if the presence of a metallic core is proven, the construction of a two-layer Lane–Emden model is possible. For example, in Figure 2 (blue dashed line), we add a two-layer model of the Earth that is in agreement with the

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Predicted densities of a Lane–Emden planet with the mass and radius of the Earth, with *p* = 2, 3, or 4, compared to the PREM density model (dotted line).</figcaption></figure>

<table><caption>**Table 2**
Quantities Computed for Various Planets and Exponents *p* Calculated from Our Model of Lane–Emden Planets According to the EoS (4)</caption><thead><tr><th><i>p</i></th><th><i>r</i><sub>c</sub><br>(km)</th><th><i>ρ</i><sub>s</sub><br>(kg m<sup>-3</sup>)</th><th><i>ρ</i><sub>c</sub><br>(kg m<sup>-3</sup>)</th><th><i>M̃</i> = <i>M</i>/<i>M</i><sub>s</sub><br>= <i>p̄</i>/<i>ρ</i><sub>s</sub></th></tr></thead><tbody><tr><td>CoRotT-7b</td><td></td><td></td><td></td><td></td></tr><tr><td>3</td><td>5761</td><td>6009</td><td>13695</td><td>1.559</td></tr><tr><td>4</td><td>8558</td><td>8092</td><td>10995</td><td>1.158</td></tr><tr><td>Earth</td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td>3113</td><td>3747</td><td>8625</td><td>1.471</td></tr><tr><td>3</td><td>4210</td><td>4191</td><td>7315</td><td>1.316</td></tr><tr><td>4</td><td>5209</td><td>4549</td><td>6693</td><td>1.212</td></tr><tr><td>Mars</td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td>3113</td><td>3612</td><td>4438</td><td>1.089</td></tr><tr><td>3</td><td>3274</td><td>3598</td><td>4422</td><td>1.093</td></tr><tr><td>4</td><td>3430</td><td>3583</td><td>4407</td><td>1.098</td></tr><tr><td>Mercury</td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td>3113</td><td>5201</td><td>5774</td><td>1.043</td></tr><tr><td>3</td><td>3709</td><td>5259</td><td>5675</td><td>1.032</td></tr><tr><td>4</td><td>4364</td><td>5303</td><td>5607</td><td>1.023</td></tr><tr><td>Ganymede</td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td>3113</td><td>1848</td><td>2088</td><td>1.051</td></tr><tr><td>3</td><td>2329</td><td>1737</td><td>2238</td><td>1.118</td></tr><tr><td>4</td><td>1953</td><td>1433</td><td>2509</td><td>1.355</td></tr><tr><td>Moon</td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td>3113</td><td>3276</td><td>3452</td><td>1.021</td></tr><tr><td>3</td><td>2900</td><td>3261</td><td>3471</td><td>1.026</td></tr><tr><td>4</td><td>2719</td><td>3243</td><td>3494</td><td>1.032</td></tr></tbody></table>

Note. The masses and radii of the planets are given in Table 1.

Earth's mass, inertia, and surface density. We use *p* = 2, as the solution is analytical, but a numerical solution for other values of *p* would not be difficult. As layered models can only be proposed for a very limited number of objects in the solar system, we will only consider homogeneous models in the following.

For the different planets considered and for various *p*, we report the length $r_c$ and the planet densities at the surface and at the center in Table 2. The last column of Table 2 gives the values of the quantity

$$ \tilde{M} = \frac{M}{M_s} = \frac{\bar{p}}{\rho_s}, \qquad (18) $$
<figure><img src="image_6.png" /><figcaption>**Figure 3.** We show the predicted radii for planets of mass *M* and surface densities *ρ*<sub>s</sub>, calculated with *p* = 3 (panel (a)) and *p* = 4 (panel (b)). The surface densities are 2000 kg m<sup>-3</sup> (red dashed line), or 6000 kg m<sup>-3</sup> (red solid line). In panels (c) and (d), we show the quantity *M*̃ which characterizes the importance of compressibility for each planet. In all panels, the observations for the Moon, Ganymede, Mercury, Mars, Earth, and CoRoT-7b are shown with black circles (see values in Tables 1 and 2). We also show two limiting cases: the case of an incompressible planet of uniform density 4000 kg m<sup>-3</sup> (black dashed in panels (a) and (b), and *M*̃ = 1 in panels (c) and (d)) and the asymptotic behavior for very large planets (magenta dashed). The prediction by Seager et al. (2007) is shown by blue dashed lines.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

which is the mass of the planet normalized by the mass *M*<sub>s</sub> of a
planet having the same radius but a homogeneous density ρ<sub>s</sub>

$$
M_s = \frac{4\pi}{3} \rho_s R^3, \qquad (19)
$$

and $\bar{\rho}$ is the planet's average density. The value of $\widetilde{M}$ therefore quantifies the extent to which density is affected by compressibility. As expected, the values of $\widetilde{M}$ increase with the planet's radius.

3.3. *Mass–Radius Relations*

We can also estimate the effect of compressibility according to
our model for any planet whose mass is known. Since for a given
mass *M*, the radius of the planet depends on its composition, we
make two assumptions, one with a rather low surface density
ρ<sub>s</sub> = 2000 kg m<sup>-3</sup>, the other with a rather large surface density
ρ<sub>s</sub> = 6000 kg m<sup>-3</sup>. In Figure 3, we show the radii *R* of a planet with

mass *M* and surface density *ρ*<sub>s</sub>, in the cases *p* = 3 (panel (a)) and *p* = 4 (panel (b)). They are very similar, with radii slightly larger in the second case, especially for masses ≥2 × 10<sup>25</sup> kg. The planet's radius should lie in the shaded area, between the two red curves: closer to the dashed red line for a planet with a light composition, and closer to the solid red line for a planet with a denser composition. For very massive planets, the radius becomes independent of the surface density and the width of the shaded area decreases. The radii of the planets in Table 1 are also shown in both panels. They lie within or close to the shaded area, closer to the solid line for Mercury and its large core, or for CoRoT-7Bb which probably also has a large core (Wagner et al. 2012), closer to the dashed line for Ganymede rich in water. As indicated in Table 2, for *p* = 3, Ganymede and CoRoT-7b have surface densities slightly outside the 2000–6000 kg m<sup>-3</sup> interval of the shaded areas. To make the effect of compressibility more visible, in Figures 3(c) and (d), we report the values of *M*̃ = *M*/*M*<sub>s</sub> as a function of the planet's mass.
The behavior of the predicted mass–radius relationship can be understood with the help of two limiting cases. The gravity of small planets is too low for compressibility to be an important factor. A black dashed line in Figure 3(a) and (b) shows the radius of an incompressible planet (i.e., $M \propto R^3$), with homogeneous density $\rho_s = 4000 \text{ kg m}^{-3}$. It lies above the curve corresponding to a compressible planet: compressibility decreases the radius for a given planet mass. An incompressible planet corresponds to $\tilde{M} = 1$ in Figure 3(c) and (d). Another limiting case is for a very large planet. From the mass expression (15) and the expression of $\rho_c$ (17), we get that for $p \neq 2$ (for $p=2$ the radius becomes independent of $M$ at large $M$)

$$M = 4\pi\rho_0 R^3 \left(\frac{R}{r_0}\right)^{\frac{2}{p-2}} S, \qquad (20)$$

where

$$S = \frac{1}{\tilde{R}^{\frac{3p-4}{p-2}}} \int_{0}^{\tilde{R}} \tilde{\rho} \tilde{r}^2 d\tilde{r}. \qquad (21)$$

When $M \to +\infty$, $\rho_s/\rho_c \to 0$ and the normalized radius $\tilde{R}$ is bounded by its value when $\tilde{\rho} = 0$ (equal to $\pi$, 1.94 or 1.53 for $p=2, 3$ and 4, see Figure 1). The asymptotic value of $S$ can be numerically calculated: $S \to 0.049$ for $p=3$ and $S \to 0.144$ for $p=4$. This result is given here with our notations, but the asymptotic relationship between $M$ and $R$ ($M \propto R^5$ and $M \propto R^4$ for $p=3$ and 4) is a well-known result (see Chandrasekhar 1958, p. 98). This relationship and that between $\tilde{M}$ and $M$ with $p=3$ are represented in Figures 3(a)–(b) by magenta dotted-dashed lines. Planets with a mass less than $10^{23}$ kg are basically incompressible, while planets with a mass greater than $10^{26}$ kg have a radius that follows the asymptotic regime. The Earth has a mass for which none of the limiting cases apply.

Several authors have studied the mass–radius relationship for large condensed planets. They have sometimes considered both a more detailed approach (with planets including silicate mantles, metallic cores, or oceans) and a more sophisticated EoS (e.g., a third-order Birch–Murnaghan equation). Their predictions are however closely similar to ours (see, e.g., Valencia et al. 2006; Sotin et al. 2007; Wagner et al. 2012) with $R \propto M^a$ with an exponent $a$ decreasing with the planet’s mass, from $1/3$ (incompressible case for small planets) to an asymptotic value that we found equal to $(p-2)/(3p-4)$, between 0.20 and 0.25. For $M_{\oplus} \le M \le 10M_{\oplus}$ we obtain a value of $\approx 0.27$, in agreement with previous findings.

Seager et al. (2007) use, like in our paper, a Lane–Emden equation with various possible compositional stratifications. They propose a generic mass–radius relation valid for more or less all compositions, up to a few terrestrial radii, on the form (with our notations)

$$\ln \tilde{M} = \left(\frac{M}{M_1}\right)^{k_3} \quad \text{or} \quad R^3 = \frac{3M}{4\pi\rho_s} \exp \left[ -\left(\frac{M}{M_1}\right)^{k_3} \right] \qquad (22)$$

where they predict an exponent $k_3 \approx 0.4$. The Equation (23) of Seager et al. (2007) is written in a more complex form as a relation between $(1/3)\log_{10}(M/\rho_s)$ and $\log_{10}R$, including an adjusted constant $k_1 \approx -0.209$ which should be logically $(1/3)\log_{10}(4\pi/3) = -0.207$ to insure a sound behavior when

$M \to 0$. The constant $M_1$ is $m_1(3k_2 \ln(10))^{-1/k_3}$ using their notations. We choose $M_1 = 1.6 \times 10^{26}$ kg which corresponds to $m_1 = 6.5M_{\oplus}$ which is in the range of values proposed in Seager et al. (2007), 4.34–10.55$M_{\oplus}$. Their predictions are also plotted in Figure 3, panels (a) and (b) (blue lines). As their models consider stratified planets in which each building layer verifies a polytropic equation with $p \approx 3–4$, their prediction remains in agreement with our simpler model for the same range of exponent $p$. However, for very large planet mass, their mass–radius parameterization diverges from the expected asymptotic behavior.

Our approach leads to a mass–radius relationships in perfect agreement with previous, more complex attempts (Valencia et al. 2006; Seager et al. 2007; Sotin et al. 2007). This suggests that the exponent parameter of the Murnaghan EoS (or Lame Emdeen EoS) controls the relationship more than the compositional details. The advantage of our simple model is that it can easily be perturbed analytically when certain conditions, such as the planet’s internal temperature or rate of rotation, change. This is the subject of the following paragraphs.

### 4. Thermal Contraction

What happens now when the temperature of a planet changes while remaining close to an adiabatic state? Although these temperature evolutions are small, how is the Earth’s radius affected by its cooling rate of 50–100 K Gyr$^{-1}$ (Herzberg et al. 2010)? The density profile that we compute does not explicitly include the temperature at each depth. However, from the EoS, the adiabatic temperature profile in the planet can be easily derived when the surface density and therefore the surface temperature are chosen (e.g., Ricard & Alboussière 2023). In this section, we perturb the surface temperature and calculate the resulting radius change. This does not change the solution without dimension $\tilde{\rho}$ of Equation (11) which is independent of any parameter. However, the solution with physical dimensions is affected by changes (denoted with $\delta$) in the quantities $\rho_s$, $\rho_c$, $r_c$, and $u$. Since the mass of the planet does not change, by perturbing (15) we obtain

$$4\pi\delta\rho_c r_c^3 \int_0^{\tilde{R}} \tilde{\rho}\tilde{r}^2 d\tilde{r} + 12\pi\rho_c r_c^2 \delta r_c \int_0^{\tilde{R}} \tilde{\rho}\tilde{r}^2 d\tilde{r} + 4\pi\rho_c r_c^3 \tilde{\rho}(\tilde{r})\tilde{R}^2 \delta\tilde{R} = 0 \qquad (23)$$

which can be reset using (12)–(15)–(19) as

$$\frac{\delta\rho_c}{\rho_c}M + 3\frac{\delta r_c}{r_c}M + 3\left(\frac{\delta R}{R} - \frac{\delta r_c}{r_c}\right)M_s = 0. \qquad (24)$$

The perturbation of (9) leads to

$$\frac{\delta r_c}{r_c} = \frac{p-2}{2} \frac{\delta \rho_c}{\rho_c}. \qquad (25)$$

If the temperature changes while the pressure at the planetary surface is unchanged, the EoS implies that the surface density variation is

$$\frac{\delta \rho_s}{\rho_s} = -\alpha_s \delta T_s, \qquad (26)$$

where $\alpha_s$ is the thermal expansion at the surface and $T_s$ the adiabatic temperature extrapolated to the surface (sometimes
called the foot of the adiabat). The definition of $\rho_c$ closes the system. Indeed, the perturbation of $\rho_s = \rho_c \tilde{p}(\tilde{R})$ leads to

$$\delta\rho_s = \tilde{p}(\tilde{R})\delta\rho_c + \rho_c \frac{d\tilde{p}}{d\tilde{r}}(\tilde{R})\delta\tilde{R} \quad (27)$$

which can be reset as

$$\frac{\delta\rho_s}{\rho_s} = \frac{\delta\rho_c}{\rho_c} + \frac{\tilde{R}}{\tilde{p}(\tilde{R})} \frac{d\tilde{p}}{d\tilde{r}}(\tilde{R}) \left( \frac{\delta R}{R} - \frac{\delta r_c}{r_c} \right). \quad (28)$$

However from (16) and the definition of gravity $g = GM/R^2$ we obtain

$$\begin{aligned} \frac{\tilde{R}}{\tilde{p}(\tilde{R})} \frac{d\tilde{p}}{d\tilde{r}}(\tilde{R}) &= -\frac{\tilde{M}\tilde{R}^2}{3\tilde{p}^{p-2}(\tilde{R})} \\ &= -\frac{\rho_s g R}{K_0} \left(\frac{\rho_0}{\rho_s}\right)^p = -\frac{\rho_s g R}{K_s}, \end{aligned} \quad (29)$$

where we introduce the surface incompressibility $K_s = K_0(\rho_s/\rho_0)^p$, leading to

$$\frac{\delta\rho_s}{\rho_s} = \frac{\delta\rho_c}{\rho_c} - \frac{\rho_s g R}{K_s} \left( \frac{\delta R}{R} - \frac{\delta r_c}{r_c} \right). \quad (30)$$

Finally using (24), (25), (26), and (30), we obtain:

$$\frac{1}{3}\alpha_s\delta T_s = \left( 1 + \frac{\rho_s g R}{3K_s} \frac{\tilde{p}\tilde{M}}{\tilde{p}\tilde{M} - 1} - \frac{\tilde{p}(\tilde{M} - 1)}{\tilde{p}\tilde{M} - 1} \right) \frac{\delta R}{R}, \quad (31)$$

with

$$\tilde{p} = \frac{p - 4/3}{p - 2}. \quad (32)$$

The parameter $\tilde{p}$ is therefore $+\infty$, $5/3$ or $4/3$ for $p=2, 3$, or 4 (Equation (31) remains valid for $\tilde{p} \to +\infty$ when $p \to 2$). In the case of an incompressible fluid, when $K_s \to \infty$ and $\tilde{M} \to 1$, are radius and temperature simply related by thermal expansivity only. Surprisingly, this is also the case for $\tilde{p} = 0$ when $p = 4/3$.

The nondimensional quantity $\rho_s g R / K_s$ can also be expressed as the ratio of two quantities well known to those working in compressible convection, namely

$$\frac{\rho_s g R}{K_s} = \frac{D}{\Gamma}, \quad (33)$$

where $D$ is the dissipation number and $\Gamma$ the Grüneisen parameter

$$D = \frac{\alpha_s g R}{C}, \quad \text{and} \quad \Gamma = \frac{\alpha_s K_s}{\rho_s C}. \quad (34)$$

In a vigorously convecting planet, $D$ controls the slope of the adiabatic temperature and $D/\Gamma$ the slope of the adiabatic density (e.g., Ricard & Alboussière 2023).

We define as the effective thermal expansion the quantity

$$\alpha_e = \alpha_s \left( 1 + \frac{\rho_s g R}{3K_s} \frac{\tilde{p}\tilde{M}}{\tilde{p}\tilde{M} - 1} - \frac{\tilde{p}(\tilde{M} - 1)}{\tilde{p}\tilde{M} - 1} \right)^{-1}. \quad (35)$$

The effective thermal expansion is also

$$\alpha_e = \alpha_s \frac{\tilde{p}\tilde{M} - 1}{\tilde{p} - 1 + \frac{4\pi\rho_s^2 GR^2}{9K_s}\tilde{p}\tilde{M}^2}, \quad (36)$$

where we use $g = (4\pi/3)G\rho_s R\tilde{M}$. When $\tilde{M} \to +\infty$, $\alpha_e \to 0$ and therefore the effective thermal expansivity is close to 0 for large planets. The effective thermal expansion is also lower than $\alpha_e$ when $\tilde{M} \approx 1$ if $\tilde{p} > 1$. In all the cases that we have considered, i.e., $p \ge 2$ which implies $\tilde{p} > 1$, the compressibility decreases the thermal contraction, $\alpha_e \le \alpha_s$. However for $0 \le p < 2$ which implies $\tilde{p} \le 2/3$, $\alpha_e$ becomes larger than $\alpha_s$ when $\tilde{M} \approx 1$ and compressibility enhances the contraction of a small cooling planet. This is why Jaupart et al. (2015) who assumed a constant incompressibility concluded that compressibility enhances the contraction of the Earth. Their expression is identical to ours when $p=0$ and the compressibility is low. However, using a Murnaghan EoS with $p < 2$ (a polytropic index $n \ge 1$) is inappropriate for a condensed planet, and with reasonable exponents $p$, compressibility always decreases the thermal contraction.

The effective thermal expansion is shown in Figure 4. We only consider the case $p=3$ and three possible surface densities: $\rho_s = 2000$, $\rho_s = 6000 \text{ kg m}^{-3}$ (same cases as in Figure 3), and $\rho_s = 4191 \text{ kg m}^{-3}$ (the value found for the Earth, see Table 1). For the various planets we used the values of $\tilde{M}$ and $\rho_s$ for $p=3$, from Table 1.

Already for the Earth, the effective thermal expansivity is significantly reduced compared to its surface value: this is due both to the reduction of $\alpha$ with depth and to the trade-off between pressure-temperature-density and gravity. For masses above $10^{26}$ kg the planet's radius becomes insensitive to temperature: the density profile is controlled solely by incompressibility.

Note that the expression (31) relates the radius to the adiabatic temperature at the surface $T_s$. Using the Eos (4), we could easily derive that the adiabatic temperature $T$ and the density are related by (see, e.g., Ricard & Alboussière 2023)

$$T = T_s \exp \left[ \Gamma \left( \frac{\rho_0}{\rho_s} - \frac{\rho_0}{\rho} \right) \right]. \quad (37)$$

The adiabatic temperature increases with depth but this increase remains moderate as bounded by $\exp(\Gamma) \lesssim 3$ (obtained when $\rho \gg \rho_0$, since $\rho_s \approx \rho_0$ and $\Gamma \approx 1$). The average temperature $\tilde{T}$ is therefore comparable to $T_s$ in small planets and less than $\approx 3$ times larger in massive planets. The relative variations of $\delta T_s / T_s$ and $\delta \tilde{T} / \tilde{T}$ are comparable. If instead of relating the radius changes to the surface temperature, we rather use the average temperature of the planet, the effective thermal expansion must be further multiplied by $T_s / \tilde{T}$, $\approx 1$ for small planets, $\lesssim 1/3$ for large planets. In this case, $\alpha_e$ decreases even faster with the size of the planet.

## 5. Change of Rotation Velocity

In this section, we discuss how a change in rotation rate affects the planet's radius. To account for planetary rotation, a centrifugal potential must be added to Poisson's Equation (7).
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Effective thermal expansivity for *p* = 3 of a planet as a function of its mass. We consider three possible surface densities for the condensed planets: that appropriate for the Earth (solid blue), and for a compositionally denser (solid red) or lighter (dashed red) planet. Heavy planets like super-Earth, have a very low effective expansivity as their density is solely related to compressibility.</figcaption></figure>

With this term, Equation (7), verified by density becomes

$$ \frac{1}{r^2} \frac{d}{dr} \left( K_0 r^2 \frac{\rho^{p-2}}{\rho_0^p} \frac{d\rho}{dr} \right) + 4\pi G \rho - 2\Omega^2 = 0. \quad (38) $$

Introducing the previously defined $\tilde{\rho}$ et $\tilde{r}$, the equation to solve with rotation, is

$$ \frac{1}{\tilde{r}^2} \frac{d}{d\tilde{r}} \left( \tilde{r}^2 \tilde{\rho}^{p-2} \frac{d\tilde{\rho}}{d\tilde{r}} \right) + \tilde{\rho} = \tilde{\Omega}^2 \tilde{\rho} (\tilde{R}) \quad (39) $$

where

$$ \tilde{\Omega}^2 = \frac{\Omega^2}{2\pi G \rho_s}. \quad (40) $$

In this case, it is not possible to find a generic solution, as previously with $\tilde{\rho}$, so that the physical solution with dimensions depends only on scaling parameters like $R/r_c$ and $\rho_c/\rho_s$. Here the solution will necessarily depend on a third quantity, $\tilde{\Omega}^2$. It therefore seems that, in general, only a numerical solution can be sought if we want to quantify the relationship between the planetary radius and its angular rotation.

Only for $p=2$ the solution can be found analytically and is

$$ \tilde{\rho} = \tilde{\rho}(\tilde{R}) \left( (1 - \tilde{\Omega}^2) \frac{\mathrm{sinc} \tilde{r}}{\mathrm{sinc} \tilde{R}} + \tilde{\Omega}^2 \right). \quad (41) $$

Using this expression, the mass of the planet is

$$ \begin{aligned} M &= 4\pi \rho_c r_c^3 \int_0^{\tilde{R}} \tilde{\rho} \tilde{r}^2 d\tilde{r} \\ &= M_s (3(1 - \tilde{\Omega}^2) F(\tilde{R}) + \tilde{\Omega}^2), \end{aligned} \quad (42) $$

where

$$ F(\tilde{R}) = \frac{\sin \tilde{R} - \tilde{R} \cos \tilde{R}}{\tilde{R}^2 \sin \tilde{R}}. \quad (43) $$

When the planet's rotation changes, the composition does not change and unlike in Section 4 the surface temperature is unaffected. So, by changing $\tilde{\Omega}$, the density will change at depth

but not the boundary conditions at the surface that maintain $\rho_s$ fixed.

Therefore, Equation (42) relates the rotation rate of a planet to its radius which appears both in $\tilde{R} = R/r_c$ and in $M_s = 4\pi\rho_s R^3/3$. We can now differentiate (42) taking into account that $\delta M = 0$, $\delta M_s = 3M_s\delta R/R$, $\delta F = (dF/d\tilde{R})\delta\tilde{R}$, $\delta\tilde{R} = \tilde{R}\delta R/R$, and get

$$ \begin{aligned} & 3\left( (1 - \tilde{\Omega}^2)(3F + \tilde{R}\frac{dF}{d\tilde{R}}) + \tilde{\Omega}^2 \right) \frac{\delta R}{R} \\ &= 2\tilde{\Omega}^2(3F - 1) \frac{\delta\tilde{\Omega}}{\tilde{\Omega}}. \end{aligned} \quad (44) $$

This expression can be simplified through a rather cumbersome algebra. First, the definition of $F$ leads to

$$ 3F + \tilde{R} \frac{dF}{d\tilde{R}} = 1 + \tilde{R}^2 F^2, \quad (45) $$

where $F$, extracted from the mass conservation (42) is

$$ F = \frac{\tilde{M} - \tilde{\Omega}^2}{3(1 - \tilde{\Omega}^2)} \quad (46) $$

with again $\tilde{M} = M/M_s$. Defining the parameter

$$ m = \frac{\Omega^2 R}{g}, \quad (47) $$

which is central in the planet's hydrostatic theory, allows us to write $\tilde{\Omega}^2$ as

$$ \tilde{\Omega}^2 = \frac{2}{3} m \tilde{M}. \quad (48) $$

Using (29) and all simplifications done, we obtain

$$ \frac{\delta R}{R} = \frac{4 m \tilde{M} (\tilde{M} - 1)}{9 \left( 1 - \frac{2}{3} m \tilde{M} + \frac{\rho_s g R}{3 K_s} \tilde{M} \left( 1 - \frac{2}{3} m \right)^2 \right) \Omega}. \quad (49) $$

The $\tilde{M} - 1$ term insures that the planetary radius is independent of the rotation rate when the planet is incompressible.

By comparison with the result of Saito (1974; see Equation (3)), our model predicts a Love number

$$ \begin{aligned} h_0^* &= \frac{2\tilde{M}(\tilde{M}-1)}{3 - 2m\tilde{M} + \frac{\rho_s g R}{K_s}\tilde{M}\left(1 - \frac{2}{3}m\right)^2} \\ &\approx \frac{2\tilde{M}(\tilde{M}-1)}{3 + \frac{\rho_s g R}{K_s}\tilde{M}} \end{aligned} \quad (50) $$

where the terms in $m$ can generally be omitted ($m=1/289$ for the Earth). The parameter $\rho_s g R/K_s = D/\Gamma$ that we discussed in Section 4 appears again.

We plot in Figure 5(a), the Love number $h_0^*$ as a function of planetary mass for a planet having three possible surface densities $\rho_s = 2000$, $\rho_s = 6000 \text{ kg m}^{-3}$ (same cases as in Figures 3 and 4), and $\rho_s = 3747 \text{ kg m}^{-3}$ (the value found for the Earth, see Table 1) and $m \ll 1$. When the mass is large, $h_0^*$ reaches an asymptotic value, but in this case $m = \Omega^2 R/g = 3\Omega^2/(4\pi G\bar{\rho})$ goes to zero and the radius becomes independent of the rotation rate (see (3)). This is more conspicuous in Figure 5(b), where we plot $h_0^*\rho_s/\bar{\rho}$ that
<figure><img src="image_9.png" /><figcaption>**Figure 5.** (a) Degree 0 rotational Love number as a function of the planetary mass. The curves depict the analytical solution for *p* = 2 (Equation (50)) for three values of the surface density. The filled circles correspond to the analytical solution for various planets, using the values of the surface density taken from Table 1 with *p* = 2 (we excluded CoRoT-7b, too massive for the case *p* = 2). The red and green curves are the numerical solution of the equations for *p* = 3 and 4, and the black square is Saito's value. (b) As the contraction of a planet during despinning contains the flattening term *m* which decreases as $1/\bar{\rho}$, we plot $h_0^* \rho_s / \bar{\rho}$ which accounts for all the terms depending on the planet's mass, as a function of the planet's mass. The prediction of Saito corresponds to a value of *p* slightly lower than 4.</figcaption></figure>

includes all the terms dependent on the mass. Small planets are incompressible and their radius is independent of the rotation rate, large planets have such a large density that $h_0^* \rho_0 / \bar{\rho}$ is very small, and the Earth and CoRoT-7b are in a mass range where their radius is most sensitive to their rotation rate.

The value we obtain for the Earth, $h_0^* = 0.23$ is twice as large as the value proposed by Saito (1974). Our analytical model, however, uses *p* = 2 which is too small a value. To confirm that this large $h_0^*$ is only due to an inappropriate choice of *p*, we take a brute-force approach and numerically solve (38) for *p* = 2, 3, or 4, and for a surface density of 3747, 4191, and 4549 kg m⁻², respectively (the values of surface densities listed in Table 1 for the Earth). The rotation rate is set successively to the rotation period of the Earth $\Omega_e$, then to $\Omega_e - \delta\Omega_e$, and to $\Omega_e + \delta\Omega_e$. Noting the corresponding radii $R^-$ and $R^+$, the Love number is approximated by (see Equation (3))

$$h_0^* = \frac{3}{4mR} \frac{R^+ - R^-}{\delta\Omega_e} \quad (51)$$

The results of these numerical experiments are also shown in Figure 5(a). The numerical estimate for *p* = 2, in perfect agreement with the analytical result, is not repeated. For *p* = 3, the numerical estimate for the Earth is closer to the value proposed by Saito (1974) and for *p* = 4, it becomes basically identical (even slightly lower). This confirms that a *p* = 2 model is too compressible and leads to a planet that is excessively sensitive to changes in temperature or rotation speed.

## 6. Discussion and Conclusion

The Lane–Emden models have been widely used in astrophysics. Their application to condensed planets has been less so although the EoS of Murnaghan’s EoS (Murnaghan 1937), often used in geophysics for silicate of metal planets, belongs to the family of polytropic EoS. An important difference between astrophysical models of stars or gas planets and telluric planets comes from a difference in the polytropic exponent *n*. Another difference is that the ratio between central and surface densities in telluric planets is never very large so that the surface boundary conditions remain crucial while for say, a giant gas planet, the surface density is zero and the density profile or the mass only depends on the central density $\rho_e$.

For the Earth, an accurate elastic model of compressibility is known, and even for several objects of the solar system, more realistic compressibility than Lane–Emden models can be proposed. The aim of this article is to propose a simple generic model for all planets whose properties are not precisely known. For any specific planet, a more detailed model could account for their composition, differentiation, and phase changes. At any rate, for the smallest objects (Moon, Mercury, and Mars), the assumption of total incompressibility makes little difference: the change in radius is mainly linked to thermal expansivity and is insensitive to rotation. For the Earth compressibility plays a minor but significant role (see Figures 4 and 5), for CoRoT-7B or larger super-Earth, it becomes a major effect.

In this paper, we have refrained from making numerical applications of radius changes for the planets that we have previously examined (Earth, Mars, Mercury, and the Moon); they are not very different from previous estimates. The radius decreases due to the cooling of the Earth (say 250 K in 3 Gyr with $\alpha_0 = 3 \times 10^{-5}$ K⁻¹) should be around 9 km using $\alpha_e / \alpha_s = 0.6$ in Figure 4. The Earth’s sidereal day was only 13 or 15 hr, during the Archean, 3.2 Gyr ago; see, e.g., Farhat et al. (2022) or Eulenfeld & Heubeck (2023), using the geological records of tidalites (Eriksson & Simpson 2001). This implies a further radius reduction of 1.4 km (as $\Omega$ varied significantly we integrate (3) and use $\delta R = (1/3)h_0^* mR(1 - (\Omega_A^2 / \Omega_e^2))$ where $\Omega_A$ is the Archean rotation rate). Of course, plate tectonics has erased all evidence of these contractions which only reach $\approx 3$ µm yr⁻¹; only on planets whose lithosphere has been frozen for billions of years can thermal contractions be observed.

Our main objective was to show that the important parameters controlling the changes in the radius are the dissipation number $D = \alpha_s gR/C$ and the Grüneisen parameter $\Gamma = \alpha_s K_s / (\rho_s C)$. The Grüneisen parameter varies little between 1 and 2 for most planet’s compositions. On the contrary, very
large dissipation numbers are specific to super-Earths since R and g increase with the planet's mass (Ricard & Alboussière 2023). From Figures 4 and 5, it is clear that Earth lies in something of a transition zone, between smaller objects affected by temperature, and larger objects which are mostly insensitive to temperature but affected by their angular rotation. For the large exoplanets that have been discovered, dissipation numbers larger than 10 are expected ($D \approx 0.6$ for the Earth) from the observed radius and masses (see, e.g., Otegi et al. 2020). In the range $10^{24}-10^{26}$ kg and for planets whose interiors are largely unknown, our approach can provide first-order estimates of density profiles and potential changes of the planetary radii through time.

## Acknowledgments

This research was founded by the French National Program of Planetology (PNP, CNRS-INSU, Proposal: Compressibility and Convection in Condensed Planets).

## ORCID iDs

Yanick Ricard https://orcid.org/0000-0002-0998-4670
Frédéric Chambat https://orcid.org/0000-0001-8967-077X

## References

* Anderson, O. L. 1979, JGR, 84, 3537

* Barros, S. C. C., Almenara, J. M., Deleuil, M., et al. 2014, A&A, 569, A74

* Bullen, K. E. 1975, The Earth's Density (London: Chapman and Hall)

* Byrne, P. K., Klimczak, C., Sengor, A. M. C., et al. 2014, NatGe, 7, 301

* Chambat, F., Ricard, Y., & Valette, B. 2010, GeoJI, 183, 727

* Chandrasekhar, S. 1958, Introduction to the Study of Stellar Structure (New York: Dover Publications), 509

* Dana, J. D. 1847, AmJSA, 3, 176

* Eriksson, K., & Simpson, E. 2001, Geo, 29, 1159

* Eulenfeld, T., & Heubeck, C. 2023, JGRE, 128, e2022JE007466

* Farhat, M., Auclair-Desrotour, P., Boue, G., & Laskar, J. 2022, A&A, 665, L1

* Fei, Y., Mao, H., & Hemley, R. 1993, JChPh, 99, 5369

* Hauck, S., Dombard, A., Phillips, R., & Solomon, S. 2004, E&PSL, 222, 713

* Herzberg, C., Condie, K., & Korenaga, J. 2010, E&PSL, 292, 79

* Horedt, G. P. 2004, Polytropes—Applications in Astrophysics and Related Fields, Vol. 306 (Dordrecht: Springer)

* Jaupart, C., Labrosse, S., Lucazeau, F., & Mareschal, J.-C. 2015, in Treatise on Geophysics, ed. G. Schubert (2nd ed.; Oxford: Elsevier), 223

* John, A. A., Collier Cameron, A., & Wilson, T. G. 2022, MNRAS, 515, 3975

* Melosh, H. J., & McKinnon, W. B. 1988, in Mercury, ed. F. Vilas, C. R. Chapman, & M. S. Matthews (Tucson, AZ: Univ. Arizona Press), 374

* Murnaghan, F. D. 1937, AmJM, 59, 235

* Nahm, A. L., & Schultz, R. A. 2011, Icar, 211, 389

* Otegi, J. F., Bouchy, F., & Helled, R. 2020, A&A, 634, A43

* Ricard, Y., & Alboussière, T. 2023, PEPI, 341, 107062

* Ricard, Y., Alboussière, T., Labrosse, S., Curbelo, J., & Dubuffet, F. 2022, GeoJI, 230, 932

* Saito, M. 1974, JPE, 22, 123

* Seager, S., Kuchner, M., Hier-Majumder, C. A., & Militzer, B. 2007, ApJ, 669, 1279

* Solomon, S. C., & Chaiken, J. 1976, LPSC, 3, 3229

* Sotin, C., Grasset, O., & Mocquet, A. 2007, Icar, 191, 337

* Stixrude, L., & Lithgow-Bertelloni, C. 2005, GeoJI, 162, 610

* Valencia, D., O'Connell, R., & Sasselov, D. 2006, Icar, 181, 545

* Wagner, F., Tosi, N., Sohl, F., Rauer, H., & Spohn, T. 2012, A&A, 541, A103

* Watters, T. R., Robinson, M. S., & Cook, A. C. 1998, Geo, 26, 991