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

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.


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 ⊕ , 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 (Lo-zovsky 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.

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) Σ g (a) = Σ 0 a −3/2 (1) where a is the radial distance in au, and Σ 0 is the value of Σ g at a =1 au. If we assume that the disk extends to 50 au and has a mass of M d = 0.02M , then we have where a 1 and a 2 are the inner and outer radii of the disk, respectively. Setting M d = 4 × 10 31 g and taking a 1 = 0.01 au and a 2 = 50 au in Equation (2) gives Σ 0 = 2 × 10 3 g cm −2 . The value of Σ 0 is insensitive to the choice of a 1 so long as a 1 a 2 .
We consider the region at ∼ 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, Σ s ≈ 0.16 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 ∆a on either side, where ∆a is related to the Hill's sphere radius of the final body. This isolation mass is given by Lissauer (1987) as where, again, a is in au, B ∼ 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 × 10 28 g or about 7 M ⊕ .

Mass Accretion Rate
If R is the body's radius, ρ 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) 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) where Θ = v 2 e /2v 2 is the so-called Safronov parameter, and Ω is the Kepler frequency. Lissauer (1987), based on a generalization of the results of Wetherill & Cox (1985), gives whereρ is the mean density of the planetary embryo in cgs, and a is again in au. If we ignore the very weak dependence onρ, we find that at 40 au the expected value is Θ = 16, 000. We assume that the total mass (planetesimals + growing embryo) in this annulus remains constant with time. In what follows, we allow Θ 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Ṁ 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θ = Θ a , than a is in au units.

Composition
As the planet grows in a region of ∼ 40 au, we assume that the planetesimals are composed of rock (SiO 2 ) and water ice (H 2 O). We do not take into account H-He gas, since the prescription we have chosen forṀ 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.

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θ, 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 = X i /[1 + e β(1−T /Tm) ] (see Prialnik & Merk 2008, for further details). When the ice is heated and eventually melts, its viscosity η(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 J r and an equal upward flux of ice/water 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 = 2 9 (ρ r − ρ i )gr 2 p /η, (Stokes flow) where ρ r and ρ 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Ṁ /R 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 1 4 of the kinetic energy obtained if the accreted material free-falls from infinity (c.f. 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, where ρu denotes the sum of weighted specific energies, F is the heat flux, F = −K∇T , K denoting the thermal conductivity, andQ 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 masṡ M (t) given by Equation 4, the surface boundary condition becomes where σ 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 H v is the corresponding latent heat absorbed. For mass, At the surface, J i = Z(T ). Finally, the bulk density profile at any time t is given by the hydrostatic equation, where P (ρ) 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 where ρ i and ρ r are the ice and rock densities, respectively, for each of which we adopt the secondorder 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. (2013Vazan et al. ( , 2015, based on the quotidian EOS (QEOS) of More et al. (1988). Hence ρ i is obtained by solving where ρ 0,i = 0.917 g/cm 3 , P 0,i = 2.07 × 10 11 dyn/cm 2 , and similarly, ρ r is obtained by solving where ρ 0,r = 2.5 g/cm 3 and P 0,r = 1.28×10 12 dyn/cm 2 . 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 sourceQ (energy per unit time per unit volume) in Eq. (8) are as follows: The accretional heating source for the interior iṡ where h(z) is an exponentially decreasing function of the normalized depth z, such that 1 0 h(z)dz = 1. The radioactive energy source iṡ where X 0,j is the initial abundance (relative to rock) of the radionuclide, H j -the energy released per unit mass, and τ j -the characteristic decay time. We consider 40 K, 232 Th, 238 U and 235 U. The heat absorbed/released by melting/refreezing is given bẏ where H 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 ≤ r ≤ R), denoted by V (0 ≤ V ≤ 4π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 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 Although the range of x is fixed, the total volume V (s) changes with time due to accretion of material, V (s) =Ṁ /ρ(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: 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, Gmdm r is calculated and we parametrize the distribution of this energy to obtain the localQ grav for the timestep. The initial and physical parameters assumed for all cases are given in Table 1.

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.

Evolutionary course of a typical case
We choose the case with normalized Safronov parameterθ = 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 ∼ 1 Gyr. Accretional heating at the surface, which surpasses the solar energy by orders of magnitude, causes the surface temperature to rise up to ∼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).
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 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).
(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.
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 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.

Evolution
We now choose a baseline case that has a normalized Safronov parameterθ 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θ between 200 and 500. Similarly, to determine the effect of initial composition we keepθ b and vary the ice mass fraction between 0.3 and 0.7. All models start with a very small embryo (∼0.00015M, where M is final mass) having the same composition as the accreted material.
In Figure 4 we show evolution results forθ 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ρ  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. (fifth panel, Figure 4) and central pressure P c decrease. The mass accretion rateṀ , 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Ṁ , as the accretion heating causes ice to evaporate, the ratio is lowered. Afteṙ 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Ṁ corresponds to the peak in the surface temperature (dashed curves in the last panel).
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θ gradually changes from 200 to 500. As can be seen from Equation 5,Ṁ will change by approximately the same factor. In addition, whenṀ is larger, the available mass is accreted more efficiently, and the accretion heating is more concentrated in time. At the highestθ, 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θ, the effect is barely noticeable.
Only the model corresponding to the highestθ 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 × 10 9 years of accretion and evolution. Profiles are shown forθ 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 7 shows how the internal structures vary for different values ofθ for the initial ice mass fraction X i,b . For the lower values ofθ, 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θ < 400. But forθ = 500, andθ = 400 we clearly see the distinct density discontinuity at near 0.8R ⊕ and 0.4R ⊕ , respectively, marking the core-mantle boundary. This is also reflected in the ice to rock ratio and the temperature profile.

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 (θ = 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θ simultaneously, while Figure  9 does the same for the various parameter profiles after  Fig.4 for different values ofθ. The models assume an ice mass fraction of 0.5.

Gyr of evolution. The figures show that a highṀ
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Ṁ .

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 oḟ M on a smaller mass body, we ran our model for the case of Σ s = 0.1 g cm −2 . This has the effect of reducing the final mass of the planet by half.
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.

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/κ, where κ is the diffusivity, should be comparable to the accretion timescale M/Ṁ . But, for the accretion rate used in Section 4.1, corresponding tô θ = 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 αGMṀ /R, where α is of order unity, we have which yields 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 σT 4 s 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.

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 en- vironment, as well as to constrain the possible interior structure based on measurements of mass and radiusby 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Ṁ (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Ṁ 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 (∼ 10 −6 M ⊕ 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;Voelkel et al. 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,Ṁ 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 ∼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 ⊕ . At 40 au such a body will have a Bondi radius of ∼ 2 × 10 10 cm. For the nebular model we are considering, a sphere of this radius would contain 2.6 × 10 18 g = 4.4 × 10 −10 M ⊕ 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 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. 2020Denman et al. , 2022Wyatt 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 γ = −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 where ρ 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, γ = −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 em- bedded 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. 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.
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,  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. 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. ice=0.2 ice=0.5 ice=0.7 Figure 13. Mass -Radius relation for models with various ice mass fractions andθ = 400. Circles denote the sample of discovered exoplanets, up to 8 M⊕ (based on NASA Exoplanet Archive (2022)). See text for discussion.
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θ, 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 influencê θ. 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.   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.