# The Effect of Accretion Rate and Composition on the Structure of Ice-rich
Super-Earths

Michael Lozovsky Department of Geosciences,
Tel-Aviv University, Tel-Aviv,
Israel  Dina Prialnik  Department of Geosciences,
Tel-Aviv University, Tel-Aviv,
Israel  Morris Podolak Department of Geosciences,
Tel-Aviv University, Tel-Aviv,
Israel

###### Abstract

It is reasonable to assume that the structure of a planet and the interior
distribution of its components are determined by its formation history. We
thus follow the growth of a planet from a small embryo through its subsequent
evolution. We estimate the accretion rate range based on a protoplanetary disk
model at a large enough distance from the central star, for water ice to be a
major component. We assume the accreted material to be a mixture of silicate
rock and ice, with no H-He envelope, as the accretion timescale is much longer
than the time required for the nebular gas to dissipate. We adopt a thermal
evolution model that includes accretional heating, radioactive energy release,
and separation of ice and rock. Taking the Safronov parameter and the ice-to-
rock ratio as free parameters, we compute growth and evolutionary sequences
for different parameter combinations, for 4.6 Gyr.

We find the final structure to depend significantly on both parameters. Low
initial ice to rock ratios and high accretion rates, each resulting in
increased heating rate, lead to the formation of extended rocky cores, while
the opposite conditions leave the composition almost unchanged and result in
relatively low internal temperatures. When rocky cores form, the ice-rich
outer mantles still contain rock mixed with the ice. We find that a
considerable fraction of the ice evaporates upon accretion, depending on
parameters, and assume it is lost, thus the final surface composition and bulk
density of the planet do not necessarily reflect the protoplanetary disk
composition.

planets and satellites: composition — planets and satellites: formation —
planets and satellites: interiors — planet–disk interactions — planets and
satellites: physical evolution

##  1 Introduction

We now have observational data on more than 5,000 exoplanets. These planets
show very large diversity in their measured properties, such as mass and
radius. While the majority of these planets have a Jupiter mass or more,
super-Earth-sized planets are also very abundant (e.g. Fulton et al., 2017;
Fressin et al., 2013; Batalha et al., 2013; Bean et al., 2021). Gas giants are
expected to have large gaseous envelopes (as their name implies), and it is
not surprising that in a given mass range the radii of these bodies show a
large variation (e.g. Laughlin, 2018; Hou & Wei, 2022). Smaller bodies,
however, with masses of only a few $M_{\oplus}$, are less likely to accumulate
large gaseous envelopes, and we would expect to see a correspondingly smaller
variation in their radii. Yet, observed radii in this mass range can vary
widely (see, e.g. Zeng & Sasselov, 2014, Figure 3).

Many factors contribute to this variation. Obvious candidates are the presence
or absence of a significant gaseous envelope, differences in internal
composition, and surface temperature as fixed by the distance from the central
star and temperature of that star. However, there are also more subtle factors
that are less commonly considered. Zeng & Sasselov (2014) have pointed out
that the slow cooling of a water-rich planet will cause noticeable changes in
radius as a function of time. Thus it is of interest to understand how the
initial temperature and composition profiles of a planet might be dependent on
the circumstances of its formation.

The connection between planetary origin and internal structure, is nontrivial,
and different formation mechanisms and birth environments can lead to a large
range of compositions and internal structures (e.g. Helled et al., 2014;
Lozovsky et al., 2017; Valletta & Helled, 2020, and references within). Even
the host stellar type might influence the planetary structure and composition
(Lozovsky et al., 2021; Silva & Valio, 2016; Adibekyan et al., 2021). While
there are numerous observations of protoplanetary disks (e.g. Andrews, 2020),
and various aspects of their evolution have been modeled theoretically, the
properties of such protoplanetary disks are still not well constrained.
Indeed, we still don’t have a complete picture of the evolution of even our
own proto-solar nebula (e.g. Villenave et al., 2021).

In this paper we address one very specific aspect of the connection between
formation processes and the resultant planetary structure: For a planet formed
by accretion of water and rock, how do the temperature and composition
profiles depend on the accretion rate and planetesimal composition? In view of
the uncertainties involved, we consider a simple model of the protoplanetary
disk and parameterize the associated accretion process. We limit ourselves to
the case where the planet is formed far from the central star (40 au) and
assume a simple model of planetary accretion. In this way the amount of solid
material for planet formation is low enough so that the planet does not reach
the stage of rapid gas accretion before the disk dissipates (Pollack et al.,
1996; Alibert et al., 2005; Frelikh & Murray-Clay, 2017). In addition, the
protoplanetary embryo remains small enough while the gas is in place, so that
migration is not an issue (Armitage, 2020; Ogihara et al., 2022). In the case
of pebble accretion things will happen more quickly and we address this in
section 5.

In Section 2 we describe the nebular model and the parameterization of the
accretion process, in Section 3 we describe the thermal evolution model used
to evolve the planet, in Section 4 we present our results, and in Section 5 we
present our conclusions.

##  2 Accretion model

###  2.1 Protoplanetary Disk

We assume a sun-like star (1 M⊙) surrounded by a protoplanetary disk with a
surface density of gas given by (Hayashi, 1981)

$\Sigma_{g}(a)=\Sigma_{0}a^{-3/2}$ (1)

where $a$ is the radial distance in au, and $\Sigma_{0}$ is the value of
$\Sigma_{g}$ at $a=$1 au. If we assume that the disk extends to 50 au and has
a mass of $M_{d}=0.02M_{\odot}$, then we have

$M_{d}=2\pi\int_{a_{1}}^{a_{2}}r\Sigma_{g}(r)dr$ (2)

where $a_{1}$ and $a_{2}$ are the inner and outer radii of the disk,
respectively. Setting $M_{d}=4\times 10^{31}$ g and taking $a_{1}=0.01$ au and
$a_{2}=50$ au in Equation (2) gives $\Sigma_{0}=2\times 10^{3}$ $\rm
g~{}cm^{-2}$. The value of $\Sigma_{0}$ is insensitive to the choice of
$a_{1}$ so long as $a_{1}\ll a_{2}$.

We consider the region at $\sim$ 40 au, far from the influence of other
massive bodies which might be forming, and at a low enough temperature so that
almost all the material other than hydrogen and helium is in the form of
solids. We further assume that the mass fraction of these solids is 0.02, in
accordance with solar composition (Lodders, 2010). This gives a surface
density of solids, $\Sigma_{s}\approx 0.16$ $\rm g~{}cm^{-2}$. In order to
parameterize the growth rate in a simple, yet physically reasonable manner, we
assume that the proto-planetary embryo grows by accreting the material in an
annulus centered at 40 au and having a width of $\Delta a$ on either side,
where $\Delta a$ is related to the Hill’s sphere radius of the final body.
This isolation mass is given by Lissauer (1987) as

$M_{iso}=1.95\times 10^{24}\left(Ba^{2}\Sigma_{s}\right)^{3/2}$ (3)

where, again, $a$ is in au, $B\sim 3-4$ is a dimensionless constant, and all
the rest of the parameters are in cgs units. Taking $B=3$ gives an isolation
mass at 40 au of $4.2\times 10^{28}$ g or about $7\,M_{\oplus}$.

###  2.2 Mass Accretion Rate

If $R$ is the body’s radius, $\rho_{s}$ is the space density of solids in the
region, and $v$ is the random velocity of the planetesimals with respect to
the growing embryo, the rate of growth is given by (Lissauer, 1987)

$\dot{M}=\frac{dM}{dt}=\pi
R^{2}\rho_{s}v\left[1+\left(\frac{v_{e}}{v}\right)^{2}\right]$ (4)

where $v_{e}$ is the escape velocity from the embryo and the term in brackets
is the gravitational enhancement factor (see Safronov, 1972). This can be
rewritten in terms of the surface density as (Lissauer, 1987)

$\dot{M}(t)=\frac{\sqrt{3}\pi}{2}R(t)^{2}\Sigma_{s}(t)\Omega\left[1+2\Theta\right]$
(5)

where $\Theta=v_{e}^{2}/2v^{2}$ is the so-called Safronov parameter, and
$\Omega$ is the Kepler frequency. Lissauer (1987), based on a generalization
of the results of Wetherill & Cox (1985), gives

$\Theta\approx 400a\left(\frac{\bar{\rho}}{4.04}\right)^{1/3}$ (6)

where $\bar{\rho}$ is the mean density of the planetary embryo in cgs, and $a$
is again in au. If we ignore the very weak dependence on $\bar{\rho}$, we find
that at 40 au the expected value is $\Theta=16,000$. We assume that the total
mass (planetesimals + growing embryo) in this annulus remains constant with
time. In what follows, we allow $\Theta$ to vary around its fiducial value in
order to explore how different accretion rates affect the temperature and
composition profiles for a given mass. It should be noted that alternative
formation mechanisms, such as pebble accretion, will yield a very different
$\dot{M}$ profile (e.g. Bitsch et al., 2015; Andama et al., 2021; Ormel et
al., 2021), and we address this issue in the discussion section. We define a
normalized Safronov parameter for this study $\hat{\theta}=\frac{\Theta}{a}$,
than a is in au units.

###  2.3 Composition

As the planet grows in a region of $\sim$ 40 au, we assume that the
planetesimals are composed of rock (SiO2) and water ice (H2O). We do not take
into account H-He gas, since the prescription we have chosen for $\dot{M}$
assures that the nebular gas will dissipate while the mass of the
protoplanetary embryo is still very small. For simplicity, heavier elements
such as metals are not considered, as their contribution to the mass-radius
relation is weak (Vazan et al., 2018). Thus, the rock is assumed to be one
single substance, and we ignore possible chemical reactions or differentiation
within the rock, such as separation of iron from silicates. The ice is also
taken as one substance, thus phase transitions between various forms of solid
ice are ignored and only melting/freezing is taken into account.

The ice mass fraction we investigate ranges from 0.3 to 0.7 out of the total
mass of the material available for accretion in the disk. After accretion, the
ice and rock distribution depends on the tendency of those two substances to
mix or separate in different regimes. A separation between ice and rock may
lead to formation of layers, based on their densities (Stevenson, 2013; Vazan
et al., 2020). Since the phase separation between rock and ice is temperature
and pressure dependent (e.g. Dorn & Lichtenberg, 2021; Vazan et al., 2020), we
investigate how different mass fractions of ice influence the formation of
different cores. The composition of the accreted material is the same
throughout the accretion process, but the radioactive content corresponds to
the accretion time, that is, diminishes with time as in the planet itself.

##  3 Thermal Evolution Model

We consider a spherical body that grows from an initial embryo of small mass
compared to the final mass. The composition of both initial embryo and
accreted material is a mixture of rock and ice. After $\hat{\theta}$, the
ratio of ice to rock is the second free parameter of this study. The mass
fraction of ice is denoted $X_{i}$ and that of rock, $X_{r}=1-X_{i}$. If
melting occurs, a fraction of $X_{i}$ will be liquid; it is taken as a
smoothed step function of the temperature around the pressure-dependent
melting temperature $T_{m}(P)$, $X_{\ell}=X_{i}/[1+{\rm
e}^{\beta(1-T/T_{m})}]$ (see Prialnik & Merk, 2008, for further details). When
the ice is heated and eventually melts, its viscosity $\eta(T,T_{m})$
decreases and the rock settles toward the center, allowing ice-rock separation
to occur. This is modeled by a downward flux of rocky material ${\bf J}_{r}$
and an equal upward flux of ice/water ${\bf J}_{i}$; we assume the relative
movement to be symmetrical with respect to the local center of mass. Thus
ice/water is emplaced by rock and thermal energy is exchanged in the process,
since the specific energies of ice/water and rock differ. The relative
velocity of rock and ice is taken to be

$v_{r}={\textstyle{2\over 9}}(\rho_{r}-\rho_{i})gr_{p}^{2}/\eta,$ (7)

(Stokes flow) where $\rho_{r}$ and $\rho_{i}$, are the densities of rock and
ice, respectively, $g$ is the gravitational acceleration and $r_{p}$ is a
typical radius of a rock particle. Since the movement is slow, it is assumed
that thermodynamic and hydrodynamic equilibrium are maintained, thus there is
a unique local temperature and the structure evolves quasi-statically.

When a body grows in mass by accretion, it also gains thermal energy from the
kinetic energy of the accreted material. Since the relative velocity of
planetesimals entering the activity sphere of the accreting body is lower than
the free-fall velocity, the rate of gain of accretion energy is only a
fraction of $GM\dot{M}/R$ (Bodenheimer & Pollack, 1986; Boujibar et al.,
2020). But even then, not all the kinetic energy turns into heat; a large
fraction goes into mechanical energy, such as excavating a crater or breaking
rocks (Ransford, 1982). As working hypothesis, we assume the fraction that
turns into heat to be ${\textstyle{1\over 4}}$ of the kinetic energy obtained
if the accreted material free-falls from infinity (c.f. Bodenheimer & Pollack,
1986; Pollack et al., 1986). We further assume that the accretion energy is
divided between a surface and a body source (Squyres et al., 1988). A fraction
$f$ is taken to be deposited at the surface, while the remaining fraction
$(1-f)$ is absorbed in the interior, the effect diminishing with depth. We
adopt $f=0.4$, but the results are not very sensitive to the exact value.

We solve the coupled energy and mass conservation equations. For energy,

$\frac{\partial(\rho u)}{\partial t}+{\nabla}\cdot{\bf
F}+{\nabla}\cdot(u_{r}{\bf J}_{r}+u_{i}{\bf J}_{i})=\dot{Q}$ (8)

where $\rho u$ denotes the sum of weighted specific energies, ${\bf F}$ is the
heat flux, ${\bf F}=-K{\nabla}T$, $K$ denoting the thermal conductivity, and
$\dot{Q}$ includes all energy sources: insolation, accretion, compression by
gravitational forces, latent heat (which may be negative), and radioactive
decay energy by long-lived radioisotopes embedded in the rock. With the rate
of growth of the planetary mass $\dot{M}(t)$ given by Equation 4, the surface
boundary condition becomes

$F(R)=\sigma T^{4}+Z(T)\mathcal{H}_{v}-(1-A)\frac{L}{4\pi
a^{2}}-{\textstyle{1\over 4}}{}f\frac{GM\dot{M}}{4\pi R^{3}},$ (9)

where $\sigma$ is the Stefan-Boltzmann constant, $A$ is the albedo, $L$ \- the
stellar luminosity, and $a$ \- the distance from the star. If the surface
temperature is sufficiently high, the ice may sublimate; $Z(T)$ is the
sublimation rate determined by the saturated vapor pressure, and
$\mathcal{H}_{v}$ is the corresponding latent heat absorbed. For mass,

$\frac{\partial(\rho X_{r})}{\partial t}+{\nabla}\cdot{\bf J}_{r}=0$ (10)
$\frac{\partial(\rho X_{i})}{\partial t}+{\nabla}\cdot{\bf J}_{i}=0.$ (11)

At the surface, $J_{i}=Z(T)$. Finally, the bulk density profile at any time
$t$ is given by the hydrostatic equation,

${\nabla}P(\rho)-\rho{\bf g}=0,$ (12)

where $P(\rho)$ is the pressure given by the equation of state.

For a given pressure $P$ and ice mass fraction $X_{i}$, the bulk density
weighted by the volume fractions of the components is given by

$\frac{1}{\rho}=\frac{X_{i}}{\rho_{i}(P)}+\frac{1-X_{i}}{\rho_{r}(P)},$ (13)

where $\rho_{i}$ and $\rho_{r}$ are the ice and rock densities, respectively,
for each of which we adopt the second-order Birch-Murnaghan (BM)
approximation. The coefficients are chosen so as to achieve a close agreement
with the tabulated equation of state (EOS) developed and used by e.g., Vazan
et al. (2013, 2015), based on the quotidian EOS (QEOS) of More et al. (1988).
Hence $\rho_{i}$ is obtained by solving

$P_{0,i}\left[\left(\frac{\rho_{i}(P)}{\rho_{0,i}}\right)^{7/3}-\left(\frac{\rho_{i}(P)}{\rho_{0,i}}\right)^{5/3}\right]-P=0,$
(14)

where $\rho_{0,i}=0.917$ g/cm3, $P_{0,i}=2.07\times 10^{11}$ dyn/cm2, and
similarly, $\rho_{r}$ is obtained by solving

$P_{0,r}\left[\left(\frac{\rho_{r}(P)}{\rho_{0,r}}\right)^{7/3}-\left(\frac{\rho_{r}(P)}{\rho_{0,r}}\right)^{5/3}\right]-P=0,$
(15)

where $\rho_{0,r}=2.5$ g/cm3 and $P_{0,r}=1.28\times 10^{12}$ dyn/cm2. In the
range of interest for densities and temperatures, the effect of temperature on
pressure according to the QEOS is very small. We have also estimated the Debye
correction to the BM EOS and found it small. Therefore, we assume the pressure
to be independent of temperature for the entire relevant temperature range.

The components of the volume energy source $\dot{Q}$ (energy per unit time per
unit volume) in Eq. (8) are as follows:

The accretional heating source for the interior is

$\dot{Q}_{\rm acc}={\textstyle{1\over 4}}(1-f)\frac{3GM(t)\dot{M}(t)}{4\pi
R^{4}(t)}h(z),$ (16)

where $h(z)$ is an exponentially decreasing function of the normalized depth
$z$, such that $\int_{0}^{1}h(z)dz=1$.

The radioactive energy source is

$\dot{Q}_{rad}=\rho_{r}\sum_{j}\tau_{j}^{-1}X_{0,j}H_{j}{\rm
e}^{-t/\tau_{j}},$ (17)

where $X_{0,j}$ is the initial abundance (relative to rock) of the
radionuclide, $H_{j}$ \- the energy released per unit mass, and $\tau_{j}$ \-
the characteristic decay time. We consider 40K, 232Th, 238U and 235U.

The heat absorbed/released by melting/refreezing is given by

$\dot{Q}_{\ell}=-\rho\frac{\partial X_{\ell}}{\partial t}H_{\ell}=-\rho
X_{i}\frac{d}{dT}\left[\frac{1}{{\rm
e}^{\beta(1-T/T_{m})}}\right]\frac{\partial T}{\partial t}H_{\ell},$ (18)

where $H_{\ell}$ is the latent heat of melting.

As we are considering a spherically symmetric body, we choose the volume
enclosed by a spherical surface of radius $r$ ($0\ \leq r\ \leq R$), denoted
by $V$ ($0\ \leq V\ \leq 4\pi R^{3}/3$) as the independent space variable.
Thus mass and energy fluxes are replaced by energy or mass crossing a
spherical surface per unit time, and

${\nabla}\ \Longrightarrow\frac{\partial}{\partial V}.$

In the numerical computations, where the equations are discretized in space,
and the spatial grid is adapted to the varying configuration of the model, we
consider a different, dimensionless space variable $x$, defined over a finite
range $[c,s]$, where $c$ and $s$ are the system’s boundaries (center and
surface). In this case $V(x)$ must be supplied as a monotonically increasing
solution of an equation, satisfying $V_{c}=0$ and $V_{s}=V_{s}(t)$ (moving
boundary). The equation is so chosen as to ensure fine zoning where steep
gradients (of temperature or other physical property) arise and has the
general form

$\frac{\partial f(V)}{\partial x}={\rm constant}.$ (19)

Although the range of $x$ is fixed, the total volume $V(s)$ changes with time
due to accretion of material, $\dot{V}(s)=\dot{M}/\rho(s)$, and consequently
$V(x)$ changes at all $x$. Since temporal derivatives are taken at constant
$V$, whereas $V=V(x,t)$, the following transformation is implemented in the
difference scheme:

$\left(\frac{\partial}{\partial t}\right)_{V}=\left(\frac{\partial}{\partial
t}\right)_{x}-\left(\frac{\partial V}{\partial
t}\right)_{x}\left(\frac{\partial}{\partial V}\right)_{t}.$ (20)

The set of two-boundary time-dependent implicit difference equations is
linearized and solved iteratively using the LINPACK library. Time steps are
not limited because the difference scheme is implicit. Hydrostatic equilibrium
is restored at each time step by solving numerically equation 12. The
corresponding change in gravitational potential energy,
$E_{g}=-\int_{0}^{M(t)}\frac{Gmdm}{r}$ is calculated and we parametrize the
distribution of this energy to obtain the local $\dot{Q}_{grav}$ for the
timestep. The initial and physical parameters assumed for all cases are given
in Table 1.

Table 1: Initial and physical parameters

Parameter | Value
---|---
Initial 40K abundance |  $1.13\times 10^{-6}$ ppm
Initial 232Th abundance |  $5.52\times 10^{-8}$ ppm
Initial 235U abundance |  $6.16\times 10^{-9}$ ppm
Initial 238U abundance |  $2.18\times 10^{-8}$ ppm
Albedo | 0.5
Ice specific energy |  $3.75\times 10^{4}T^{2}+9.0\times 10^{5}T$ erg g-1
Water specific energy |  $4.187\times 10^{7}T$ erg g-1
Rock specific energy |  $1.3\times 10^{7}T$ erg g-1
Water thermal conductivity |  $5.5\times 10^{4}$ erg cm-1 s-1 K-1
Ice thermal conductivity |  $5.67\times 10^{7}/T$ erg cm-1 s-1 K-1
Rock thermal conductivity |  $2\times 10^{5}$ erg cm-1 s-1 K-1
Water dynamic viscosity |  $2.939\times 10^{-4}\exp{\left(\frac{507.88}{T-149.3}\right)}$ dyn cm-2 s
| [ $5.05\exp{(-5.71T/T_{m})}$ for high $T_{m}$]
Latent heat of melting |  $3.34\times 10^{9}$ erg g-1
Latent heat of sublimation |  $A_{0}-A_{1}T+A_{2}T^{2}-A_{3}T^{3}$ erg g-1,
| $A_{0}=3.714\times 10^{10},\ A_{1}=7.823\times 10^{7},$
| $A_{2}=1.761\times 10^{5},\ A_{3}=1.902\times 10^{2}$

Note: The abundances of radioactive nuclides are those of chondritic
meteorites scaled back 4.5 Gyr (Anders & Grevesse, 1989); for ice, the
specific energy was derived by Klinger (1980) based on data from Giauque &
Stout (1936) and the thermal conductivity is from Klinger (1980) as well; for
the thermal properties of rock, see Neumann et al. (2020) and also the recent
paper by Bagheri et al. (2022), based on new data from New Horizons; the
viscosity of ice at high pressure is taken from Pigott et al. (2010); the
thermal conductivity of rock, from Seipold (1998), ignoring the relatively
weak temperature dependence.

##  4 Results of evolutionary calculations

We start by describing the evolution of a prototype model in some detail. We
then turn to discuss the effect of the two leading parameters on the evolution
outcome.

###  4.1 Evolutionary course of a typical case

We choose the case with normalized Safronov parameter $\hat{\theta}=400$ and a
composition of 0.3 ice and 0.7 rock by mass. The initial embryo’s mass is
about 0.00015 of the final mass (roughly 0.1 lunar masses). The temperature is
uniform and equal to the ambient temperature and the structure is in
hydrostatic equilibrium. Accretion starts slowly, but picks up quickly, as
shown in the lower right-panel of Figure 1; the central pressure rises with
the growing mass, as does the bulk density (see upper left-panel of Figure 1).
The accretion rate peaks at about $10^{8}$ yr, declines and ceases at $\sim 1$
Gyr. Accretional heating at the surface, which surpasses the solar energy by
orders of magnitude, causes the surface temperature to rise up to $\sim$180 K.
Although most of this heat is reradiated, a fraction is absorbed in ice
sublimation, which is significant at this temperature (see lower left-panel of
Figure 1). When mass accretion is over, the planet settles into equilibrium,
reradiating the stellar energy absorbed.

Sublimation occurs both at the surface and in a porous subsurface layer, the
vapor flowing to the surface and escaping. The loss of ice affects the bulk
ice to rock ratio, which decreases during the phase of intense heating, but
rises again when the surface temperature drops with the declining accretion
rate, and most of the accreted ice is retained again. This effect is
illustrated in the upper panels of Figure 2, which show the variation of the
bulk ice mass fraction with time (left) and its profile throughout the planet
at the end of evolution (right).

Figure 1: Evolution of various properties for a model with $\hat{\theta}=400$
and $X_{i}=0.3$. Note that time is on a logarithmic scale, emphasizing changes
in evolutionary timescales. Upper left \- the growing planetary mass, the bulk
(average) density, and the central pressure. Upper right \- the maximal
temperature, attained off-center, the central temperature and the surface
temperature. Lower left \- energy per unit time supplied by accretion and
insolation, absorbed by sublimation, and emitted by thermal radiation, where
the small balance is conducted to/out of the interior. Lower right \- the
ratio of energy released by accretional heating and by radioactive decay, and
the accretion rate (for guidance).

Accretional heating of the interior, mainly of the outer layers, surpasses
heating by radioactive decay, therefore the maximal temperature is attained
above the center. The heat wave penetrates inward and when it reaches the core
and the ice reaches the melting temperature (above $10^{3}$ K at these
pressures), the rock settles to the center and separation between ice and rock
begins, forming a rocky core and an increasingly ice-rich mantle. This occurs
at 1.8 Gyr and goes on to the present. The rocky core heats rapidly by
radioactive decay and compression; the onset of core formation is marked by a
rise in central (as well as in maximal) temperature and in central pressure,
as seen in the upper panels of Figure 1. Between the rocky core and the mantle
there is a layer of liquid water mixed with rock, an internal ocean, about 150
km thick, but only a thin layer compared to the planet’s size. This layer
moves out as the core grows.

Figure 2: Characteristics of the model with $\hat{\theta}=400$ and $X_{i}=0.3$
(same as Figure 1): Upper-left \- evolution of the bulk ice to rock ratio,
reflecting the loss of ice at the peak of accretion. Other panels \- profiles
of various characteristics throughout the planet at the end of evolution. In
all profiles the core boundary can be clearly seen. The water line marks the
radial distance below which the temperature exceeds the melting temperature.
Figure 3: Evolution of planetary temperature, density, ice mass faction as a
function of planetary radius and time. The time-frame is chosen to be such
that the planetary radius is constant. $T_{melt}$ is melting temperature of
water. The model assumed is $\hat{\theta}=400$ and ice mass fraction of 0.3.

The continuous post-accretion evolution is shown in Figure 3. We can see the
advance of the surface heat wave toward the center and the marked effect of
radioactive heating of the rock, once ice is separated from it. We also note
that the core forms quickly at 2 Gyr and then expands slowly. At any given
depth, the temperature rises with time.

The emerging structure after 4.5 Gyr of evolution is shown in Figure 2 by
profiles of the main characteristics. The rocky core boundary is clearly seen
by the sharp change in the temperature, density and ice mass fraction profiles
(only the pressure is continuous). The core boundary coincides with the point
where the temperature crosses the melting temperature curve. Although
radioactive heating is stronger in the core, where the rock mass fraction is
higher, the temperature peak occurs away from the center, even after the end
of accretion. This is due to net heat transport outwards by the water, both
because its heat capacity is higher than that of rock, and because it
transports latent heat - absorbing heat by melting in deeper regions and
releasing it by freezing in outer regions. We note that the ice content of the
mantle is much higher than the original. The core extends to about 2/3 of the
planet’s radius.

We show in what follows that the extent of the core depends critically on the
parameters that we investigate, to the point that it does not form at all. The
history of ice accretion remains imprinted in the composition of the planet,
where the ice content has a minimum at roughly half the mass at the end of the
accretion phase. This memory is eventually erased when the core forms, but it
remains present in models that do not form a core.

###  4.2 Parameter Study

####  4.2.1 Evolution

We now choose a baseline case that has a normalized Safronov parameter
$\hat{\theta}_{b}=400$ and an ice mass fraction $X_{i,b}=0.5$ (ice to rock=1).
To determine the effect of the accretion rate and correspondingly, of the
accretion heating, we keep $X_{i,b}$ and vary $\hat{\theta}$ between 200 and
500. Similarly, to determine the effect of initial composition we keep
$\hat{\theta}_{b}$ and vary the ice mass fraction between 0.3 and 0.7. All
models start with a very small embryo ($\sim$0.00015M, where M is final mass)
having the same composition as the accreted material.

In Figure 4 we show evolution results for $\hat{\theta}_{b}$, and varying ice
mass fraction. In this set of models, the final mass is found to be almost
independent of the planetesimal composition, while other properties depend on
it strongly. As expected, the radius gradually increases with the ice mass
fraction, while the bulk density $\bar{\rho}$ (fifth panel, Figure 4) and
central pressure $P_{c}$ decrease. The mass accretion rate $\dot{M}$,
increases with the ice mass fraction, because the radius of the planet
increases (see equation 5). In addition, the peak becomes somewhat narrower
for higher ice mass fractions because the larger planetary radius allows it to
empty out the feeding zone more quickly.

Several interesting effects can be seen. In the seventh panel showing the ice
to rock ratio as a function of time, all models exhibit the same pattern:
initially this ratio is equal to that of the accreted material, but around the
time of maximum $\dot{M}$, as the accretion heating causes ice to evaporate,
the ratio is lowered. After $\dot{M}$ decreases again, the loss of water
decreases as well, but since a large fraction of the planet has been accreted
by this time, the ice to rock ratio remains visibly below its original value.
The amount of mass lost by evaporation is roughly up to 10% of the final mass
of the planet. The peak in $\dot{M}$ corresponds to the peak in the surface
temperature (dashed curves in the last panel).

Figure 4: Various planetary properties versus time for different ice mass
fraction. The models assume $\hat{\theta}=400$. Top row panels: the planetary
mass $M(t)$; the mass of the escaping vapor $M_{\rm esc}(t)$; the mass
accretion rate $\dot{M}(t)$; the planetary radius $R(t)$. Bottom row panels:
the average (bulk) density $\bar{\rho}(t)$; the central pressure $P_{c}(t)$,
the ice-to-rock mass ratio; the central and surface temperatures, $T_{c}(t)$
and $T_{s}(t)$ respectively. Figure 5: Various planetary properties versus
time - as in Fig.4 for different values of $\hat{\theta}$. The models assume
an ice mass fraction of 0.5.

Another effect can be seen in the plots of the central pressure $P_{c}$ and
central temperature (solid curves in the last panel). There is a sharp change
of slope in the curves for the lower two values of the ice mass fraction. This
is due to radioactive heating after core formation, as mentioned in the
previous section. The effect is not seen in models of high ice content. A high
ice content has two consequences, which combined, prevent the formation of a
rocky (ice-depleted) core. First, the heat capacity of ice is higher than that
of rock and secondly, the relatively low content of rock means lower
abundances of radioactive species and hence less radioactive heating. As a
result, the internal temperatures do not rise sufficiently for the ice to melt
and for rock to separate from the ice. The core formation causes a small
decrease in the mean density of the body, and this is also visible in the plot
of bulk density as a function of time.

In Figure 5 we show results for models where the ice mass fraction is kept at
$X_{i,b}$, while $\hat{\theta}$ gradually changes from 200 to 500. As can be
seen from Equation 5, $\dot{M}$ will change by approximately the same factor.
In addition, when $\dot{M}$ is larger, the available mass is accreted more
efficiently, and the accretion heating is more concentrated in time. At the
highest $\hat{\theta}$, nearly 20% of the total mass is lost to water
evaporation and as a result, the final ice to rock ratio is just over 70% of
its original value. For the lowest value of $\hat{\theta}$, the effect is
barely noticeable.

Only the model corresponding to the highest $\hat{\theta}$ value forms a large
core, and the core starts forming relatively late. We have seen that bodies
with ice mass fractions lower than 0.5 are more likely to form cores. In the
present case, ice loss during accretion lowered the ice content sufficiently
for a core to form.

####  4.2.2 Final structure - after 4.5 Gyr

Figure 6 shows the planetary internal structure after $4.5\times 10^{9}$ years
of accretion and evolution. Profiles are shown for $\hat{\theta}_{b}$ and
different ice mass fractions. For ice mass fractions of 0.5 and above, the
density is smoothly decreasing as a function of radius. For lower ice mass
fractions, where there is sufficient radioactive heating from the rock, so
that melting occurs and some of the rock settles towards the center to form a
core, the ice to rock ratio is zero in the central part. The extent of the
rocky core is larger for the lower ice content. In bodies that do not form a
core, one can see in the ice to rock profile a depression corresponding to the
accretion phase where ice loss was strongest. The central pressure is
inversely proportional to the ice mass fraction, since rock is more
compressible than ice.

Figure 6: Various planetary properties versus radius in the final stage of
planetary evolution, for different compositions. The models assume
$\hat{\theta}=400$.  Figure 7: Various planetary properties versus radius in
the final stage of planetary evolution, for different normalized Safronov
parameters. The models assume 0.5 ice.

Figure 7 shows how the internal structures vary for different values of
$\hat{\theta}$ for the initial ice mass fraction $X_{i,b}$. For the lower
values of $\hat{\theta}$, the internal structure is quite similar. The effect
of lower peak accretion rate, resulting in weaker accretional heating and less
sublimation, is seen in the ice to rock profiles. Recall that the 0.5 ice case
does not have enough radioactive heating to form a distinct core for
$\hat{\theta}<400$. But for $\hat{\theta}=500$, and $\hat{\theta}=400$ we
clearly see the distinct density discontinuity at near $0.8R_{\oplus}$ and
$0.4R_{\oplus}$, respectively, marking the core-mantle boundary. This is also
reflected in the ice to rock ratio and the temperature profile.

####  4.2.3 Combined effect of parameters

We now present the full picture, resulting from the combined effect of the
accretion rate—through the Safronov parameter—and the composition—through the
initial ice to rock ratio, represented by a set of models corresponding to
three values of the normalized Safronov parameter ($\hat{\theta}=200,400,500$)
and three values of the initial ice mass fraction (0.3, 0.5, 0.7) and all 9
combinations between them.

Figure 8 shows the combined effect of varying the water ice mass fraction and
$\hat{\theta}$ simultaneously, while Figure 9 does the same for the various
parameter profiles after 4.5 Gyr of evolution. The figures show that a high
$\dot{M}$ can help initiate the process of core formation if the rock mass
fraction is marginally too low. Another signature of the accretion rate is
seen in the dip in the ice to rock ratio around the radius corresponding to
the peak in $\dot{M}$.

Figure 8: Various planetary properties versus time during planetary evolution,
for different compositions and Safronov parameters.  Figure 9: Various
planetary properties versus radius in the last stage of planetary evolution,
for different compositions and $\hat{\theta}$.

###  4.3 Effect of fixed parameters

####  4.3.1 Results for a smaller planetary mass

In all of the cases presented, the final mass of the planet was almost the
same. This is because that mass is essentially determined by the solid mass
available in the disk via Equation 3. In order to gauge the effect of
$\dot{M}$ on a smaller mass body, we ran our model for the case of
$\Sigma_{s}=0.1$ $\rm g~{}cm^{-2}$. This has the effect of reducing the final
mass of the planet by half.

Figure 10: Various planetary properties versus time during planetary
evolution, for different ice fractions. The isolation mass is half of
original. Can be compared with Figure 4. Figure 11: Various planetary
properties versus radius in the last stage of planetary evolution, for
different compositions, for $\hat{\theta}=400$. The isolation mass is half of
original.

The results are shown in Figures 10 and 11. We note that although for these
lower-mass models, the densities and pressures are much lower than those of
the more massive ones, the temperature profiles are similar, because the
radiogenic heat source is the same per unit mass. Similarly, the occurrence
and extent of the rocky core are the same. The marked difference is that
surface heating by accretion is not sufficient for ice loss, and thus the ice
to rock ratio does not differ from the initial value.

####  4.3.2 Accretional heating

It is often assumed that accretional heating should raise the interior
temperatures of planets to very high values (Zeng & Sasselov, 2014; Boujibar
et al., 2020), higher than those reached by our models. This implies that the
accretion energy (or most of it) has been absorbed by the object during
formation. We argue that this assumption is not realistic. In order for the
accretion energy, which is essentially a surface source, to be turned into
internal thermal energy, the thermal diffusion timescale $R/\kappa$, where
$\kappa$ is the diffusivity, should be comparable to the accretion timescale
$M/\dot{M}$. But, for the accretion rate used in Section 4.1, corresponding to
$\hat{\theta}=400$, the diffusion timescale is over four orders of magnitude
larger than the accretion timescale, and the discrepancy increases for more
rapid accretion.

The question then is, where does the accretion energy go? Some of it goes into
ice sublimation (if there is ice in the accreted material), but there is not
enough ice to absorb all of it. Most of the energy is simply reradiated. The
surface temperature adjusts itself to radiate it, emitting accretion
luminosity. If the rate of energy release by accretion is $\alpha
GM\dot{M}/R$, where $\alpha$ is of order unity, we have

$\alpha\frac{GM\dot{M}}{R}=4\pi R^{2}\sigma T_{s}^{4},$ (21)

which yields

$T_{s}=\left(\frac{\alpha G\bar{\rho}}{3\sigma}\right)^{1/4}\dot{M}^{1/4}.$
(22)

Thus, $T_{s}$ only needs to be of the order of a few hundred Kelvins in order
to accomplish this. If there is ice close to the surface it is controlled by
the ice, in which case $\sigma T_{s}^{4}$ on the right-hand side of Equation
21 is replaced by $Z(T_{s})H_{v}$. Our detailed results, shown in Figure 1,
confirm this estimate.

Large gas planets that form by collapse are different: there, the energy
source is a body source that heats the interior. Even before collapse, if the
core grows quickly enough, it can attract a hydrogen-helium envelope from the
surrounding gas disk before that disk dissipates (e.g. Pollack et al., 1996).
In that case the heating comes from the accretional energy of infalling
planetesimals and that energy will be deposited inside a massive envelope. Any
water vaporization at the core surface due to this accretional energy will be
mixed with the H-He envelope and provide heat throughout a large volume, not
just at the surface.

Changing the value of $f$—the division of accretion energy between surface and
interior—has the expected effect, but does not change our conclusions: a
higher $f$ results in more ice loss upon accretion, and converges to a case of
an initially lower ice to rock ratio, while a lower $f$, which leaves more ice
in the planet, results in lower internal temperatures, similar to the case of
an initially higher ice to rock ratio. In both cases, the surface temperature
changes only slightly.

##  5 Discussion and Conclusions

The formation of a planet, assuming that it is a continuous process, rather
than a cataclysmic event (or a chain of such events), is determined by two
major parameters: the accretion rate and the composition of the accreted
material. These are the parameters considered in the present study, assuming
the composition to be a simple mixture of water ice and silicate rock, and
considering continuous spherical accretion. Beside understanding evolutionary
processes of such planets, our goal is to link the interior structure to the
planet’s birth environment, as well as to constrain the possible interior
structure based on measurements of mass and radius - by no means a
straightforward task (e.g. Valencia et al., 2006; Rogers & Seager, 2010; Dorn
et al., 2015), because different formation scenarios can lead to similar
values for the same observables (e.g. Mordasini et al., 2012; Lozovsky et al.,
2018).

For this purpose, we have run a large number of models for different
combinations of the two leading parameters. The accretion rate is determined
by the initial protoplanetary disc structure, in particular its density. Since
in any scenario, the accretion rate is proportional to the local density of
material and to the growing planetary embryo, the former decreasing with time
and the latter, increasing, the rate is time dependent and has a maximum. We
find that the functional shape of $\dot{M}(t)$ has a crucial effect on the
final structure of the planet.

It is important to note that different assumptions regarding the mechanism of
formation and accretion (such as pebble accretion) lead to very different
$\dot{M}$ profiles (e.g. Ida et al., 2016; Lambrechts & Johansen, 2012; Bitsch
et al., 2015; Schneider & Bitsch, 2021; Alibert et al., 2018; Grishin et al.,
2020). These studies assume that the growth is rapid enough so that the core
can attract a significant H-He envelope from the surrounding nebula. It was
found, using the pebble accretion paradigm, that sub-Neptune sized planet
formation is preferable for relatively low solid accretion rates ($\sim
10^{-6}M_{\oplus}$ yr-1 ), which corresponds to low-metallicity environments
(Venturini & Helled, 2017). However, some claim that pebble accretion plays
little or no role at all in terrestrial planet formation (Mah et al., 2021) or
in planets formed more than a few au from the central star (Voelkel et al.,
2021a, b, 2022).

We also note that migration during the formation process, which we have not
taken into account, can also affect the accretion rate, as well as the
composition of accreted material (Alibert et al., 2005; Shibata et al., 2020;
Lega et al., 2021; Bitsch et al., 2019). In the scenario we investigate,
$\dot{M}$ is small enough so that the nebular gas dissipates well before a
significant core can grow and therefore migration is not important in our
case.

Studies of young stellar clusters indicate that the protostellar gaseous disk
dissipates after $\sim$10 Myr (Haisch et al., 2001). At this time, for the
accretion rates we assume, the protoplanetary embryo has a mass less than
$0.05M_{\oplus}$. At 40 au such a body will have a Bondi radius of $\sim
2\times 10^{10}$ cm. For the nebular model we are considering, a sphere of
this radius would contain $2.6\times 10^{18}\ {\rm g}=4.4\times
10^{-10}M_{\oplus}$ of nebular gas. Even if all this gas was accreted onto the
embryo, it would comprise a minuscule fraction of the total planet, and
therefore we have neglected H-He accretion.

An important effect of the accretion rate is due to accretion heating, which
is proportional to it. We have shown that if heating is intense enough to
raise the surface temperature to the sublimation temperature of water ice,
then a fraction of the accreted ice escapes ($M_{\rm esc}$). This effect was
also found by Bierson & Nimmo (2020). Therefore the ice to rock ratio of a
planet does not necessarily reflect the composition of the environment where
it was formed. The fate of the escaping vapor into the cold environment is not
clear and is not included in our planetary model. While some of it may
condense on dust grains, or else drift to other regions of the disc, some may
form a temporary water vapor atmosphere.

It is important to note that the water vapor escapes while the planet grows
and therefore it is in constant interaction with the infalling planetesimals.
Several studies have shown that planetesimal impacts lead to atmospheric mass
loss (Schlichting et al., 2015; Biersteker & Schlichting, 2021; Denman et al.,
2020, 2022; Wyatt et al., 2020). Significant loss may be caused by large
impactors (Denman et al., 2022), but the cumulative effect of a large number
of small and particularly medium-size impactors may lead to the loss of the
atmosphere as well (Schlichting et al., 2015). Assuming a discrete power law
size distribution of impacting planetesimals (radii $a_{i}$) with the commonly
assumed power index $\gamma=-3.5$ (Wyatt et al., 2020), a mass range of 1 km
to a few 1000 km and a total mass equal to the planetary mass $M$, the number
of impactors with radius $a_{i}$ is

$N_{i}=\frac{3M}{4\pi\rho_{p}}\frac{a_{i}^{-\gamma}}{\sum_{i}a_{i}^{3-\gamma}},$
(23)

where $\rho_{p}$ is the density of a planetesimal. The total number of impacts
is of the order of $10^{12}$, the number of large impactors ($a_{i}>1000$ km)
is on the order of 10. Extending the sizes distribution to lower radii reduces
this number, while a shallower size distribution, say, $\gamma=-2.5$,
increases the number of large impactors by more than a factor of 10. We note
that the energy required for the vapor to escape is accounted for in the
calculation, since only a fraction of the nominal accretion energy is assumed
to be absorbed by the growing planet (see Section 3), thus in this respect our
assumptions are self-consistent. A thorough calculation of the effect of
impacting planetesimals is beyond the scope of the present study and will be
considered in future work, but it seems reasonable to assume that the
atmosphere or most of it should be stripped away. We have also tested other
processes that can dissipate the atmosphere, such as Jeans escape,
hydrodynamic escape or photoevaporation, but found their effect negligible at
this distance from the central star. In any case, a detailed atmospheric model
is necessary to follow the temperature profile and to compute the structure
and thermal evolution of such a temporary atmosphere, and to estimate its
escape rate.

Although we have not considered additional volatiles in our study, they will
be influenced by the same physics as the water. Öberg & Wordsworth (2019) and
Bosman et al. (2019) have suggested that the uniform enrichment of volatiles
in Jupiter’s envelope could be due to formation at the outer edge of the solar
system where such volatiles could be accreted as ices. This scenario needs to
be revisited taking surface heating into account to see if these ices are
indeed retained.

We find that the ice to rock ratio in the environment where the planet grows,
determines whether or not a rocky core will form, by rock settling and
separation from the ice. This is induced by radiogenic heating by the common
long-lived radioisotopes assumed to be embedded in the rock. If the rock
content is less than 0.5, the temperature does not rise sufficiently for the
ice to melt and allow rock settling, keeping in mind that the melting
temperatures are above $10^{3}$ K at the internal pressures of the planet. The
result is an almost homogeneous configuration. Even when a core forms, the
displaced ice increases the ice to rock ratio in the mantle, but there is no
pure ice layer. Between the rocky core and the ice-rich mantle, a relatively
thin layer is found, where superheated water is mixed with the rock. The
temperature is however too high for the water to be chemically absorbed by the
rock.

Figure 12: Left: Bulk density of the final structure of models with
$\hat{\theta}=400$ as a function of the initial ice to rock ratio. Red points
mark models that have developed a rocky core. Right: Relative change of the
ice to rock ratio due to loss of ice during accretion, as a function of the
initial ratio.

In principle, if the mass of a planet and its radius are known, so that the
bulk density can be derived, models like those presented here may be used to
derive the planet’s ice to rock ratio and to determine whether the planet
possesses a rocky core. Furthermore, the same models will reveal what was the
ice to rock ratio in the environment where the planet formed. This is
illustrated in Figure 12, which shows—for a given Safronov parameter—the
dependence of bulk density and the relative change in ice to rock ratio as a
function of initial ice to rock ratio.

Figure 13:  Mass - Radius relation for models with various ice mass fractions
and $\hat{\theta}=400$. Circles denote the sample of discovered exoplanets, up
to 8 $M_{\oplus}$ (based on NASA Exoplanet Archive (2022)). See text for
discussion.

Our results may be seen in the light of exoplanets. In Figure 13 we show the
$R-M$ relation for the baseline model on the background of observed
exoplanets. While the scatter of exoplanets in the $R-M$ plane is
considerable, there are numerous objects in the region of our R(M) relation,
indicating an ice-rich composition (in the absence of a H-atmosphere).
Admittedly, the observed exoplanets are much closer to the host star than our
models, but if these bodies were formed in debris disks after the nebular gas
has dissipated, the accretion mechanism could be similar to the one we
consider. In addition, the orbital distance of 40 au used in our study was
coupled with the assumption of a solar-like central star; the same boundary
condition may be obtained much closer to a fainter star. Moreover, the stellar
radiation contributes only a small fraction of the energy for a protoplanet,
as shown in Figure 1, and therefore, the formation distance may be smaller
still. There is however a lower limit, imposed by the snow line (demanding
that the local equilibrium temperature be lower than 170-180 K). Planets that
are now closer to the host star, but have densities indicating an ice-rich
composition, may have migrated inward since formation (Bean et al., 2021, and
references therein). For now, the snow line is beyond the detection limit of
small planets, but we may expect more ice-rich super-earths to be discovered
in future.

A very interesting study of the Trappist-1 system by Coleman et al. (2019)
investigated how the ice content of these planets differed for planets formed
via pebble accretion as compared to planets formed via planetesimal accretion.
This study should be revisited taking the physics of evaporation at the
surface into account. More careful treatment of the details of the accretion
process can shed new light on the modeling of exoplanet evolution.

Finally, one of the interesting results of this work is the dependence of the
final composition on the Safronov parameter (see e.g., Figure 9)—for
relatively high values of $\hat{\theta}$, the ice mass fraction in the planets
is lower than the assumed ice mass fraction of the planetesimals. However, one
should keep in mind that the Safronov parameter indicates the efficiency of
accretion and might depend on the mass spectrum (e.g. Namouni et al., 1996).
Therefore, the ice to rock ratio of the accreted material and the solid
surface density of the disk might influence $\hat{\theta}$. Hence it is
possible that the two parameters considered here as independent, may be
linked. This would make the present study much more conclusive by eliminating
some of the parameter combinations.

We acknowledge the support of the Israeli Science Foundation for this research
through the ISF grant 566/17.

##  Appendix A List of models

ice mass fraction | $\hat{\theta}$ | Fig.4, 6 | Fig. 5, 7 | Fig. 8, 9 | Fig. 10, 11 | M ($M_{\oplus}$) | R ($R_{\oplus}$)
---|---|---|---|---|---|---|---
Full Mass
0.3 | 200 |  |  | x |  | 8.29 | 2.30
0.5 | 200 |  | x | x |  | 8.90 | 2.48
0.7 | 200 |  |  | x |  | 8.91 | 2.61
0.5 | 300 |  | x |  |  | 8.82 | 2.48
0.3 | 400 | x |  | x |  | 8.29 | 2.30
0.4 | 400 | x |  |  |  | 8.32 | 2.36
0.5 | 400 | x | x | x |  | 8.37 | 2.43
0.6 | 400 | x |  |  |  | 8.43 | 2.50
0.7 | 400 | x |  | x |  | 8.49 | 2.58
0.3 | 500 |  |  | x |  | 7.82 | 2.24
0.5 | 500 |  | x | x |  | 7.69 | 2.36
0.7 | 500 |  |  | x |  | 7.92 | 2.52
Half Mass
0.3 | 400 |  |  |  | x | 4.45 | 2.01
0.4 | 400 |  |  |  | x | 4.45 | 2.06
0.5 | 400 |  |  |  | x | 4.45 | 2.11
0.6 | 400 |  |  |  | x | 4.45 | 2.17
0.7 | 400 |  |  |  | x | 4.45 | 2.23
Table 2: List of models used in this study and figures there they presented. M
and R are the final mass and radius of the bodies. "Half Mass" is referring to
models with a isolation mass set to be half of the original ("Full Mass").

In Table 2.A we list the models used in this study, with the corresponding
figures numbers.

## References

  * Adibekyan et al. (2021) Adibekyan, V., Dorn, C., Sousa, S. G., et al. 2021, Science, 374, 330, doi: 10.1126/science.abg8794
  * Alibert et al. (2005) Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343, doi: 10.1051/0004-6361:20042032
  * Alibert et al. (2018) Alibert, Y., Venturini, J., Helled, R., et al. 2018, Nature Astronomy, 2, 873, doi: 10.1038/s41550-018-0557-2
  * Andama et al. (2021) Andama, G., Ndugu, N., Anguma, S. K., & Jurua, E. 2021, Monthly Notices of the Royal Astronomical Society, 510, 1298, doi: 10.1093/mnras/stab3508
  * Anders & Grevesse (1989) Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197, doi: 10.1016/0016-7037(89)90286-X
  * Andrews (2020) Andrews, S. M. 2020, Annual Review of Astronomy and Astrophysics, 58, 483, doi: 10.1146/annurev-astro-031220-010302
  * Armitage (2020) Armitage, P. J. 2020, Astrophysics of Planet Formation, 2nd edn. (Cambridge University Press), doi: 10.1017/9781108344227
  * Bagheri et al. (2022) Bagheri, A., Khan, A., Deschamps, F., et al. 2022, Icarus, 376, 114871, doi: 10.1016/j.icarus.2021.114871
  * Batalha et al. (2013) Batalha, N. M., Rowe, J. F., Bryson, S. T., et al. 2013, ApJS, 204, 24, doi: 10.1088/0067-0049/204/2/24
  * Bean et al. (2021) Bean, J. L., Raymond, S. N., & Owen, J. E. 2021, Journal of Geophysical Research (Planets), 126, e06639, doi: 10.1029/2020JE006639
  * Bierson & Nimmo (2020) Bierson, C. J., & Nimmo, F. 2020, ApJ, 897, L43, doi: 10.3847/2041-8213/aba11a
  * Biersteker & Schlichting (2021) Biersteker, J. B., & Schlichting, H. E. 2021, MNRAS, 501, 587, doi: 10.1093/mnras/staa3614
  * Bitsch et al. (2015) Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112, doi: 10.1051/0004-6361/201526463
  * Bitsch et al. (2019) Bitsch, B., Raymond, S. N., & Izidoro, A. 2019, A&A, 624, A109, doi: 10.1051/0004-6361/201935007
  * Bodenheimer & Pollack (1986) Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391, doi: 10.1016/0019-1035(86)90122-3
  * Bosman et al. (2019) Bosman, A. D., Cridland, A. J., & Miguel, Y. 2019, A&A, 632, L11, doi: 10.1051/0004-6361/201936827
  * Boujibar et al. (2020) Boujibar, A., Driscoll, P., & Fei, Y. 2020, Journal of Geophysical Research (Planets), 125, e06124, doi: 10.1029/2019JE006124
  * Coleman et al. (2019) Coleman, G. A. L., Leleu, A., Alibert, Y., & Benz, W. 2019, A&A, 631, A7, doi: 10.1051/0004-6361/201935922
  * Denman et al. (2022) Denman, T. R., Leinhardt, Z. M., & Carter, P. J. 2022, MNRAS, doi: 10.1093/mnras/stac923
  * Denman et al. (2020) Denman, T. R., Leinhardt, Z. M., Carter, P. J., & Mordasini, C. 2020, MNRAS, 496, 1166, doi: 10.1093/mnras/staa1623
  * Dorn et al. (2015) Dorn, C., Khan, A., Heng, K., et al. 2015, A&A, 577, A83, doi: 10.1051/0004-6361/201424915
  * Dorn & Lichtenberg (2021) Dorn, C., & Lichtenberg, T. 2021, The Astrophysical Journal Letters, 922, L4, doi: 10.3847/2041-8213/ac33af
  * Frelikh & Murray-Clay (2017) Frelikh, R., & Murray-Clay, R. A. 2017, AJ, 154, 98, doi: 10.3847/1538-3881/aa81c7
  * Fressin et al. (2013) Fressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766, 81, doi: 10.1088/0004-637X/766/2/81
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, The Astronomical Journal, 154, 109.  http://stacks.iop.org/1538-3881/154/i=3/a=109
  * Giauque & Stout (1936) Giauque, W. F., & Stout, J. W. 1936, Journal of the American Chemical Society, 58, 1144, doi: 10.1021/ja01298a023
  * Grishin et al. (2020) Grishin, E., Rozner, M., & Perets, H. B. 2020, The Astrophysical Journal, 898, L13, doi: 10.3847/2041-8213/aba266
  * Haisch et al. (2001) Haisch, Karl E., J., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153, doi: 10.1086/320685
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35, doi: 10.1143/PTPS.70.35
  * Helled et al. (2014) Helled, R., Bodenheimer, P., Podolak, M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 643, doi: 10.2458/azu_uapress_9780816531240-ch028
  * Hou & Wei (2022) Hou, Q., & Wei, X. 2022, arXiv e-prints, arXiv:2201.07008.  https://arxiv.org/abs/2201.07008
  * Ida et al. (2016) Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72, doi: 10.1051/0004-6361/201628099
  * Klinger (1980) Klinger, J. 1980, Science, 209, 271, doi: 10.1126/science.209.4453.271
  * Lambrechts & Johansen (2012) Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32, doi: 10.1051/0004-6361/201219127
  * Laughlin (2018) Laughlin, G. 2018, Mass-Radius Relations of Giant Planets: The Radius Anomaly and Interior Models, ed. H. J. Deeg & J. A. Belmonte (Cham: Springer International Publishing), 1357–1373, doi: 10.1007/978-3-319-55333-7_1
  * Lega et al. (2021) Lega, E., Nelson, R. P., Morbidelli, A., et al. 2021, A&A, 646, A166, doi: 10.1051/0004-6361/202039520
  * Lissauer (1987) Lissauer, J. J. 1987, Icarus, 69, 249, doi: 10.1016/0019-1035(87)90104-7
  * Lodders (2010) Lodders, K. 2010, Astrophysics and Space Science Proceedings, 16, 379, doi: 10.1007/978-3-642-10352-0_8
  * Lozovsky et al. (2018) Lozovsky, M., Helled, R., Dorn, C., & Venturini, J. 2018, The Astrophysical Journal, 866, 49.  http://stacks.iop.org/0004-637X/866/i=1/a=49
  * Lozovsky et al. (2021) Lozovsky, M., Helled, R., Pascucci, I., et al. 2021, A&A, 652, A110, doi: 10.1051/0004-6361/202140563
  * Lozovsky et al. (2017) Lozovsky, M., Helled, R., Rosenberg, E. D., & Bodenheimer, P. 2017, The Astrophysical Journal, 836, 227, doi: 10.3847/1538-4357/836/2/227
  * Mah et al. (2021) Mah, J., Brasser, R., Bouvier, A., & Mojzsis, S. J. 2021, Monthly Notices of the Royal Astronomical Society, doi: 10.1093/mnras/stab3766
  * Mordasini et al. (2012) Mordasini, C., Alibert, Y., Georgy, C., et al. 2012, A&A, 547, A112, doi: 10.1051/0004-6361/201118464
  * More et al. (1988) More, R. M., Warren, D. A., Young, D. A., & Zimmerman, G. B. 1988, Physics of Fluids, 31, 3059
  * Namouni et al. (1996) Namouni, F., Luciani, J. F., & Pellat, R. 1996, A&A, 307, 972
  * NASA Exoplanet Archive (2022) NASA Exoplanet Archive. 2022, Confirmed Planets Table, IPAC, doi: 10.26133/NEA1
  * Neumann et al. (2020) Neumann, W., Jaumann, R., Castillo-Rogez, J., Raymond, C. A., & Russell, C. T. 2020, A&A, 633, A117, doi: 10.1051/0004-6361/201936607
  * Öberg & Wordsworth (2019) Öberg, K. I., & Wordsworth, R. 2019, AJ, 158, 194, doi: 10.3847/1538-3881/ab46a8
  * Ogihara et al. (2022) Ogihara, M., Kokubo, E., Nakano, R., & Suzuki, T. K. 2022, arXiv e-prints, arXiv:2201.08840, doi: 10.1051/0004-6361/202142354
  * Ormel et al. (2021) Ormel, C. W., Vazan, A., & Brouwers, M. G. 2021, A&A, 647, A175, doi: 10.1051/0004-6361/202039706
  * Pigott et al. (2010) Pigott, J. S., Reaman, D. M., & Panero, W. R. 2010, in AGU Fall Meeting Abstracts, Vol. 2010, DI51B–1865
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62, doi: https://doi.org/10.1006/icar.1996.0190
  * Pollack et al. (1986) Pollack, J. B., Podolak, M., Bodenheimer, P., & Christofferson, B. 1986, Icarus, 67, 409, doi: 10.1016/0019-1035(86)90123-5
  * Prialnik & Merk (2008) Prialnik, D., & Merk, R. 2008, Icarus, 197, 211, doi: 10.1016/j.icarus.2008.03.024
  * Ransford (1982) Ransford, G. A. 1982, Physics of the Earth and Planetary Interiors, 29, 209, doi: 10.1016/0031-9201(82)90012-7
  * Rogers & Seager (2010) Rogers, L. A., & Seager, S. 2010, ApJ, 712, 974, doi: 10.1088/0004-637X/712/2/974
  * Safronov (1972) Safronov, V. 1972, Israel program for scientific translations
  * Schlichting et al. (2015) Schlichting, H. E., Sari, R., & Yalinewich, A. 2015, Icarus, 247, 81, doi: 10.1016/j.icarus.2014.09.053
  * Schneider & Bitsch (2021) Schneider, A. D., & Bitsch, B. 2021, A&A, 654, A71, doi: 10.1051/0004-6361/202039640
  * Seipold (1998) Seipold, U. 1998, Tectonophysics, 291, 161, doi: 10.1016/S0040-1951(98)00037-7
  * Shibata et al. (2020) Shibata, S., Helled, R., & Ikoma, M. 2020, A&A, 633, A33, doi: 10.1051/0004-6361/201936700
  * Silva & Valio (2016) Silva, D., & Valio, A. 2016, Proceedings of the International Astronomical Union, 12, 301–304, doi: 10.1017/S1743921317004501
  * Squyres et al. (1988) Squyres, S. W., Reynolds, R. T., Summers, A. L., & Shung, F. 1988, J. Geophys. Res., 93, 8779, doi: 10.1029/JB093iB08p08779
  * Stevenson (2013) Stevenson, D. J. 2013, Gas and Ice Giant Interiors, ed. T. D. Oswalt, L. M. French, & P. Kalas (Dordrecht: Springer Netherlands), 195–221, doi: 10.1007/978-94-007-5606-9_4
  * Valencia et al. (2006) Valencia, D., O’Connell, R. J., & Sasselov, D. 2006, Icarus, 181, 545 , doi: https://doi.org/10.1016/j.icarus.2005.11.021
  * Valletta & Helled (2020) Valletta, C., & Helled, R. 2020, ApJ, 900, 133, doi: 10.3847/1538-4357/aba904
  * Vazan et al. (2015) Vazan, A., Helled, R., Kovetz, A., & Podolak, M. 2015, ApJ, 803, 32, doi: 10.1088/0004-637X/803/1/32
  * Vazan et al. (2013) Vazan, A., Kovetz, A., Podolak, M., & Helled, R. 2013, MNRAS, 434, 3283, doi: 10.1093/mnras/stt1248
  * Vazan et al. (2018) Vazan, A., Ormel, C. W., Noack, L., & Dominik, C. 2018, ApJ, 869, 163, doi: 10.3847/1538-4357/aaef33
  * Vazan et al. (2020) Vazan, A., Sari, R., & Kessel, R. 2020, arXiv e-prints, arXiv:2011.00602.  https://arxiv.org/abs/2011.00602
  * Venturini & Helled (2017) Venturini, J., & Helled, R. 2017, The Astrophysical Journal, 848, 95, doi: 10.3847/1538-4357/aa8cd0
  * Villenave et al. (2021) Villenave, M., Ménard, F., Dent, W. R. F., et al. 2021, A&A, 653, A46, doi: 10.1051/0004-6361/202140496
  * Voelkel et al. (2021a) Voelkel, O., Deienno, R., Kretke, K., & Klahr, H. 2021a, A&A, 645, A131, doi: 10.1051/0004-6361/202039214
  * Voelkel et al. (2021b) —. 2021b, A&A, 645, A132, doi: 10.1051/0004-6361/202039245
  * Voelkel et al. (2022) Voelkel, O., Klahr, H., Mordasini, C., & Emsenhuber, A. 2022, Astronomy & Astrophysics, doi: 10.1051/0004-6361/202141830
  * Wetherill & Cox (1985) Wetherill, G. W., & Cox, L. P. 1985, Icarus, 63, 290, doi: 10.1016/0019-1035(85)90012-0
  * Wyatt et al. (2020) Wyatt, M. C., Kral, Q., & Sinclair, C. A. 2020, MNRAS, 491, 782, doi: 10.1093/mnras/stz3052
  * Zeng & Sasselov (2014) Zeng, L., & Sasselov, D. 2014, ApJ, 784, 96, doi: 10.1088/0004-637X/784/2/96