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

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 ρ s gR/K s (ρ 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 Gruneïsen 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.


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 which 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 Ricard@ens-lyon.frFrederic.Chambat@ens-lyon.fr to a reduction in the planet's radius.Conversely, the absence of obvious signs of planet compression or extension has been used to put constrains on planet's thermal evolution.The amount of planet radius variation has been discussed in the case on Mercury (3 to 7 km of contraction, e.g.Byrne et al. 2014), Mars (0 to 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 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 exemple, crystallization of a liquid metallic core to form a denser solid inner arXiv:2404.14278v1[astro-ph.EP] 22 Apr 2024 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 δR of its radius R is often used, for example which takes into account the variations of thermal expansivity α and temperature T , with depth (Solomon & Chaiken 1976), or even where ᾱ and δT are uniform (Hauck et al. 2004).For the Earth, an estimate of the cooling can be obtained from petrologic observation of ancien 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 α = 3 × 10 −5 K −1 would imply a radius reduction of 10-20 km in 3 Gyr.These relations (1)-( 2) are in fact problematic.Choosing an appropriate thermal expansivity for a quantity that varies with time and depth is a 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 that 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 fonction of temperature by perturbing the planet's elasto-gravitational 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 Ω, Saito (1974) has applied this method of perturbing the elastogravitational equations to obtain where g being Earth's surface gravity and h * 0 is a degree0 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.

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 ρ.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, hightemperature experiments on silicates and metals, as well as with the radius properties of the Earth (Ricard et al. 2022;Ricard & Alboussière 2023): In this expression α 0 , K 0 , ρ 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 non-dimensional empirical exponent.This Eos leads to simple expressions for thermal expansivity and isothermal incompressibility, which are related only to density 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 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 litterature since Chandrasekhar (1939) where the focus is often in stars and gaseous planets, and the index n is large (n = +∞ for a perfect gas).Neutron stars are well modeled with polytropic index close to, but lower than 1.In a solid planet p ≈ 3 − 4 or n ≈ 1/3 − 1/2 (Stixrude & Lithgow-Bertelloni 2005).The exponent p is also ∂K T /∂P and therefore it quantifies the increase of incompressibility with pressure.

Density equation
If we consider a non-rotating 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) In this expression, the quantity K is defined by If the planet is convecting, its state is close to an adiabatic state with an isentropic incompressibilityK 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 + ΓαT where the Grüneisen parameter Γ = αK T /ρC V is of order 1 for solid silicates, liquid metal (and even for ideal gases), and decreases slightly with depth, while αT is a small quantity of order 10 −2 (in planets, T increases with depth but α decreases more, see equations ( 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 that we later note C.

Non dimensional variables
We now define a dimensionless density ρ = ρ/ρ c and a dimensionless radius r = r/r c where ρ 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: where the length r 0 is Using these variables, Equation ( 7) becomes This equation must be solved with ρ(0) = 1 and dρ/dr(0) = 0, until the outer non dimensional radius r = R, where The resolution of ( 11) is analytical for p = 2, in which case the solution is the first arc of a sinc function, and for p = ∞ 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 r are shown in Figure 1 for p = 2, 3 and 4.These curves are computed until ρ = 0 but the density with dimensions are only related to these curves for r ≤ R, i.e. for ρ(r) ≥ ρ s /ρ c where ρ s is the surface density.Notice that, following Chandrasekhar (1939) the density of stellar objects in astrophysical literature is often sought in the form ρ = θ n = θ 1 p−1 (Horedt 2004).In which case (11) becomes the well known Lane-Emden equation ∇2 θ + θ n = 0 (13) when the radius in the Laplace operator ∇ 2 is normalized by r c / √ 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 ρ than θ.

Input parameters
To obtain a solution with physical dimensions, various quantities need to be specified.We need to know the central density ρ 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 ρ 0 .Therefore, four quantities need to be specified, (ρ c , R, K 0 , ρ 0 ).Alternatively we can choose another set of four independent parameters from which (ρ c , R, K 0 , ρ 0 ) can be derived.For example, we will calculate the mass of planets when (ρ s , R, K 0 , ρ 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 , ρ 0 ) and calculate the appropriate values of the central and surface densities.Finally, in the paragraph where we discuss the effects of rotation, (ρ s , M, K 0 , ρ 0 ) will be fixed values, and R will vary with the rate of rotation.When (ρ s , R, K 0 , ρ 0 ) are known, the definition of r c = R/ R, (see ( 9)) using Similarly when (M, R, K 0 , ρ 0 ) are known, we start from which can also be written, using (11), as To get again an implicit equation for R, we can eliminate ρ c from this equation using ( 9) and ( 12), to get Finally when (ρ s , M, K 0 , ρ 0 ) are known, we eliminate R from ( 17) using ( 14) which leads to a third implicit equation of R. In each of the three cases, when R is determined, the central and surface densities or the mass are readily obtained.
In the remainder of this article, we take K 0 = 130GPa which is a typical incompressibility near the surface of the Earth and ρ 0 = 4000 kg m −3 (for planets made of silicates and metal, we thought 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 exemple, a seismologic model of the Earth like PREM indicates that mantle density increases with depth by a factor of 1.7 (from ≈ 3200 kg m −3 to ≈ 5400 kg m −3 ) while incompressibility increases by a factor 5.0 (to ≈ 650 GPa).Using (5), this implies p = ln 5.0/ ln 1.7 = 3.03 which again confirms that p ≈ 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 ρ 0 appear in the equations only by the ratio K 0 /ρ p 0 in the definition (9) of r c .Compositionally denser (resp.lighter) materials have often a larger (resp.smaller) incompressibility, for example, K 0 /ρ 3 0 for silicates and water are similar (2.03 and 2.1 Pa m kg 3 , using K 0 = 2.1 GPa and ρ 0 = 1000 kg m −3 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.

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 (dotted line) is shown for comparison.Of course, the large dis- continuity between 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, that 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 ρ = ρ c sinc(r/r c ) (sinc(x) is the function (sin x)/x), and r c is independent of ρ 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 ≤ π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 ρ 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 of 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 layers 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 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 model 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 which is the mass of the planet normalized by the mass M s of a planet having the same radius but a homogeneous density ρ s and ρ is the planet's average density.The value of M therefore quantifies the extent to which density is affected by compressibility.As expected, the values of M increase with the planet radius.

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 ρ s = 2000 kg m −3 , the other with a rather large surface density ρ s = 6000 kg m −3 .In Figure 3, we show the radii R of a planet with mass M and surface density ρ s , 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 25 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 light composition, and closer to the solid red line for a planet with a denser composition.For very massives planets, the radius becomes independent of the surface density and the width of the shaded area decreases.The radii of the planets of 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 2000-6000 kg m −3 interval of the shaded areas.To make the effect of compressibility more visible, in Figure 3c and 3d, we report the values of M = M/M s 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 thin black dashed line in Figure 3a shows the radius of an incompressible planet (i.e.M ∝ R 3 ), with homogeneous density ρ s = 4000 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 M = 1 in Figure 3b.Another limiting case is for a very large planet.From the mass expression (15) and the expression of ρ c (17), we get that for p ̸ = 2 (for p = 2 the radius becomes independent of M at large M ) where When M → +∞, ρ s /ρ c → 0 and the normalized radius R is bounded by its value when ρ = 0 (equal to π, 1.94 or 1.53 for p = 2, 3 and 4, see Figure 1).The asymptotic value of S can be numerically calculated: S → 0.049 for p = 3 and S → 0.144 for p = 4.This result is given here with our notations, but the asymptotic relationship between M and R (M ∝ R 5 and M ∝ R 4 for p = 3 and 4) is a well-known result (see Chandrasekhar 1939, p. 98).This relationship and that between M and M with p = 3 are represented in Figures 3a-b by thin magenta dotted-dashed lines.Planets with mass less than 10 23 kg are basically incompressible, while planets with 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).They predictions are however closely similar to ours (see e.g., Valencia et al. 2006;Sotin et al. 2007;Wagner et al. 2012) with R ∝ 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 .25.For M ⊕ ≤ M ≤ 10M ⊕ we obtain a value of ≈ 0.27, in agreement with previous findings.Seager et al. (2007) uses, like in our paper, a Lane-Emden equation with various possible compositionnal 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)  1 and 2).We also show two limiting cases: the case of an incompressible planet of uniform density 4000 kg m −3 (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) are shown by blue dashed lines.
where they predict an exponent k 3 ≈ 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/ρ s ) and log 10 R, including an adjusted constant k 1 ≈ −0.209 which should be logically (1/3) log 10 (4π/3) = −0.207 to insure a sound behavior when M → 0. The constant M 1 is m 1 (3k 2 ln(10)) −1/k3 using their notations.We choose M 1 = 1.6 × 10 26 kg which corresponds to m 1 = 6.5M ⊕ which is in the range of values proposed in Seager et al. (2007), 4.34 − 10.55M ⊕ .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 ≈ 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 parametrization 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;Sotin et al. 2007;Seager 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.

THERMAL CONTRACTION
K s → ∞ and M → 1, are radius and temperature simply related by thermal expansivity only.Surprisingly, this is also the case for p = 0 when p = 4/3.
The non-dimensional quantity ρ s gR/K s can also be expressed as the ratio of two quantities well-known to those working in compressible convection, namely where D is the dissipation number and Γ the Grüneisen parameter In a vigorously convecting planet, D controls the slope of the adiabatic temperature and D/Γ the slope of the adiabatic density (e.g., Ricard & Alboussière 2023).
We define as the effective thermal expansion the quantity The effective thermal expansion is also where we use g = (4π/3)Gρ s R M .When M → +∞, α e → 0 and therefore the effective thermal expansivity is close to 0 for large planets.The effective thermal expansion is also lower than α e when M ≈ 1 if p > 1.
In all the cases that we have considered, i.e., p ≥ 2 which implies p > 1, the compressibility decreases the thermal contraction, α e ≤ α s .However for 0 ≤ p < 2 which implies p ≤ 2/3, α e becomes larger than α s when M ≈ 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 ≥ 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: ρ s = 2000, ρ s = 6000 kg m −3 (same cases as in Figure 3) and ρ s = 4191 kg m −3 (the value found for the Earth, see Table 1).For the various planets we used the values of M and ρ 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 α 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) The adiabatic temperature increases with depth but this increase remains moderate as bounded by exp(Γ) ≲ 3 (obtained when ρ ≫ ρ 0 , since ρ s ≈ ρ 0 and Γ ≈ 1).
The average temperature T is therefore comparable to T s in small planets and less than ≈ 3 times larger in massive planets.The relative variations of δT s /T s and δ T / 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 / T , ≈ 1 for small planets, ≲ 1/3 for large planets.In this case, α e decreases even faster with the size of the planet.

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).With this term, equation ( 7), verified by density becomes Introducing the previously defined ρ et r, the equation to solve with rotation, is where In this case, it is not possible to find a generic solution, as previously with ρ, so that the physical solution with dimensions depends only on scaling parameters like R/r c and ρ c /ρ s .Here the solution will necessarily depend on a third quantity, Ω 2 .It therefore seems that, in general, only a numerical solution can be sought if we want to 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 compositionnally 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.
quantify the relationship between the planetary radius and its angular rotation.
Only for p = 2 the solution can be found analytically and is Using this expression, the mass of the planet is where When the planet's rotation changes, the composition does not change and unlike in the previous section the surface temperature is unaffected.So, by changing Ω, the density will change at depth but not the boundary conditions at the surface that maintain ρ s fixed.Therefore, equation ( 42) relates the rotation rate of a planet to its radius which appears both in R = R/r c and in M s = 4πρ s R 3 /3.We can now differentiate (42) taking into account that δM = 0, δM s = 3M s δR/R, δF = (dF/d R)δ R, δ R = RδR/R, and get This expression can be simplified through a rather cumbersome algebra.First the definition of F leads to where F , extracted from the mass conservation ( 42) is with again M = M/M s .Defining the parameter which is central in the planet's hydrostatic theory, allows us to write Ω 2 as Using ( 29) and all simplifications done, we obtain The 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 Eq. 3), our model predicts a Love number where the terms in m can generally be omitted (m = 1/289 for the Earth).The parameter ρ s gR/K s = D/Γ that we discussed in the previous section appears again.We plot in Figure 5a, the Love number h * 0 as a function of planetary mass for a planet having three possible surface densities ρ s = 2000, ρ s = 6000 kg m −3 (same cases as in Figure 3 and Figure 4) and ρ s = 3747 kg m −3 (the value found for the Earth, see Table 1) and m ≪ 1.When the mass is large, h * 0 reaches an asymptotic value, but in this case m = Ω 2 R/g = 3Ω 2 /(4πGρ) goes to zero and the radius becomes independent of the rotation rate (see ( 3)).This is more conspicuous in Figure 5b, where we plot h * 0 ρ s /ρ that includes all the terms dependent of 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 ρ 0 /ρ is very small, 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 3 , 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 Ω e , then to Ω e − δΩ e and to Ω e + δΩ e .Noting the corresponding radii R − and R + , the Love number is approximated by (see Eq. 3) The results of these numerical experiments are also shown in Figure 5a.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.

DISCUSSION AND CONCLUSION
The Lane-Emdden 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 ou 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 depend on the central density ρ c .
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 prop- for p = 2 (Eq.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/ρ, we plot h * 0 ρs/ρ which accounts all the terms depending of 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.
erties 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 a 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); there are not very different from previous estimates.The radius decrease due to the cooling of the Earth (say 250 K in 3 Gyr with α 0 = 3 × 10 −5 K −1 ) should be around 9 km using α e /α s = 0.6 in Figure 4.The Earth's sideral day was only 13 or 15 hours, 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 Ω varied significantly we integrate (3) and use δR = (1/3)h * 0 mR 1 − (Ω 2 A /Ω 2 e ) where Ω A is the Archean rotation rate).Of course, plate tectonics has erased all evidence of these contractions which only reach ≈ 3 µm yr −1 ; only on planets whose lithosphere has been frozen for billions years can thermal contractions be observed.
Our main objective was to show that the important parameters controlling the changes of radius are the dissipation number D = α s gR/C and the Grüneisen parameter Γ = α s K s /(ρ 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 for 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 ≈ 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 those interiors are largely unknown, our approach can provide first order estimates of density profiles and potential changes of the planetary radii through time.

Figure 2 .
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).

Figure 3 .
Figure 3.We show the predicted radii for planets of mass M and surface densities ρs, calculated with p = 3 (panel a) and p = 4 (panel b).The surface densities are 2000 kg m −3 (red dashed line), or 6000 kg m −3 (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 Tables1 and 2).We also show two limiting cases: the case of an incompressible planet of uniform density 4000 kg m −3 (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 bySeager et al. (2007) are shown by blue dashed lines.
Figure4.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 compositionnally 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.

Figure 5
Figure 5. a) Degree 0 rotational Love number as a function of the planetary mass.The curves depict the analyticalfor p = 2 (Eq.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 Table1with 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/ρ, we plot h * 0 ρs/ρ which accounts all the terms depending of 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.

Table 1 .
Mass, radius and average densities of various planets.

Table 2 .
Quantities computed for various planets and exponents p calculated from our model of Lane-Emden planets according to the EoS (4).The masses and radii of the planets are given in Table1.A