1 Introduction

The two ice giants in our solar system are similar in mass and radius, but differ markedly in their thermal evolution: the specific internal luminosity of Uranus is 10 times smaller than that of Neptune (Hubbard et al. 1995). The origin of this difference is unknown, but is of fundamental importance for understanding the thermal evolution of these bodies, as well as for the large number of exoplanets that have similar mass and radius, and which may also be ice-rich (Zeng et al. 2019).

Many previous studies have argued that the internal luminosity of Neptune can be explained by a model of its interior originally developed in the context of the gas giants (Fortney et al. 2011; Nettelmann et al. 2013; Linder et al. 2019). In this view, the interior is fluid and adiabatic throughout the radius of the planet and over the entire course of its evolution. However, a more recent study of Neptune has questioned this view, based on improved equations of state, finding cooling times that are somewhat too short (by 0.9 Gyr) (Scheibe et al. 2019). In the case of Uranus it has been known for some time that fluid, adiabatic models are inadequate as these produce cooling times that are much too long (Nettelmann et al. 2013).

These results have focused attention on alternative models of the interiors of the ice giants. It has long been realized that the low luminosity of Uranus could be explained if some mechanism were found to trap heat deep in its interior (Hubbard et al. 1995) A possible mechanism is the presence of compositional gradients that compete with thermal buoyancy and reduce the efficiency of convective heat transport (Podolak et al. 2019; Vazan and Helled 2020). This possibility has received renewed attention from the discovery of compositional gradients in Saturn revealed by ring seismology (Fuller 2014; Mankovich and Fuller 2021). A number of models now exist that invoke compositional gradients in the ice giants and which can reproduce the cooling times of Uranus and Neptune (Vazan and Helled 2020; Scheibe et al. 2021).

Such models however do not account for tidal dissipation in the interiors of the ice giants because they assume that the interior is inviscid. The evidence for tidal dissipation (\(1/Q\), where \(Q\) is the quality factor) comes from the orbital evolution of moons. Tidal dissipation in the planet causes the moons’ orbits to grow with time (for those moons that are pro-grade and lie outside co-rotation). The rate of recession depends on the mass and orbital radius of the moon, so that different moons recede at different rates. Moons may therefore encounter resonances with each other. On the basis of current orbital characteristics it is possible to infer that certain resonances were encountered in the past, and others avoided. Tittemore and Wisdom (1990) argued that the tidal quality of Uranus is bounded \(Q_{U}\): [11000,39000]; while Zhang and Hamilton (2008) argued that of Neptune is bounded \(Q_{N}\): [9000,36000]. Various mechanisms have been proposed for producing dissipation in a fluid planet (Goldreich and Nicholson 1977; Guenel et al. 2014; Terquem 2023). While none of these mechanisms have yet been quantitatively evaluated for the ice giants, at least some do not appear to be capable of producing dissipation of the magnitude required; for example turbulent viscosity yields \(Q=10^{13}\), which is orders of magnitude too large.

Recent experiments show that the fluid assumption for the ice giants may not be justified. The freezing temperature of H2O at high pressure has been found experimentally to be much higher than previously thought (Millot et al. 2018, 2019), in good agreement with earlier theoretical predictions (Cavazzoni et al. 1999; Mattsson and Desjarlais 2006; Redmer et al. 2011). Most of the interior of Uranus and Neptune fall well below the freezing temperature according to many thermal models of these planets, and within the stability field of a solid, superionic phase. It was argued (Redmer et al. 2011) that the viscosity of the superionic phase is very similar to that of the fluid phase, in which case thermal evolution and tidal dissipation would be no different from the inviscid case. However, more recent ab initio molecular dynamics simulations, combined with proper application of the Green-Kubo relations, show that the superionic phase has a high viscosity (Stixrude et al. 2021).

The experimental results along with a need to account for observational constraints on both thermal and tidal evolution, point towards a new class of models of the ice giant interiors, ones that feature frozen cores. The high viscosity of the frozen core limits the efficiency of convective heat transport, trapping heat in the interior. At the same time, the core provides a source of tidal dissipation.

We have shown that a model of Uranus with a growing, frozen core can account for observations of this planet, including its low internal luminosity, and the evolution of the orbits of its moons (Stixrude et al. 2021). Here we investigate whether the growing frozen core model can account for the properties of Neptune. We review the growing frozen core model, and applications of it to the computation of thermal evolution, tidal dissipation, and the evolution of orbital radii of the moons.

2 Methods

2.1 The Growing Frozen Core Model

At any stage in the evolution of the planet, the frozen core occupies that portion of the interior that lies below the freezing curve (Fig. 1). The outer boundary of the frozen core is set by the condition

$$ T(c)=T_{f}(c)=T_{c} $$
(1)

where \(c\) is the value of the radius \(r\) at the top of the core, \(T\) is the temperature in the planet, and \(T_{f}\) is the freezing temperature. As the planet cools, the frozen core grows

$$ \frac{\partial c}{\partial t} = - \frac{1}{\rho _{c} g_{c} (\Gamma - \Gamma _{a})} \frac{\partial T_{c}}{\partial t} $$
(2)

where we have assumed hydrostatic equilibrium, \(t\) is time, \(\rho _{c}\) and \(g_{c}\) are the density and gravitational acceleration at the top of the core, \(\Gamma \) is the Claussius-Clapeyron slope of the freezing curve, and

$$ \Gamma _{a}=\frac{T}{P}\nabla $$
(3)

is the slope of the planetary temperature distribution. For values of \(\nabla \) typical of ice-giant materials (\(\nabla \)=0.24-0.29) \(\Gamma >\Gamma _{a}\) and \(\partial c / \partial t>0\) on cooling. The planetary temperature profile at any time is

$$ T(P,t) = \textstyle\begin{cases} T_{1} (\frac{P}{P_{1}})^{\nabla} & P < P_{c} \\ T_{i} (\frac{P}{P_{c}})^{\nabla} & P \geq P_{c} \end{cases} $$
(4)

with \(T_{i}=T_{c}+\Delta T\), \(P_{c}\) is the pressure at the top of the core and the subscript 1 refers to the near-surface reference condition with \(P_{1}=1\) bar.

Fig. 1
figure 1

The growing frozen core model showing the frozen core (brick pattern) with radius \(c\), at the top of which develops a thermal boundary layer of temperature contrast \(\Delta T\) and thickness \(\delta \), and overlain by a fluid envelope (blue). The black arrow indicates the growth of the frozen core with time. The fluid layer is assumed to be composed predominantly of H2O over the range of radii occupied by the frozen core, and also includes an outermost gas-rich (H/He) portion

Full size image

The temperature contrast across the thermal boundary layer at the top of the core \(\Delta T\) develops because convection in the frozen core must overcome viscous stresses (Chandrasekhar 1961). If we assume that the Nusselt number of the core is proportional to its Rayleigh number (\({ \mathrm{Ra}}\)) to the one-third power, then \(\Delta T\) is related to the core luminosity

$$ L_{\mathrm{core}}=4 \pi c k \Delta T \left ( \frac{{\mathrm{Ra}}}{{\mathrm{Ra_{cr}}}} \right )^{1/3} $$
(5)

and the thickness of the thermal boundary layer

$$ \delta = c \left (\frac{\textrm{Ra}_{\textrm{cr}}}{\textrm{Ra}} \right )^{\frac{1}{3}} $$
(6)

where \(k\) is the thermal conductivity, and \({\mathrm{Ra_{cr}}}\) is the critical Rayleigh number. The Rayleigh number is a dimensionless quantity that expresses the vigor of convection via the balance between thermal buoyancy and viscous stresses

$$ \textrm{Ra} = \frac{\rho _{c} \alpha g_{c} \Delta T h^{4}}{\eta \kappa c} $$
(7)

where \(\alpha \) is the thermal expansivity, \(\kappa =k \rho C_{P}\) is the thermal diffusivity, \(C_{P}\) is the isobaric heat capacity, \(\eta \) is the viscosity, and \(h=\partial \ln \eta / \partial r\) is the viscous scale height. It is worth emphasizing that in our approach, the properties of the thermal boundary layer, including the temperature contrast and the thickness (Eq. (6)) arise naturally from the physics of high viscosity convection, and are computed self-consistently. Our approach therefore contrasts with that of Scheibe et al. (2021) who assume thermal boundary layers of arbitrary and constant thickness.

The finite viscosity of the frozen core is important because it limits the core luminosity (Eqs. (5), (7)), and also provides a source of tidal dissipation. For the viscosity, we adopt the homologous temperature scaling (Moore 2006)

$$ \eta = \eta _{0} \exp \left [A \left (\frac{T_{f}}{T} - 1 \right ) \right ] $$
(8)

and the freezing temperature increases with pressure

$$ T_{f} = T_{0} \left (1 + \frac{P - P_{0}}{a} \right )^{b} $$
(9)

where \(\eta _{0}\), \(A\), \(a\) and \(b\) are material-specific constants. The viscous scale height is then

$$ h = \frac{T_{c}}{\rho _{c} g_{c} A (\Gamma - \Gamma _{a})} $$
(10)

by combining Eqs. (8), (9).

2.2 Thermal Evolution

The presence of a frozen core means that we must generalize the equations that govern the thermal evolution of purely fluid bodies (Stixrude et al. 2021)

$$ 4\pi R^{2} \sigma (T_{\textrm{eff}} - T_{\textrm{eq}}) = L_{\mathrm{fluid}} + L_{\mathrm{core}} $$
(11)

with

$$ L_{\mathrm{fluid}}=-\int _{M_{c}}^{M} C_{P} \frac{\partial T}{\partial t} dm $$
(12)

and

$$ L_{\textrm{core}} = -\int _{0}^{M_{c}} C_{\textrm{P}} \frac{\partial T}{\partial t}dm - \int _{S} C_{\textrm{P}} \Delta T \frac{\partial c}{\partial t} \rho _{c} dS $$
(13)

where \(T_{\textrm{eff}}\) is the effective temperature of the planet, \(T_{\textrm{eq}}\) is the effective temperature it would have in the absence of internal luminosity, \(R\) and \(M\) are the radius and mass of the planet, \(M_{c}\) is the mass of the core, and \(S\) is the surface of the core. Early in the planet’s evolution, the interior temperature \(T>T_{f}\) everywhere, \(L_{\mathrm{core}}=0\), and the thermal evolution equations reduce to that of a purely fluid body (Fortney and Nettelmann 2010). A relationship between \(T_{\mathrm{eff}}\) and \(T_{1}\) closes the equations (Guillot et al. 1995)

$$ T_{1} = B T_{\textrm{eff}}^{1.244} $$
(14)

with \(B = 0.44838\) reproducing the present-day observed value of \(T_{1} = 72\) K for \(T_{\textrm{eff}} = 59.3\) K. We assume that the adiabatic gradient (away from thermal boundary layers) and the heat capacity are homogeneous. We have shown that a homogeneous value of \(\nabla \) captures very well the temperature structure of detailed multilayer models in which the adiabatic gradient is given at every point by an assumed composition and equations of state of rock, ice, and gas components (Stixrude et al. 2021). We do not consider the effects of the latent heat of freezing because we showed in our study of Uranus that its effects on the thermal evolution of the planet are negligible (Stixrude et al. 2021).

2.3 Tidal Evolution

The quality factor of the planet

$$ Q = \frac{\left | k_{2} \right | }{\mathrm{Im} ( k_{2} ) } $$
(15)

where \(k_{2}\) is the (complex) tidal Love number. We compute \(k_{2}\) and \(Q\) at each stage of the planet’s evolution from the six coupled differential equations that govern the visco-elastic response of the planet to imposed tidal stress, including the radial and tangential displacement, the radial and shear stresses, and the gravitational potential and its radial derivative. This set of equations is widely used in studies of viscoelastic planetary bodies (Alterman et al. 1959; Sabadini and Vermeersen 2004; Matsuyama et al. 2018; Henning and Hurford 2014), and our implementation is discussed in detail in Stixrude et al. (2021). In order to compute the tidal response from these equations, we must specify the density and the complex shear modulus as a function of radius within the planet. The complex shear modulus depends on the viscosity already specified and the infinite frequency shear modulus \(\mu _{0}\). For fluid layers, we take \(\mu _{0}\) to be asymptotically small; for the frozen core, we adopt the value from ab initio molecular dynamics simulations (Hernandez and Caracas 2016).

$$ \mu _{0}(P,T)=101+P/2.41-(T-1650)/23.4 $$
(16)

with pressure and shear modulus in GPa, and temperature in K. The complex shear modulus according to the standard linear solid is then (Nowick and Berry 1972)

$$ \mu = \mu _{0} - \frac{\delta \mu}{1 + i \omega \tau _{\epsilon}} $$
(17)

where

$$ \delta \mu = \mu _{0} (1 - \frac{\mu _{1}}{\mu _{0} + \mu _{1}}) $$
(18)

and the characteristic material timescales are

$$ \tau _{\epsilon} = \frac{\eta}{\mu _{0} + \mu _{1}} \qquad \text{and} \qquad \tau _{\sigma} = \frac{\eta}{\mu _{1}} $$
(19)

where \(\eta \) is the viscosity Eq. (8). The shear modulus is evaluated at the frequency \(\omega =2 \left | \chi -n \right |\), where \(\chi \) is the rotation rate of the planet and \(n\) is the mean motion of the moon.

2.4 Orbital Dynamics

The planet’s moons migrate radially as the planet evolves according to (Tittemore and Wisdom 1989)

$$ \frac{\dot{a}}{a} = \pm 3 \frac{k_{2}}{Q} n \frac{m}{M} \left ( \frac{R}{a}\right )^{5} $$
(20)

where \(a\) is the orbital radius of the moon, \(\dot{a}\) is the time derivative of \(a\), \(n\) is the mean motion, \(k_{2}\) is the tidal love number, \(m\) is the mass of the moon, \(M\) is the mass of the planet, \(R\) is the radius of the planet, and the sign is that of \(\omega \) for prograde moons and \(-\omega \) for retrograde moons. We integrate Eq. (20) backward in time from present-day observed values of \(a\) (Jacobson 2009). For some of the moons of Neptune, the mass is unknown. For these moons, we combine the observed volumes with an estimate of their density \(\rho _{\text{moon}}=0.8\) g/cm3 (Buratti and Thomas 2014).

Our application differs from many previous studies of migration in that the values of \(Q\) and \(k_{2}\) change with time as the planet thermally evolves. In the earliest stages of planetary evolution, the planet is entirely fluid so that \(1/Q \rightarrow 0\), and the migration is negligible. Only as the frozen core nucleates and grows do the attenuation and the orbital migration become significant.

3 Application to Neptune

We use a model of the interior of Neptune that agrees with all relevant observational data (mass, radius, \(J_{2}\), \(J_{4}\)): the so-called empirical model of Helled et al. (2011). For comparison, we also consider another model that includes a dense core and matches observations equally well (Kaspi et al. 2013). We do not account explicitly for the possibility of compositional layering in Neptune, but assume that over the range of radii that the frozen core exists, the composition is dominated by H2O. As in our previous study of Uranus, we have considered the effect of a possible inner-most rocky core, with a much higher freezing temperature and find that this rocky core is too small to affect our results. Other parameters of our model are shown in Table 1. Material properties are the same as in our study of Uranus (Stixrude et al. 2021) and references therein, and are those of H2O, except for \(C_{P}\) and \(\nabla \), for which we consider a range of possible values. The range of values of \(C_{P}\) and \(\nabla \) that we consider is identical to that in our study of Uranus (Stixrude et al. 2021).

Table 1 Values of parameters used in the thermal and tidal evolution
Full size table

We find that the growing frozen core model is able to account for the thermal evolution of Neptune (Fig. 2). We gain additional insight into the role of the frozen core in thermal evolution by comparing with the purely fluid case, which we compute using the same equations, but with the value of \(\eta _{0}\) assumed to be arbitrarily small. We find that early in the planet’s history, when the core is either absent or too small to significantly affect heat transport, thermal evolution is identical to that of the purely fluid case. Frozen core and fluid cases begin to diverge significantly in their thermal evolution near 3 Gyr, as the core retains an increasing proportion of internal heat.

Fig. 2
figure 2

Thermal evolution of Neptune with the growing frozen core model (blue) and fluid model (red) with parameters \(\nabla = 0.281\) and \(C_{\textrm{P}} = 5000\). The top inset shows the core radius and the bottom inset shows the thickness (dashed curve) and temperature contrast (solid curve) across the thermal boundary layer

Full size image

Initially, the planet is entirely fluid with \(T>T_{f}\) throughout the radius: in our thermal evolution calculation, the initial temperature is arbitrarily hot (Figs. 2, 3). The frozen core nucleates at 2.2 Gyr when \(T(0)=T_{f}(0)\). The rate of growth is initially rapid because the rate of cooling is more rapid at earlier times and because the Clapeyron slope \(\Gamma \) is smaller at high temperature (Eq. (9)). The frozen core grows to occupy 2/3 of the planet’s radius. The core radius and \(\Delta T\) change little upon subsequent evolution after 4.5 Gyr. As the core grows, the temperature contrast and thickness of the thermal boundary layer also grow, to values of 1100 K, and 16 km, respectively at 4.5 Gyr.

Fig. 3
figure 3

The temperature profiles of Neptune in the growing frozen core model at various snapshots in time as indicated

Full size image

We find that a wide range of values of the heat capacity and adiabatic gradient can account for the cooling time of Neptune (Fig. 4). The composition of the interior of Neptune cannot be uniquely constrained by current observations. We have therefore considered a range of possible values \(C_{P}\)=3000-6000 J/g/K, which encompasses the Dulong-Petit values of pure H2O, a solar mix of water, methane, and ammonia (Asplund et al. 2009), and the possible effects of admixtures of hydrogen, helium, and/or heavier elements, including “synthetic Uranus” (Hubbard et al. 1991). Our assumed range of adiabatic gradient \(\nabla \)=0.24-0.29 encompasses detailed multi-layer models of the ice giants (Nettelmann et al. 2013; Bethkenhagen et al. 2017). Higher values of the heat capacity and the adiabatic gradient increase the thermal inertia of the planet, leading to slower calculated cooling times. The frozen core model cools more quickly than the fluid case because the presence of the core traps heat at depth, allowing for the fluid envelope to cool more quickly.

Fig. 4
figure 4

Cooling time of Neptune as a function of \(\nabla \) for various values of \(C_{P}\) (color bar in units of J/g/K) for frozen core (solid lines) and for fluid (dashed lines) cases

Full size image

The tidal quality of Neptune varies significantly over the course of its evolution (Fig. 5). The reason is that the planet starts in a completely fluid state with very large \(Q\). As the frozen core nucleates and grows, \(Q\) drops by many orders of magnitude, reaching a minimum value of 5000 at 3.5 Gyr. The value of \(Q\) subsequently recovers slightly as the frozen core grows to its maximum size and an increasingly large proportion becomes so cold (\(T/T_{f} \ll 1\)) that it no longer contributes significantly to dissipation.

Fig. 5
figure 5

The tidal evolution of Neptune with the tidal love number \(k_{2}\) (blue, left-hand axis) and the tidal quality \(Q\) (red, right-hand axis)

Full size image

The tidal Love number also varies with time because \(k_{2}\) is sensitive to the variation of shear modulus with depth (Fig. 5). In the initial, fluid, stage of evolution \(k_{2}=0.41\). As the frozen core grows, the tidal response of the planet is reduced by the rigidity of the frozen core, and the value of the tidal love number diminishes significantly to \(k_{2}=0.29\).

The orbits of the moons of Neptune respond to the evolving dissipation in a way that explains their orbital characteristics (Fig. 6). The orbital evolution of the moons is very slow at first as the planet is entirely fluid and \(Q\) is very high. Only once the frozen core begins to form do the orbital radii begin to vary significantly. Moons outside the co-rotation limit recede with time, while those within approach the planet. Triton approaches Neptune with time because its orbit is retrograde. The moons encounter many resonances, and our computed evolution is consistent in all respects with the scenario envisioned by Zhang and Hamilton (2008). For example, Proteus and Galatea encounter a 2:1 resonance that Zhang and Hamilton (2008) argue explains the inclination of Larissa. Naiad and Despina encounter a 9:11 resonance that Banfield and Murray (1992) argue explains the large inclination of Naiad. Proteus and Larissa do not encounter the co-rotation limit and avoid a 3:2 resonance with each other, which Zhang and Hamilton (2008) argue would have caused the inclination of Proteus to be larger than observed. We predict present day rates of change of orbital radii of the moons as summarized in Table 2. The rate of evolution of the moons’ orbits depends on their density, which is uncertain. For our fiducial case, we have assumed a uniform density of all the small moons \(\rho _{\text{moon}} = 0.8\) g cm−3, consistent with the study of Buratti and Thomas (2014). We have also considered a range of other possible values, and find that the set of resonances that appear is unchanged for values \(0.4 < \rho _{\text{moon}} < 1.0\) g cm−3, while for \(\rho _{\text{moon}} > 1.0\) g cm−3, the 3:2 resonance between Proteus and Larissa appears.

Fig. 6
figure 6

Orbital evolution of Neptune’s moons. The dashed line indicates co-rotation. The vertical black lines indicate resonances discussed in the text: Proteus-Galatea 1:2 and Despina-Naiad 9:11. The orbital radii are shown in units of Neptune’s radius \(R_{N}\)

Full size image
Table 2 Present-day Orbital Migration
Full size table

4 Discussion

In our model, the growing frozen core arises naturally from the thermal history of the planet. The frozen core influences thermal evolution because of its high viscosity and therefore inefficient heat transport, while at the same time providing a time-dependent source of tidal dissipation. The presence of frozen cores within one or more of the outer planets has been proposed before, and it has been recognized that these frozen cores may be an important source of tidal dissipation (Dermott 1979; Lainey et al. 2017; Mazevet and Guyot 2019). However, these previous studies have not placed the frozen core within the context of a self-consistent prediction of thermal and tidal evolution as we have done.

An interesting consequence of the growing frozen core is that we find that the rate of migration of the orbits of the moons changes markedly with time (Fig. 6). The present day rate of migration is much greater than the rate averaged over the entire history of the planet. This presents a very different perspective than many analyses of migration. In many studies, \(Q\) of the planet is assumed to be constant with time, and the rate of migration therefore varies only slightly with time according to Eq. (20) (Zhang and Hamilton 2008; Banfield and Murray 1992). Instead, in the growing frozen core model, the value of \(Q\) changes by many orders of magnitude over the planet’s history. In the case of Neptune, over most of the evolution, the value of \(Q\) is so high, that the rate of migration is negligible. It is only in the last 1.5 Gyr that the rate of migration increased due to growth of the frozen core.

Time-dependent \(Q\), which is a natural result of a growing frozen core, could provide a new perspective on problems related to the migration of the moons of other planets. For example, in the case of Saturn, the present day migration rate of many of its moons is so rapid that, if one assumes that \(Q\) is independent of time, the moons must be billions of years younger than Saturn (Lainey et al. 2012). Another possibility though is that \(Q\) of Saturn is much lower today than it was in the past. Such would be the case if dissipation originated in a frozen core that has grown with time as in our model, although the core in Saturn would make up a much smaller fraction of the planet’s radius (Lainey et al. 2017). Another proposed mechanism (resonance locking) envisions a different source of the dissipation in Saturn, but also emphasizes the importance of considering the planet’s thermal evolution in understanding the migration of the moons (Fuller et al. 2016).

The presence of a frozen core in the ice giants is important for understanding the generation of their magnetic fields. The magnetic field geometry of Uranus and Neptune is known to be significantly non-dipolar (Connerney 1993). It has been suggested that the complex field geometry may be more readily produced if the dynamo generating region is either thin and/or close to the surface (Soderlund and Stanley 2020). The frozen core is far too viscous to generate a magnetic field. Moreover, the dynamo generating region must be electrically conductive, further limiting its thickness. In our model then, the dynamo is generated in the relatively thin and shallow fluid ice layer that lies between the frozen core and the H-rich envelope. Taking the boundaries of this layer to be the top of the frozen core and the base of the H-rich layer in model N1 Nettelmann et al. (2013), the dynamo is generated in the layer: \(0.68-0.86 R_{N}\), with a total thickness of 4000 km. While the frozen core is too viscous to generate a magnetic field of its own, its electrical conductivity is comparable to that of fluid ice (Grasselli et al. 2020). The frozen core may therefore affect the geometry of the magnetic field generated in the layer above.

We have shown here and in our previous study (Stixrude et al. 2021) that the growing frozen core model can explain the thermal and tidal evolution of Uranus and Neptune. For both planets, we find a range of values of \(C_{P}\) and \(\nabla \) that match their respective cooling times, and constraints on their tidal evolution, and which fall within the range of previous models of these planets (Nettelmann et al. 2013; Scheibe et al. 2019; Bethkenhagen et al. 2017). The ranges of values that match the properties of the two planets do not overlap. For example, for the range of heat capacity values that we consider (3000-6000 J/g/K), the range of values of \(\nabla \) that match the cooling time are 0.257-0.263 for Uranus and 0.276-0.294 for Neptune. We interpret this difference as arising from a contrast in composition between the two planets, as seen in many previous models of Uranus and Neptune (Nettelmann et al. 2013; Scheibe et al. 2019). While presently available observations cannot pin down uniquely the nature of these compositional differences, they may include, for example, differences in the relative proportions of gas to ice to rock, or the composition of the ice (relative proportions of H2O, NH3, CH4). Our results therefore point to directions of future study, such as including the effects of a growing frozen core in thermal evolution models that consider a detailed representation of the equations of state of gas, ice, and rock (Scheibe et al. 2019), computing the tidal evolution from such models, and future missions that may be able to narrow the range of allowable compositions in the two planets (Movshovitz and Fortney 2022).

The growing frozen core model, inspired in part by the experimental observation that the freezing temperature of H2O is much higher than previously thought, emphasizes the importance of experimental measurements and quantum mechanical simulations for understanding planetary evolution, while also highlighting the need for further research into material behavior at extreme conditions. For example, while our assumption that the ice region of the planet is dominated by H2O may be reasonable, a better understanding of the effects of ammonia and methane on the freezing curve and the superionic transition will be important for future studies (Bethkenhagen et al. 2015; Robinson and Hermann 2020; Cheng et al. 2023). It will also be important to build on predictions of physical properties, including heat capacity, adibatic gradient, thermal and electrical conductivity, and shear modulus of crystalline phases to include a broader range of pressure and temperature and the effects of admixtures of ammonia and methane (Bethkenhagen et al. 2017; Grasselli et al. 2020; Millot et al. 2019; Hernandez and Caracas 2016).

Our results are subject to observational tests by a future space mission. The tidal Love number \(k_{2}\) is diagnostic of the presence of a frozen core because, in contrast to the moments of the gravitational field \(J_{2n}\), it is sensitive to the shear modulus. The value of \(k_{2}\) has now been measured for Jupiter and Saturn, but not yet for Uranus or Neptune. We find that the value of \(k_{2}\) with and without a frozen core, respectively differ by 0.12 or 30%; we previously found a difference of similar magnitude for Uranus (Stixrude et al. 2021). This difference is much larger than the uncertainty on the value of \(k_{2}\) for Jupiter (3%) (Durante et al. 2020) or Saturn (6%) (Lainey et al. 2017). We also predict the present-day rate of migration of the orbits of the moons, which could also be measured by a future space mission (Table 2).

5 Conclusion

We have shown that the growing frozen core model can account for the thermal and tidal evolution of the ice giants of our solar system in this study of Neptune and our previous study of Uranus (Stixrude et al. 2021). It is often assumed that the ice giants are entirely fluid, but no currently available observation requires this to be the case. In fact, inviscid models are difficult to reconcile with experimental data and with observations of the orbits of the ice giant moons. Our results can be tested by future space missions, including the proposed Uranus Orbiter and Probe and by a possible future mission to Neptune.