Radial Drift and Concurrent Ablation of Boulder-Sized Objects

TL;DR

This study investigates how meter-sized bodies in protoplanetary disks drift inward, crossing the water snowline and releasing water and volatiles, which impacts planetary composition and water distribution.

Contribution

It provides a detailed analysis of water transport via radial drift of boulder-sized objects, including a simple sublimation model and thermal evolution coupling.

Findings

Bodies up to 100 meters drift inward past the snowline.

Polluted region extends up to 10% closer to the star than the snowline.

Water transport by meter-sized bodies is significant in disk chemistry.

Abstract

Context. The composition of a protoplanetary disk at a given location does not only depend on temperature and pressure but also on the time dependent transport of matter, such as radial drift of solid bodies, which could release water and other volatile species before disintegration or accretion onto a larger body with potentially considerable implications for the composition of planets. Aims. We perform a parameter study focused on the water depletion of different sized bodies able to cross the water snowline by gas induced radial drift. Methods. Either the analytical Hertz-Knudsen-Langmuir sublimation formula assuming equilibrium temperature within the body or a more involved, numerical model for the internal thermal evolution is coupled with an alpha-disk model. Different properties of the disk and the embedded body are explored. Results. Bodies with radii up to 100 meter drift faster towards the central star than the water snowline, hence, cross it. The region that can be reached before complete disintegration - and is therefore polluted with H2O ice - extends to 10 percent closer to the star than the snowline location. The extent of this polluted region could be multiple times larger in the presence of a dust mantle, which is, however, unlikely to form due to frequent collisions with smaller-than centimeter sized objects. Conclusions. Given a significant abundance of meter sized boulders in protoplanetary disks, the transport of water by radial drift of these bodies towards regions closer to the star than the snowline is not negligible and this flux of volatiles can be estimated for a given distribution of solid body sizes and compositions. A simple expression for surface sublimation is applicable for a homogeneous body consisting of only dust and water ice without a dust mantle.

Tables

Table 1. Table 1: Physical parameters for the nominal disk initial structure and evolution

Parameter	Value
Stellar mass	$1M_{\odot }$
Stellar radius	$2.416R_{\odot }$(a)
Stellar effective temperature	$4377K$(a)
Helium fraction	$0.24$
Power law slope $\beta$	$0.9$(b)
Cut-off radius $R_{\text{out}}$	$50\mathrm{AU}$(b)
Inner boundary radius	$0.03\mathrm{AU}$
Surface density at $5.2\mathrm{AU}$ ${\Sigma }_0$	$268.5g{\mathrm{cm}}^{-2}$(c)
Shakura-Sunyaev $\alpha$-viscosity	$2\text{×}{10}^{-3}$
Photo-evaporation parameter ${\dot{M}}_{\text{wind}}$	${10}^{-7}{M}_{\odot }{\mathrm{yr}}^{-1}$
References.
(a) Baraffe et al. (2015) at $1\mathrm{Myr}$;
(b) Andrews et al. (2010);
(c) MMSN (e.g. Weidenschilling 1977)

Table 2. Table 2: Physical parameters of the cometary nucleus

Parameter	Value
Initial nucleus porosity	$70\%$
Dust mantle porosity	$70\%$
Tortuosity	${\sqrt{2}}^{\left(a\right)}$
Initial dust/ice mass ratio	$1$(b)
Water ice bulk density	$920\mathrm{kg}{m}^{-3}$(c)
Dust bulk density	$3000\mathrm{kg}{m}^{-3}$(c)
Enthalpy of sublimation	$\left[51983.9-20.0904T\right]$ $J{\mathrm{mol}}^{-1}$(d)
Heat conductivity	$\left[0.0028+1.3/T\right]$ $W{m}^{-1}{K}^{-1}$(c)
Volumetric heat capacity	$\left[1582\left(114.8+T\right)\right]$ $J{K}^{-1}{m}^{-3}$(c)
References.
(a) Carman (1956); Mekler et al. (1990); Kossacki & Szutowicz (2006);
(b) Marboeuf et al. (2014);
(c) Marboeuf et al. (2012);
(d) Washburn (1928); Delsemme & Miller (1971)

Equations

$$\sigma T^4=\frac{9}{8}\left(\frac{3}{8}{\tau }_R+\frac{1}{2{\tau }_p}\right)\Sigma \nu {\Omega }_{\text{K}}^2+\sigma T_l^4,$$

$${\tau }_R={\kappa }_R\Sigma .$$

$${\kappa }_R={\kappa }_i{\rho }_0^a{T}^b,$$

$$T_l=T_{\ast }{\left[\frac{2}{3\pi }{\left(\frac{R_{\ast }}{r}\right)}^3+\frac{1}{2}{\left(\frac{R_{\ast }}{r}\right)}^2\frac{H}{r}\left(\frac{d\ln H}{d\ln r}-1\right)\right]}^{1/4},$$

$$\frac{d\Sigma }{dt}=\frac{3}{r}\frac{\partial }{\partial r}\left(r^{1/2}\frac{\partial }{\partial r}\left(r^{1/2}\nu \Sigma \right)\right)+{\dot{\Sigma }}_w.$$

$$\{\begin{array}{ll}{\dot{\Sigma }}_w=0,& \text{if }r\le r_g\\ {\dot{\Sigma }}_w\propto r^{-1},& \text{otherwise},\end{array}$$

$$\eta \equiv 1-v_g/v_K\approx -\frac{r}{2v_K^2{\rho }_0}\frac{dP}{dr}$$

$$\frac{dr}{dt}\approx \{\begin{array}{ll}-\frac{2r\eta }{t_{\text{stop}}},& \text{(Quadratic regime)}\\ -\frac{2r\eta }{t_{\text{stop}}}\left(\frac{s^2}{1+s^2}\right),& \text{(Epstein/Stokes regime)},\end{array}$$

$$s=t_{\text{stop}}{\Omega }_{\text{K}}$$

$$t_{\text{stop}}=\{\begin{array}{ll}\frac{{\rho }_{s}R}{{\rho }_0{v}_{\text{therm}}}& \text{(Epstein)}\\ \frac{2{\rho }_s{R}^2}{3{\rho }_0\lambda v_{\text{therm}}}& \text{(Stokes)}\\ \frac{6{\rho }_{s}R}{{\rho }_0(v_K-v_g)}& \text{(Quadratic)},\end{array}$$

$$r_c=\frac{3}{8}\frac{C_{D\text{,Kn}}{m}_{{\text{H}}_2\text{O}}{\varphi }_{{\text{H}}_2\text{O}}{v}_{{\text{H}}_2\text{O}}}{{\rho }_{\text{grain}}G\frac{M_{\text{nucleus}}}{R_{\text{nucleus}}^2}},$$

$$\varphi(T)=\frac{P^{s}(T)}{\sqrt{2\pi m_{\text{H}_{2}\text{O}}R_{g}T}}\quad\left($\text{\,}\mathrm{mol}\text{\,}{\mathrm{m}}^{-2}\text{\,}{\mathrm{s}}^{-1}$\right)\,,$$

$$\frac{dm}{dt}\bigg{|}_{\text{H}_{2}\text{O}}=\varphi(T)\,m_{\text{H}_{2}\text{O}}\,4\pi R^{2}\,.$$

$$\frac{dR}{dt}=\phi (T)\frac{m_{{\text{H}}_2\text{O}}}{{\rho }_{{\text{H}}_2\text{O}\text{,matrix}}},$$

$$\Sigma(r)=\Sigma_{0}\left(\frac{r}{$5.2\text{\,}\mathrm{AU}$}\right)^{-\beta}\exp{\left[-\left(\frac{r}{R_{\text{out}}}\right)^{(2-\beta)}\right]}\,,$$

$$(T_s^{\text{eq}}{)}^4\approx T_g^4+\frac{C_D{\rho }_0}{32\sigma {\epsilon }_s}|v_g-v_{\text{K}}{|}^3,$$

$$\frac{d{v}_{r,s}}{dt}=\frac{v_{\theta ,s}}{r}-{\Omega }_{\text{K}}^{2}r-\frac{{\Omega }_{\text{K}}}{t_{\text{stop}}}\left(v_{r,s}-v_{r,g}\right),$$

$$\frac{d}{dt}\left(r{v}_{\theta ,s}\right)=-\frac{v_{\text{K}}}{t_{\text{stop}}}\left(v_{\theta ,s}-v_{\theta ,g}\right),$$

$${⟨{\Gamma }_{\text{col}}⟩}_{ti}=h_{ti}^2{a}^2{\Omega }_K{n}_s(m_i)d{m}_i{⟨P_{\text{col}}⟩}_{ti},$$

$$h_{ti}={\left(\frac{m_t+m_i}{3M_{\ast }}\right)}^{1/3}.$$

$$⟨P_{\text{col}}⟩=\min \left({⟨P_{\text{col}}⟩}_{\text{med}},{\left({⟨P_{\text{col}}⟩}_{\text{high}}^{-2}+{⟨P_{\text{col}}⟩}_{\text{low}}^{-2}\right)}^{-1/2}\right),$$

$${⟨P_{\text{col}}⟩}_{\text{high}}=\frac{{\tilde{r}}_p^2}{2\pi }\left(\mathcal{F}(I^{\ast })+\frac{6}{{\tilde{r}}_p}\frac{\mathcal{G}(I^{\ast })}{({\tilde{e}}^{\ast }{)}^2}\right),$$

$$\mathcal{F}(I^{\ast })\equiv 8{\int }_0^1\frac{{I^{\ast }}^{2}E[\sqrt{3(1-{\lambda }^2)}/2]}{[{I^{\ast }}^2+(1-{I^{\ast }}^2){\lambda }^2{]}^2}d\lambda$$

$$\mathcal{G}(I^{\ast })\equiv 8{\int }_0^1\frac{K[\sqrt{3(1-{\lambda }^2)}/2]}{[{I^{\ast }}^2+(1-{I^{\ast }}^2){\lambda }^2]}d\lambda ,$$

$${⟨P_{\text{col}}⟩}_{\text{med}}=\frac{{\tilde{r}}_p^2}{4\pi {\tilde{i}}^{\ast }}\left(17.3+\frac{232}{{\tilde{r}}_p}\right)$$

$${⟨P_{\text{col}}⟩}_{\text{low}}=11.3\sqrt{{\tilde{r}}_p},$$

$${\tilde{e}}^{\ast }\equiv e^{\ast }/h_{ti},{\tilde{i}}^{\ast }\equiv i^{\ast }/h_{ti}$$

$${\tilde{r}}_p\equiv \frac{R_t+R_i}{h_{ti}a}.$$

$${\Gamma }_{\text{col,ti}}=n_V(m_i)\pi (R_t+R_i{)}^2\Delta v\left(1+\frac{v_{\text{esc}}^2}{\Delta v^2}\right),$$

$$v_{\text{esc}}^2=2G\frac{m_t+m_i}{R_t+R_i}.$$

$$f_{\text{solid}}=\frac{0.01}{\sqrt{\frac{{\alpha }_Z}{s+{\alpha }_Z}}},$$

Full text

*Aims. *We perform a parameter study focused on the water depletion of different sized bodies able to cross the water snowline by gas induced radial drift.

*Methods. *Either the analytical Hertz-Knudsen-Langmuir sublimation formula assuming equilibrium temperature within the body or a more involved, numerical model for the internal thermal evolution is coupled with an $\alpha$-disk model. Different properties of the disk and the embedded body are explored.

*Results. *Bodies with radii up to $100\text{\,}\mathrm{m}$ drift faster towards the central star than the water snowline, hence, cross it. The region that can be reached before complete disintegration – and is therefore polluted with H2O ice – extends to $10\text{\,}\mathrm{\char 37\relax}$ closer to the star than the snowline location. The extent of this polluted region could be multiple times larger in the presence of a dust mantle, which is, however, unlikely to form due to frequent collisions with smaller-than centimeter sized objects.

*Conclusions. *Given a significant abundance of meter sized boulders in protoplanetary disks, the transport of water by radial drift of these bodies towards regions closer to the star than the snowline is not negligible and this flux of volatiles can be estimated for a given distribution of solid body sizes and compositions. A simple expression for surface sublimation is applicable for a homogeneous body consisting of only dust and water ice without a dust mantle.

1 Introduction

Recent years of observations and theoretical work on planet formation have stressed the importance of the physics at the various snow- or icelines, i.e. the regions in the protoplanetary disk where a volatile species reaches its condensation temperature. Rather than defining snowlines using the condensation temperature, it is more relevant for planet formation to focus on the presence of water ice in building blocks of planets, which might be different due to dynamical processes. However, we will keep the notion of snowline to refer to the ”classic” snowlines based on temperatures and pressures only.111Here, we generally ignore the fact that a snowline is a surface with a strong dependence on height, but instead only consider the snowline position at the midplane.

Thanks to the continuous improvement of radio-astronomical facilities, such as the Atacama Large Millimeter/Submillimeter Array (ALMA), observations of the carbon monoxide (CO) snowline in certain disks are now possible (Qi et al. 2013, 2015; Schwarz et al. 2016; Nomura et al. 2016; Guidi et al. 2016). The CO snowline is the most accessible snowline to observation because of the low freezing point ($30\text{\,}\mathrm{K}$-$40\text{\,}\mathrm{K}$ (Fray & Schmitt 2009)), implying a large distance to the star, and the high abundance of CO in the disk gas. Unfortunately, the H2O snowline is harder to observe owing to the higher condensation temperature of water. So far, observations were limited to a disk heated by a stellar outburst (Cieza et al. 2016).

The main interest on the water snowline stems from emerging compositional studies of terrestrial planets, motivated by increased precision on measured radii and masses of planets using radial velocity measurements combined with Kepler transit data (e.g. Marcy et al. 2014), or transit timing variation (e.g. Gillon et al. 2017; Grimm et al. 2018, for the TRAPPIST-1 system). Theoretical models of planet formation may help break the degeneracy between planets covered by oceans and those containing H, He atmospheres, while also constraining the mantle composition, if they can put reliable constraints on the volatile content of planets (Adams et al. 2008; Rogers & Seager 2010; Dorn et al. 2015). To achieve this, the compositions of solids and gas in the disk, which are accreted by the (migrating) planets, have to be well constrained over a large region. This will be ultimately necessary to assess the habitability of observed exoplanets.

Finally, the recent studies of Ida & Guillot (2016); Dra̧żkowska & Alibert (2017); Schoonenberg & Ormel (2017) show increased planetesimal formation rates by streaming instability at the H2O snowline, the latter two taking released water vapor into account (see also Ros & Johansen 2013). Those results show the need for proper treatment of all occurring physical processes at the snowline.

A redefinition of the snowline for asteroids in the solar system was explored by Schorghofer (2008) using similar means to what we will present. Schorghofer (2008) found that ice can persist on asteroids of kilometer size up to temperatures of at least $145\text{\,}\mathrm{K}$ over the solar system lifetime and calculated for multiple parameters where these conditions occur. Other recent works aimed at determining the disk composition used chemically evolved (Visser et al. 2009; Eistrup et al. 2016) or equilibrium chemistry (Thiabaud et al. 2015) disks as a basis. However, in these works, no transport of solids, such as radial drift (Weidenschilling 1977) or diffusion processes were modeled. In addition to the relevant chemical evolution, these dynamical effects need to be considered (see Pudritz et al. 2018, for a recent review). Some studies that do address radial transport of different species by pebbles and vapor are Stevenson & Lunine (1988); Dra̧żkowska & Alibert (2017); Schoonenberg & Ormel (2017). Here, we will investigate the potential impact of boulder-sized bodies (sizes from $\text{\,}\mathrm{cm}$ to $100\text{\,}\mathrm{m}$), which is a size regime not treated in the aforementioned works. Such bodies might efficiently transport water ice to regions starwards of the classical water snowline by drifting through the disk faster than the snowline is moving towards the star. In general, this fast drift leads to fast removal of these bodies, which is usually used as an argument not to include those sizes in models. Furthermore, coagulation processes are not efficient enough to let a body grow directly to this size range under nominal conditions (see Blum 2018, for a recent review). However, bodies in this range of sizes are present in the current asteroid population (e.g. Bottke et al. 2005) and models suggest they are naturally produced in collisions of larger bodies (e.g Benz & Asphaug 1999; Bottke et al. 2015).

In this work, we postulate the presence of meter-sized objects and perform a parameter study to determine the region that can be reached by fast drifting bodies crossing the water snowline. This will let us gauge the importance of modeling solid ice transport of fast drifting bodies in the disk. For this, we will identify the dominating processes contributing to thermal heating of the bodies and the parameters and properties influencing the process. Furthermore, we test a simplified, analytic model for the sublimation of water on a boulder-sized body against a more complex, cometary nucleus model (Marboeuf et al. 2012) adjusted to account for the presence of protoplanetary disk gas in the vicinity. The application of a simple analytical expression instead of a more complex numerical model for the sublimation of water ice would help to substantially reduce the computational cost and complexity of future works.

In Sect. 2, we describe the different parts and modes (numerical/analytic) of the model. The results are presented in Sect. 3 and discussed in Sect. 4, where we also describe the validity of our models with regards to all physical processes, which – to our knowledge – influence the results. We conclude in Sect 5.

2 Model description

The model is built on two components: The first component consists of a protoplanetary disk model, including a single solid body embedded in the disk midplane. Its radial drift is calculated and its radially inward motion followed (Sect. 2.1.3). Once the temperature reaches a threshold value (set to 150K), the local disk temperature and radius of the body are used to calculate the body’s thermal and compositional evolution during one timestep of the disk model. Two different modes for this second component are investigated: (a) a numerical model based on the cometary nucleus model by Marboeuf et al. (2012) (Sect. 2.2) or (b) an analytical expression treating the surface ablation (Sect. 2.3). In both cases, mass loss from the body changes its radius, which in turn affects the drift speed.

2.1 Disk

Our $\alpha$-disk model assumes axis-symmetry, vertical hydrostatic equilibrium, flatness ($z\ll r$) and no self-gravity ($M_{\text{Disk}}\ll M_{*}$). The surface density $\Sigma\equiv\int_{-\infty}^{\infty}\rho(z)dz=\rho_{0}H\sqrt{2\pi}$, where $\rho_{0}$ is the midplane density and $H$ is the vertical scale height, is evolved in time and the isothermal sound speed $c_{s}$ is frequently used to abbreviate the ideal gas law as $P=c_{s}^{2}\rho$ (see Sect. 2.1.2 or e.g. Armitage 2019, for more details).

We would like to point out that the purpose of the disk model in this paper is to simulate typical conditions and timescales in the disk midplane only, thus a simple model is sufficient. Of particular importance for this work is the thermal part of the disk model (see Sect. 2.1.1).

2.1.1 Midplane temperature

The midplane temperature $T$ is calculated, as in the model of Hueso & Guillot (2005), by assuming that the disk is geometrically thin, heated viscously and by irradiation from the star, and is optically thick in the radial direction. Instead of solving the radiative transfer numerically, analytic expressions derived by Nakamoto & Nakagawa (1994) are used222Eq. 3.10 in Nakamoto & Nakagawa (1994). In their work, the mid-plane temperature is approximated as a sum of terms accounting for optically thick, optically thin, and stellar contributions:

[TABLE]

where $\sigma$ is the Stephan-Boltzmann constant, $\nu$ is the viscosity described in Sect. 2.1.2 and $\tau_{P}$ is the Planck mean optical depth which is assumed to be $\tau_{P}=2.4\tau_{R}$ as in Nakamoto & Nakagawa (1994). $\tau_{R}$ is the Rosseland mean optical depth, which is defined in terms of the Rosseland mean opacity $\kappa_{R}$ as

[TABLE]

The Rosseland mean opacity is calculated using the modified Alexander/Cox/Stewart opacities by Bell & Lin (1994),

[TABLE]

where the exponents $a$ and $b$ and the factor $\kappa_{i}$ differ in different temperature regimes, i.e. depend on the gas state. The values for $a$, $b$ and $\kappa_{i}$ can be found in the appendix of Bell & Lin (1994).

$T_{l}$ is an effective temperature which includes effects of the irradiation by the star

[TABLE]

where the first term is due to irradiation onto the inner part of a flat disk by a finite sized star (with radius $R_{*}$) and the second term is accounting for irradiation onto the flared outer part. At all radii, we fixed $d\ln H/d\ln r=9/7$ as in Hueso & Guillot (2005). In contrast to their work, however, we did not include a molecular cloud temperature in equation (4) and instead use a fixed floor value of $10\text{\,}\mathrm{K}$ for the temperature.

2.1.2 Disk evolution

For the disk, we assume an $\alpha$-viscosity $\nu=\alpha c_{s}H$, where $\alpha$ is a numerical factor on the order of ${10}^{-4}\text{\,}$ to ${10}^{-2}\text{\,}$ (Shakura & Sunyaev 1973). This viscosity, together with mass and angular momentum conservation, and approximating the orbital velocity of the gas to be Keplerian, lead to a diffusion equation for the disk (Lynden-Bell & Pringle 1974; Pringle 1981)

[TABLE]

We have added an external photo-evaporation term

[TABLE]

where the gravitational radius $r_{g}$ is taken to be equal to $5\text{\,}\mathrm{AU}$ (Veras & Armitage 2004). The mass loss parameter $\dot{M}_{\text{wind}}$, corresponds to the mass that a disk extending to $1000\text{\,}\mathrm{AU}$ would lose due to photo-evaporation. The actual mass loss due to photo-evaporation is approximately $1\text{\,}\mathrm{\char 37\relax}$ of this value because the typical disk only extends to $\sim$100\text{,}\mathrm{AU}. $\dot{M}_{\text{wind}}$ can be chosen to reproduce reasonable lifetimes (Ribas et al. [2014](#bib.bib87)), which is the case for values $\sim${10}^{-7}\text{\,}\mathrm{M_{\odot}}\text{\,}{\mathrm{yr}}^{-1}.

In our disk model, the disk evolution equation (5) is solved numerically on a one dimensional, logarithmically spaced grid in radial direction (Alibert et al. 2005, 2013).

2.1.3 Radial drift

Solid bodies in the disk feel a drag force caused by the different velocities of the gas ($v_{g}$) and the particles ($v_{K}$) which would move with Keplerian speed, in the absence of gas, whereas the former move with a slower velocity due to pressure support (Weidenschilling 1977; Whipple 1972).

To quantify this difference,

[TABLE]

is defined, where the density $\rho_{0}$ and the pressure $P$ are the values taken at the location of the body.

The resulting radial drift depends on the drag regime. To discriminate between the different regimes, the radius $R$ of the solid body is compared to the mean free path $\lambda=m_{\text{H}_{2}}/(\pi d^{2}_{\text{H}_{2}}\rho_{0})$. Here, $d_{\text{H}_{2}}$ and $m_{\text{H}_{2}}$ are the kinetic diameter and mass, respectively, of the hydrogen molecule, i.e. the dominant species in the disk. With that, the drag regime is determined by the following conditions:

• •

If $R<3\lambda/2$, we use the Epstein drag law

• •

Else:

• –

if $R_{e}<27$, where $R_{e}=3(v_{K}-v_{g})R/(\lambda v_{\text{therm}})$ is the microscopic Reynolds number and $v_{\text{therm}}$ the midplane mean thermal velocity, the Stokes drag is used (Rafikov 2004),333This is an approximation introduced by Rafikov (2004), in the literature there is often an additional, intermediate drag regime used between the Quadratic and the Stokes regime (Weidenschilling 1977; Whipple 1972). Additionally, we use here the thermal velocity instead of approximating it as the sound speed. Furthermore, the definition of $R_{e}$ differs by a factor of two compared to Whipple (1972).

• –

if $R_{e}>27$, the bodies drag is governed by the quadratic drag law.

In this work, we use the radial drift formula that is used in Chambers (2008) and is based on the solutions found by Adachi et al. (1976) and similarly by Nakagawa et al. (1986); Takeuchi & Lin (2002):

[TABLE]

where

[TABLE]

is the Stokes number. When switching from Stokes regime to the quadratic regime, the Stokes number should be large, such that there will be no discontinuity of the drift speed. This is the case in our application ($s\sim 5000$ for radii $\sim$100\text{,}\mathrm{m}$$ for which the drag regime changes, see also Fig. 7(b)).

The stopping times differ in the three drag regimes and are given by (Chambers 2008; Takeuchi & Lin 2002; Whipple 1972)

[TABLE]

where $R$ and $\rho_{s}$ are the radius and the density of the solid body.444As in Chambers (2008), the stopping time in the quadratic regime includes a factor of 6, in agreement with Whipple (1972), but in disagreement to the factor 5 in Rafikov (2004). Correspondingly the quadratic drag regime boundary is set to $R_{e}=27$ instead of $R_{e}=20$. It can be easily verified that the drag regimes are chosen such that the stopping time is continuous.

The orbits of the bodies are assumed to be circular, thus the effect of eccentricity and inclination on the drift are omitted. This assumption is reasonable because the eccentricity and inclination get damped by gas drag on shorter timescales than the semi-major axis (Adachi et al. 1976).

We note, that the drift formula (8) does not include the radial gas flow (e.g. Desch et al. (2017)), which does not have a large influence during the fast drifting phases and for large bodies (i.e. high Stokes numbers).

2.2 Cometary nucleus model

We describe here the cometary nucleus model (Marboeuf et al. 2012) that we apply to solid bodies embedded in the protoplanetary disk. The model is able to include multiple volatile species, clathrates and amorphous water ice structures in a 1D (Marboeuf et al. 2012) or 3D model (Marboeuf & Schmitt 2014). Here, we use the 1D model and pure crystalline water ice. Therefore, the equations can be simplified and, for completeness, are presented in that form here, with emphasis on the changes due to the presence of a disk.

The presence of the disk influences the energy sources available to the nucleus, which is discussed in Sect. 2.2.1. Then, the structure and physical model (Sect. 2.2.2) are summarized. Finally, the possible formation of a dust mantle is described in Sect. 2.2.3 and has a large impact on the resulting evolution of snowline crossing bodies.

2.2.1 Energy sources

In the comet related literature, there are usually three main sources of energy considered: solar radiation which heats the surface with energy propagating inward, internal heat release by radioactive isotopes contained in the cometary dust, and the release of latent heat originating from different possible phase changes (crystalline/amorphous ice, clathrates, sublimation) (Klinger 1981).

When considering bodies smaller than $100\text{\,}\mathrm{m}$, radioactive heating is negligible: The simple estimates and values from Merk & Prialnik (2003) 555See their equation 2 with $KT\approx$5\text{,}\mathrm{W}\text{,}{\mathrm{m}}^{-1}, $\rho=$0.5\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}, $t=0$ for heating by 26Al yield a net cooling for radii below $1\text{\,}\mathrm{km}$, due to conduction and radiation from the surface. Because the heat produced by radioactive decay scales $\sim R^{3}$, whereas the cooling scales $\sim R^{2}$, radioactive heating becomes more important for larger objects with radii $>$10\text{,}\mathrm{km}$$.

The crystallization of amorphous water ice is dependent on the temperature and becomes efficient above $100\text{\,}\mathrm{K}$ (Schmitt et al. 1989). It is not clear whether the initial structure of water ice inside comets is crystalline or amorphous. We assume an initial purely crystalline and clathrate-free water ice structure, or crystallization to have happened before the start of the calculation.

For a body embedded in the protoplanetary disk, additional energy sources are available: heat transfer (mainly by isotropic thermal radiation) from the gaseous disk and frictional heating. In contrast, direct irradiation from the sun is suppressed, as the disk is opaque in the midplane (see also Sect. 2.2.2 for implementation).

Frictional heating is caused by the different azimuthal velocities of the disk gas and the solid body. This process is negligible for the main part of this study and discussed in Sect. 4.3.

Hence, the main source of energy for a small body relates to the thermal bath in which the body is moving. This justifies the use of a one dimensional model because of the invariance of heating on the orientation of the body in the disk. We use the local disk gas temperature as a surface temperature for our numerical model. This assumption is discussed in Sect. 4.2.

2.2.2 Structure and physical model

The body is composed of grains, consisting of refractory material, that are enclosed by a mantle of icy water following the model by Greenberg (1988). This structure is drawn schematically in Fig. 1.

For each layer, the diffusion of water vapor through the solid structure of grains is solved using the mass conservation equation. The only processes that can release or bind gas in our crystalline and clathrate-free model are water ice sublimation and condensation.

The flow of gas through the solid matrix can be in different flow regimes; free molecular (Knudsen) flow ($K^{n}>1$) or viscous flow ($K^{n}\ll 1$), depending on the Knudsen number $K^{n}=\lambda/(2r_{p})$, where $\lambda$ is the mean free path of the molecules and $r_{p}$ is the radius of a pore (Knudsen 1909). In addition, we include a transition flow regime for Knudsen numbers $${10}^{-2}\text{,}$<K^{n}<1$ following Fanale & Salvail (1987). One important quantity appearing in the expressions for the flow (Marboeuf et al. 2012), due to the influence of the structure of the pores, is tortuosity. Here, the arc-chord ratio definition of tortuosity is used, i.e. tortuosity is the ratio of the length of a pore to the distance between the endpoints (see Fig. 2).

Energy is conserved at each point inside the nucleus, where heat conduction is modeled using an empirical Hertz factor to account for porosity (Davidsson & Skorov 2002; Prialnik et al. 2004). For detailed equations and explanations refer to Marboeuf et al. (2012).

In Marboeuf et al. (2012) the thermal boundary condition is given by a balance between the solar energy, sublimation of water ice at the surface (if no dust mantle is present), and thermal emission at the outermost layer of the nucleus. In the midplane of a gaseous disk, there is no direct irradiation, since the disk is opaque. Instead, the midplane temperature is used as a boundary condition at the surface of the nucleus and we do not solve the energy balance equation at the surface (this assumption is discussed in Sect. 4.2). The surface sublimation rate follows expression (12).

The gas pressure of the disk in the vicinity of the body is neglected, i.e. the partial pressure of water that is relevant for the sublimation rate is set to zero and the potential influence of the disk gas on diffusion of gaseous species in the interior is not considered. In our model the total amount of gas in the interior is given by the tracked gas flow of the cometary nucleus model and is not including disk gas. We discuss the impact of disk material on the dominating surface sublimation rate in Sect. 4.4.

2.2.3 Dust mantle formation

The solid dust grains that are freed from the rigid structure by sublimation of water ice in the interior can either be ejected from the nucleus or they can accumulate at the surface. The mechanisms for this accumulation are reviewed in Prialnik et al. (2004). To summarize, there are multiple drivers of dust mantle formation that simultaneously appear in a body that undergoes sublimation.

Firstly, gas drag pulls the freed grains outward, but gravity counteracts this process. The magnitude of the gas drag force depends on the grain size, hence there is a critical radius $r_{c}$ of grains that can be ejected. For a slow spinning nucleus the centrifugal force can be neglected, thus (Rickman et al. 1990, equation 7)

[TABLE]

where $C_{D\text{,Kn}}\sim 2$ is the drag coefficient in the free molecular (Knudsen) flow regime (Prialnik et al. 2004), $R_{\text{nucleus}}$ and $M_{\text{nucleus}}$ are the radius and the mass of the whole nucleus, $m_{\text{H}_{2}\text{O}}$ the molar mass, $\phi_{\text{H}_{2}\text{O}}$ the molar flow, and $v_{\text{H}_{2}\text{O}}$ the velocity of water vapor.666This expression differs from the one in Rickman et al. (1990), since, in our case, the gas flow is numerically modeled throughout the nucleus and can be used directly instead of analytically estimating it. Grains with radii larger than $r_{c}$ do not get ejected but instead settle on the nucleus’ surface, already depleted of ice by sublimation. Hence, a porous dust mantle forms on the surface. Huebner et al. (2006) remark that $r_{c}$ only gives an upper size limit, for escaping grains, but smaller grains are not necessarily escaping, as the flow of gas thins above the surface and the grain might fall back onto the nucleus. Furthermore, already in early studies investigating this process, e.g. Brin & Mendis (1979), it was noted that for large dust-to-volatile mass ratio it is impossible to blow off all the freed dust, even though the particles might have radii smaller than $r_{c}$.

The second process comes into play if a dust mantle already exists. The accumulated grains on the surface will interfere with the liberated grains, such that they can no longer pass through the less porous mantle. Hence, they are trapped within the nucleus and further increase the size of the mantle (Shul’man 1972; Rickman et al. 1990).

Furthermore, the dust mantle can break under the gas flow, or its cohesive strength can be large enough to trap not only the dust, but the gas as well (Huebner et al. 2006). In our model, no cohesive forces between the grains are taken into account (Marboeuf et al. 2012). Therefore, we test in section 3.3.3 three cases: (a) the nominal case for which no initial dust mantle is present nor is it allowed to form subsequently, i.e. all the freed dust is lost, (b) an unstable dust mantle case, for which no cohesion forces are taken into account but particles larger than $r_{c}$ are assumed to fall back onto the surface after ejection thereby forming a dust mantle over time, and (c) a constant dust mantle case, with a fixed thickness over the full evolution of the body. In case (c) most of the ejected dust is still lost, but a fraction is kept to keep the artificial constant mantle thickness. These cases differ compared to the work of Schorghofer (2008) who assumed that no dust is lost. This is essentially related to the size of the bodies considered. Small bodies with sizes below hundreds of meters undergoing sublimation (considered in this work) can lose their dust, but larger bodies (considered in Schorghofer (2008)) will be able to keep their dust due to the increased gravity (i.e. $r_{c}\propto R_{\text{nucleus}}^{2}/M_{\text{nucleus}}$ is smaller than all the typical grain sizes).

2.3 Analytical surface ablation model

Instead of invoking the full model from Marboeuf et al. (2012) which is described in Sect. (2.2), an analytic model for the sublimation of water ice from the surface, i.e. ablation, is outlined here, which can be tested against similar models from the literature, e.g. D’Angelo & Podolak (2015), or our full model that includes the very same surface sublimation term.

For a body without a mantle, ablation follows the kinetic theory expression, also known as the Hertz-Knudsen-Langmuir formula, for a free sublimation rate (e.g. Hertz 1882; Delsemme & Miller 1971; Schorghofer 2008; Marboeuf et al. 2012)

[TABLE]

where $P^{s}$ is the water vapor sublimation pressure ($\text{\,}\mathrm{Pa}$), $m_{\text{H}_{2}\text{O}}$ is the molar weight of water and $R_{g}$ is the ideal gas constant ($\text{\,}\mathrm{J}\text{\,}{\mathrm{mol}}^{-1}\text{\,}{\mathrm{K}}^{-1}$). Equation (12) is valid assuming zero partial pressure of water in vicinity of the body. We discuss this approximation in Sect. 4.4. For non-zero pressure with the same temperature, the difference between pressures replaces $P^{s}$ in the equation.

If this amount of water is removed from a layer with thickness $\delta\ll R$ at the surface, the total water mass loss is

[TABLE]

For the refractory part (i.e. dust) of the structure, we assume that the grains are freed in the surface layer and matter gets released immediately adding their contribution to the total mass loss. This can be compared to the case without dust mantle formation of the cometary nucleus model (see Sect. 2.2.3).

Expressed as a decrease in radius, we can write

[TABLE]

where $\rho_{\text{H}_{2}\text{O}\text{,matrix}}$ is the macroscopic water density in the matrix (taking into account porosity). The initial conditions shown in table 2 yield a $\rho_{\text{H}_{2}\text{O}\text{,matrix}}=$276\text{,}\mathrm{kg}\text{,}{\mathrm{m}}^{-3}$$. At fixed porosity, increasing the amount of refractory components reduces the available water, hence reducing $\rho_{\text{H}_{2}\text{O}\text{,matrix}}$ and increasing the total mass loss. The assumption that the dust is freed with the sublimation of the ice is not valid for high dust to water ratios, since the cohesive forces between dust particles would become relevant.

It is noteworthy that if the temperature is kept constant and the body has a homogeneous structure, expression (14) is independent of the body’s total radius, leading to a constant decrease in radius over time. Fig. 3 shows the result of a cometary nucleus model run of an initially $10\text{\,}\mathrm{m}$ sized body, which exhibits this behavior and motivates the analytic sublimation formula. This is true, as long as (i) the total radius is much larger than the radial extent of the surface layer ($R\gg\delta$) and (ii) there are no interior temperature gradients that would change the surface temperature.

Furthermore, we would like to point out that the analytic surface sublimation model is identical to the full cometary nucleus model (Sect. 2.2) assuming: (i) no heat transport of any sort, i.e. the temperature inside the body’s structure is the same as in the disk, (ii) no other species than pure crystalline water ice and dust to be present and (iii) no mantle at the surface to form.

2.4 Initial conditions

2.4.1 Disk

The initial gas surface density of the disk is given by a power law with exponential outer cut-off boundary (as proposed by Andrews et al. 2010) and a normalization constant $\Sigma_{0}$, corresponding to the surface density at approximately $5.2\text{\,}\mathrm{AU}$, which determines the total disk mass.

[TABLE]

where $R_{\text{out}}$ is a constant exponential cut-off radius, $\beta$ is the power law exponent, determining the slope of the surface density profile. The disk evolution is then given by the Shakura-Sunyaev $\alpha$ parameter and the photo-evaporation (Sect. 2.1.2). We leave $\alpha$ fixed and run simulations with (nominal) and without photo-evaporation. The values, which are fixed in all results in this paper, can be seen in table 1. The initial total gas mass in the disk is accordingly $0.05\text{\,}\mathrm{M_{\odot}}$, which is the disk mass that Weidenschilling (1977) uses for the minimum mass solar nebula (MMSN). The star was assumed not to evolve during the disk’s lifetime and the temperature and radius values are taken at a time of $1\text{\,}\mathrm{Myr}$ of stellar evolution according to Baraffe et al. (2015).

In order to gauge the influence of the initial parameters, we varied the total mass, lifetime and heating mode of the disk and the results and changes to the nominal parameters can be found in Sect. 3.3.2.

2.4.2 Solid body

To reduce complexity, we chose to model a body consisting only of water ice and dust, without any other volatile species. Water is the main volatile component (see Marboeuf et al. 2014, about the composition of planetesimals in disks) and the last one to sublimate. Using the parameters listed in table 2, the resulting total density of the body is $\sim$0.42\text{,}\mathrm{g}\text{,}{\mathrm{cm}}^{-3}$$ which is of the same order of magnitude as the recently found bulk density of $0.533\pm 0.006\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$ (Pätzold et al. 2016) and the previous value of $0.470\pm 0.045\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$ by Sierks et al. (2015) of the comet 67P/Churyumov-Gerasimenko. Addition of other volatile species would increase the density to values even closer to these measurements. We chose to represent realistic dust to ice mass ratios ($\sim 1$) (Marboeuf et al. 2014) instead of tuning the ratio to represent measured bulk densities.

The initial location of the body in the disk is set to a distance $10\text{\,}\mathrm{\char 37\relax}$ further away from the star than the snowline, unless otherwise stated. This starting position allows the body to relax to the environment so that initial conditions are forgotten by the time we start computing evaporation.

For heat capacities and conductivities of dust and water ice, we adopted the values listed in Marboeuf et al. (2012) and references therein.

3 Results

We first (Sect. 3.1) present the test cases comparing the two different sublimation models described in Sects. 2.3 and 2.2. Then, we study which bodies are able to cross the snowline (Sect. 3.3). Finally, the results of simulated bodies crossing the snowline that were mainly obtained with the full cometary nucleus model for different varied quantities are presented in Sect. 3.3.

3.1 Comparision between the two sublimation models

Fig. 4 shows the results for a test case, in which we placed a body with an initial radius of $10\text{\,}\mathrm{m}$ and the composition shown in table 2 into the nominal disk (see table 1). The initial semi-major axis is set to $6\text{\,}\mathrm{AU}$ at time zero of the disk evolution. We find almost indistinguishable outcomes between the analytical surface sublimation model and the cometary nucleus model for this particular test.

For the larger, $100\text{\,}\mathrm{m}$ radius case (Fig. 5), the drift timescale is much larger. To save computation time, the body is initially positioned closer to the star than the initial snowline, namely at $4.3\text{\,}\mathrm{AU}$. The initial bulk temperature for the analytical sublimation is, by construction, assumed to be equal to the local disk gas temperature. To estimate the influence of the initial temperature of the body, we ran two different cases with the cometary nucleus model: one with a pre-heated body, i.e. the initial bulk temperature is set to $170\text{\,}\mathrm{K}$, which is the local gas temperature, and one without pre-heating, i.e. an initial bulk temperature of $20\text{\,}\mathrm{K}$. The results are discussed in Sect. 3.4.

For equal initial conditions and under the assumption of initial homogeneous temperature in the nucleus, the analytical solution for the sublimation, i.e. equation (12), and the cometary nucleus model do agree well in the tested size range of bodies (i.e. meters to $100\text{\,}\mathrm{m}$). The agreement worsens with increasing size even in the absence of a dust mantle (see Sect. 3.3.3).

3.2 Snowline versus drift velocity

Temperature and pressure determine the classical snowline position during the evolution of the disk. In our nominal disk (see table 1), the classical water ice line, i.e. the snowline, was determined (see Fig. 6). Due to external photo-evaporation, the disk vanishes almost completely after $2.8\text{\,}\mathrm{Myr}$. As the inner disk surface density decreases, direct irradiation from the central star can invert the cooling of the disk to a heating (via the direct irradiation included in $T_{l}$ in equation 1) in the inner region. Thus, the snowline motion reverts as well. This would not happen in a disk without photo-evaporation, where the disk gradually thins out as a result of the viscous evolution only, i.e. depending solely on $\alpha$, and the snowline motion never changes direction.

In the nominal disk model we calculated the drift speed (see Sect. 2.1.3) of solid bodies in the size range from ${10}^{-2}\text{\,}\mathrm{m}$ to ${10}^{5}\text{\,}\mathrm{m}$ over time, as well as the change of the snowline position. For objects bigger than ${10}^{5}\text{\,}\mathrm{m}$, gas drag is not the relevant source of migration, but the torque exerted by density waves (type I migration) (Goldreich & Tremaine 1979; Ward 1997). The ratio of the body’s drift speed to the snowline speed is shown in Fig. 7(a). Important for our goal is the size range where the transition from bodies moving slower than the snowline to faster than snowline speed lies. In Fig. 7(a) the color code is chosen such that this transition lies in the white region. We found that planetesimals with $R\gtrsim$100\text{,}\mathrm{m}$$ will no longer drift towards the star fast enough to cross the snowline, thus the water ice on these bodies will never sublimate. To help interpret the figure, the drag regime of the different sized bodies is plotted in Fig. 7(b). A size of roughly $100\text{\,}\mathrm{m}$ happens to coincide with the transition from Stokes to quadratic drag regime emphasizing the need to take into account the different drag regimes.

3.3 Parameter study of snowline crossing bodies

In this part of the results section we present the evolution of drifting solid bodies in the protoplanetary disk. These results – obtained using the cometary nucleus model – are presented in Sect. 3.3.1, 3.3.2 and 3.3.3, in which the radius, disk conditions and dust mantle properties are varied.

3.3.1 Initial radius dependence

Figs. 8(a) and 8(b) show the innermost locations reached by different sized bodies, drifting from outer regions of the disk. In the following, we call this position the location of complete disintegration. To get the results, the full cometary nucleus model mode was used. The composition of the different sized bodies was assumed to be equivalent and corresponds to the values given in Sect. 2.4 and table 2. No dust mantle is present in all shown cases, i.e. dust mantle formation is excluded.

When the body reaches high enough temperatures it undergoes ablation and thus loses mass. After shrinking to a radius of $10\text{\,}\mathrm{cm}$ the location is marked as a dot in Fig. 8(a). This location is considered to be the the location of complete disintegration, since a centimeter sized icy body at those temperatures and pressures has a very short lifetime (e.g. Lichtenegger & Kömle 1991). $20\text{\,}\mathrm{b}$odies are modeled starting at different times over the disk lifetime for each evaluated size. Initially, the bodies are located $10\text{\,}\mathrm{\char 37\relax}$ further away from the star than the snowline location at the specific starting time.

The lowest included initial radius is $0.5\text{\,}\mathrm{m}$. Due to numerical and physical assumptions of the model, such as not tracking single grains, lower initial radii are excluded and these pebble sized objects are the main subject of other studies (e.g. Dra̧żkowska & Alibert 2017; Schoonenberg & Ormel 2017).

It can be seen in Fig. 7(a) that bodies with radii on the order of one meter drift the fastest in the protoplanetary disk (see also Weidenschilling 1977; Adachi et al. 1976). Bodies with radii lower than one meter drift slower and thus only cover a small distance after crossing the snowline. The smallest body in our dataset, with an initial radius of $0.5\text{\,}\mathrm{m}$, loses all of its mass and stops very close to the snowline due to its relatively slow drift. The tabulated snowline position values and the sublimation are calculated independently. Therefore, the good agreement of the snowline location in Fig. 8 with the location of complete disintegration of the $0.5\text{\,}\mathrm{m}$ sized body shows that the tabulated snowline position is a reasonable choice of reference.

A larger than meter-sized body with identical composition will also undergo sublimation and thus lose mass. As a consequence, it gradually speeds up until it reaches maximal drift speed at a radius of one meter. This size will be reached closer to the star than the equilibrium position of the snowline. Thus, the object will ultimately drift farther than an initially smaller body until complete disintegration. The difference in the position of complete disintegration will decrease with increasing size, because the initial drift slows down and thus only little distance is traveled before reaching a smaller size. However, the difference will stay bigger than zero, and therefore bigger bodies always cross a larger distance before they completely disintegrate. This asymptotic behavior can be seen in Figs. 8 and 8(b). Since the difference of the position of disintegration between a five meter sized and a ten meter sized body is negligible compared to other effects (see e.g. Sect. 3.3.3) and due to the numerical cost of simulating a large body, no larger sizes were included. There is no reason to expect the position of disintegration of larger bodies to change significantly compared to the one of ten meter sized bodies up to the $100\text{\,}\mathrm{m}$ size boundary, where the bodies can no longer cross the snowline by radial drift (see Sect. 3.2).

To show that the results can be well decoupled from the disk evolution, Fig. 8(b) shows the same results as Fig. 8(a), but instead of measuring the distance from the central star in units of $\text{\,}\mathrm{AU}$, it is measured in units of the snowline position $r_{\text{snowline}}$ ($1\text{\,}$ corresponds to the snowline location, $0\text{\,}$ to the central star). The ratio of the location of complete disintegration to the snowline position stays approximately constant in time.

We would like to point out that this behavior is only found if the bodies are not covered by dust mantles. For bodies with dust mantles, a similar size-dependent behavior is only recovered for exactly equal mantle thicknesses using the constant dust mantle mode described in Sect. 3.3.3. However, the scaling of the mantle thickness depending on size would be the dominant factor but is to our knowledge not well constrained.

The aforementioned time-decoupling of the effect by using units of snowline distance can be used to tentatively explore the overall mass fraction of drifting bodies that should be found at a given distance from the star – measured in units of the snowline position – at all times in the disk (Fig. 9). The mass fraction value shown is an integral over the assumed distribution of bodies (see Sect. 4.1), which was cut such that all included sizes do cross the snowline at all times of the nominal disk evolution (i.e. $1\text{\,}\mathrm{kg}$ to $1\text{\times}{10}^{9}\text{\,}\mathrm{kg}$ corresponding to $8\text{\,}\mathrm{cm}$ to $76\text{\,}\mathrm{m}$). For simplicity, the density was fixed to the nominal value of $0.422\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$ and the analytical sublimation model was used. To help interpret the results, we note that the largest bodies which are abundant for flat slopes do drift the furthest (see Fig. 8), but do not lose a lot of their mass starwards of the snowline. The most efficient transport of mass is achieved by meter sized bodies who are most abundant in the -1.83 slope case, where the location at which $50\text{\,}\mathrm{\char 37\relax}$ of solids remain in the disk is moved from the snowline to two percent starwards of the snowline.

3.3.2 Disk influence

In addition to the nominal disk with values given in table 1, we repeated the calculations for bodies with a radius of $10\text{\,}\mathrm{m}$ embedded in disks for which we modified one parameter compared to the nominal case: a light disk ($M_{\text{disk}}=$0.01\text{,}\mathrm{M_{\odot}}, i.e. $\Sigma_{0}=$53.704\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}), a massive disk ($M_{\text{disk}}=$0.1\text{,}\mathrm{M_{\odot}}, i.e. $\Sigma_{0}=$537.046\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}), a long-lived disk ($\dot{M}_{\text{wind}}=$3\text{\times}{10}^{-9}\text{,}\mathrm{M_{\odot}}\text{,}{\mathrm{yr}}^{-1}$$), and a disk without heating by irradiation of the central star ($T_{l}=0$ in equation 1). As before, the cometary nucleus model is started multiple times in all the different disk evolution calculations. Initially, the body is separated from the star by a distance $10\text{\,}\mathrm{\char 37\relax}$ larger than the classical snowline distance. To cover the full evolution of the disk, the starting times of the individual calculations are scaled with the lifetimes of the different disks. The markers labelled ”no mantle” in Fig. 10 show the temporal mean of all these calculations for the different disk cases with indicated standard deviations ($1\sigma$ error-bars). As in Fig. 8(b), the distance is measured in units of the classical snowline.

We find, that most of the different tested disks have influences on the locations of complete disintegration in units of classical snowline distances on the percent level only. In Fig. 10 it is shown that different disk masses and lifetimes (controlled by photo-evaporation) do not change the result significantly.

However, the non-irradiated disk has a very different temperature profile once the viscous heating is no longer dominating. Thus, the snowline location is altered to a large extent moving, at late times, very close to the star (i.e. to $0.14\text{\,}\mathrm{AU}$ during the calculation of the latest datapoint). Due to the proximity of the snowline to the star, the slope of the surface density and thus the pressure gradient starts to decrease, hence reduces the drift speed. Therefore, the location of complete disintegration is located closer to the snowline, i.e. only $3\text{\,}\mathrm{\char 37\relax}$ below it, which is significantly different from the other cases.

Overall the resulting disintegration locations are robust for the contrasting tested disk cases at many different times. However, the local pressure gradient has a strong influence.

3.3.3 Dust mantle influence

As discussed in Sect. 2.2.3 the formation of a dust mantle on a cometary nucleus is likely. We assume here the same for a disk-embedded body and quantify its potential influence. An important factor is the size of the body, since the process of dust mantle formation depends on the gravitational force. In general, bigger grains are more easily ejected from small nuclei. We consider here thin dust mantles to initially exist on bodies in the gas disk with radii of $10\text{\,}\mathrm{m}$ and evaluate the dust mantle evolution and the influence on the location of complete disintegration. Our model includes dust formation and removal (described in Sect. 2.2.3 and Marboeuf et al. (2012)) without cohesive strength (in the following called the ”unstable” model).

To estimate the extreme case, where the dust mantle cannot be removed, simulating an infinitely large cohesive strength, we artificially set the dust mantle to a constant thickness (called the ”constant” model).

In Fig. 11, it can be seen that using the unstable model the mantle is removed very quickly and the position of complete disintegration differs only slightly from the one without mantle. However, if the dust mantle cannot be removed due to strong cohesive strength and surface sublimation is thus always suppressed, the disintegration location is up to $1\text{\,}\mathrm{AU}$ closer to the star for a $10\text{\,}\mathrm{cm}$ thick mantle. In Figs. 10 and 12, the less extreme case of a $5\text{\,}\mathrm{cm}$ thick constant dust mantle is shown. In the former figure, the temporal mean of the location of complete disintegration for different disks is depicted and in the latter its temporal evolution for the nominal disk is shown.

These results demonstrate the importance of the cohesive strength and thickness of the mantle in determining the thermal evolution of the body. The thickness of the dust mantle is not well constrained by observations, since data is very sparse. The permittivity probe SESAME-PP of the Rosetta mission showed that the first meter is more compact than the rest of the comet 67P (Lethuillier et al. 2016). However, no estimate on the total thickness can be made from this single data point and it is not clear what the composition (possible volatile content) and porosity of this compact layer is. Furthermore, it is not clear how to scale mantle properties from an object with dimensions on the order of kilometer to one with a radius of ten meter.

The large influence of the dust mantle is caused by the change of the sublimation process because free sublimation at the surface is no longer possible if the object is covered by a mantle. Sublimation in the interior still happens, it is however suppressed by the relatively slow diffusion of the released water vapor through the dust mantle, since the small pore radius of the dust mantle limits diffusion.

By analyzing the interior structure of the numerically modeled body, we found that the low thermal conductivity of the dust mantle and porous matrix does not play a dominant role. The body’s interior is heated on short timescales on the order of years for the size range we are interested. This can be seen in Fig. 13, where almost no radial gradient in terms of temperature is visible. This behavior is found for all small body cases with radii of up to $100\text{\,}\mathrm{m}$. For larger bodies or much thicker dust mantles, the picture can change.

The fact that we do not remove the dust mantle by some process is representative of infinite cohesive strength. Thus, the results for a body without dust mantle and the one with constant mantle should be interpreted as lower and upper boundaries for a realistic physical result and the results in Figs. 10 and 12 should be interpreted as such. Measuring the distance from the body to the central star in units of snowline distances again, the location of complete disintegration without dust mantle is at $\sim$0.9\text{,}$$, whereas the one with a dust mantle goes down to $0.5\text{\,}\mathrm{o}$f the snowline distance to the star. However, for a more realistic result, a dust mantle formation and removal model that takes the cohesive strength of the material into account would be needed.

3.4 Internal thermal evolution

To analyze the importance of the internal thermal evolution of the body, we first take a look at the results of the comparision of the analytical surface ablation model and the cometary nucleus model (Sect. 3.1) without a dust mantle.

The underlying assumption of the analytical ablation model is an already equilibrated temperature throughout the body’s full structure and the gas. Hence, as expected, the pre-heated cometary nucleus model results are closer to the analytical model. This good agreement between the two models shows that a numerical treatment of the internal evolution is not necessary for bodies composed mainly of dust and water ice with sizes smaller than $100\text{\,}\mathrm{m}$ and with initially equilibrated temperatures. For many applications of pure water sublimation it is not necessary to invoke a full model keeping track of the internal structure and temperature because heat conduction - and thus temperature equilibration throughout the body - happens at timescales of years. E.g. for a ten meter sized body, the thermal timescale $\tau_{T}\simeq R^{2}\rho c/K$, where $\rho c$ is the density times the heat capacity and $K$ is the heat conductivity (see table 2), is approximately $0.3\text{\,}\mathrm{yr}$. Thus, the internal temperature of a meter sized body spiraling towards the star on timescales of thousands of years is expected to be in thermal equilibrium with the disk.

In the case of a $100\text{\,}\mathrm{m}$ sized body within the snowline but without pre-heating, the body shrinks faster than the heat is transported to the interior. However, the cold interior acts as a heat sink. Thus, heat is conducted to the interior which leads to cooler surface temperatures and slower sublimation (see Fig. 5). This behavior is only reproduced at relatively high temperature regions closer to the central star than the snowline where sublimation is more efficient than heat conduction.

For bodies covered with a dust mantle, the internal temperatures that are reached are significantly higher than for bodies without a dust mantle. Similar is the observed fast heat conduction: for a $10\text{\,}\mathrm{m}$ sized body very little variability in the radial direction is visible (Fig. 13), indicating that sublimation does not lead to a faster shrinking than heat can be conducted to the interior. Thus, the body first becomes isothermal before it disintegrates. No significant thermal insulation increase by the mantle is found: as in the case without a dust mantle, heat conduction acts on timescales of years.

We remark that a numerical treatment of thermal conduction is required to track changes on timescales of years, i.e. on timescales on the order of the orbital period. Hence, assuming an isothermal interior is only valid for objects on almost circular orbits and should not be applied to bodies on eccentric orbits (such as comets).

We conclude that the interior of drifting small bodies ($R<$100\text{,}\mathrm{m}), composed of water ice and dust grains, with zero eccentricity and inclination can be assumed to be isothermal. With that, our analytical model reproduces well the results of the cometary nucleus model. We note that for planetesimals with radii $>$10\text{\,}\mathrm{km}, differentiation due to heating by 26Al (Sect. 2.2.1) can occur (Lichtenberg et al. 2016). For a differentiated body with an ice layer on the surface, sublimation is not hindered by a dust mantle and the analytical sublimation formula becomes appropriate again if no heat is lost to the interior, e.g. if the body is in thermal equilibrium.

4 Discussion

A number of simplifications and assumptions were made to obtain the presented results. These require discussion and some additional calculations that we describe in this section. In addition to that, a successful test of the radial drift formula is presented in appendix A.

4.1 Collisions

In a protoplanetary disk, collisions are a key evolution factor. In this section, we calculate collision rates between our test body – called the target – and a population of other bodies present in the disk – the impactors – and compare them to the timescale of sublimation, which we broadly estimate to be \sim$$1\text{\times}{10}^{3}\text{\,}\mathrm{yr}.

Both, the target and the population of impactors undergo radial drift. Even though radial drift timescales can be as short as $1000\text{\,}\mathrm{o}\mathrm{r}\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{a}\mathrm{l}\,\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{d}\mathrm{s}$ (Armitage 2019), they always remain much larger than the orbital period. We calculate collision rates due to two different processes: (a) caused by coupling to the gas (difference in radial and azimuthal velocities) of different sized bodies and (b) due to eccentricity and inclination distributions induced by gravitational stirring and assuming no radial drift. The latter, which we call the orbital collision rate, is applicable for larger bodies that no longer drift significantly, while the former is applicable for smaller bodies and we call it Stokes collision rate to emphasize the coupling to the gas which is quantified by the Stokes number.

In order to compute the collision rates, a statistical approach using a prescribed distribution function of solids is required. For that, the two body (or ”particle in a box”) approximation (Safronov 1969), i.e. to neglect the influence of the central star, was used to estimate collision rates until Nakazawa et al. (1989a, b); Ida & Nakazawa (1989), and independently Greenzweig & Lissauer (1990, 1992) treated collisions using Hills approximation (Hill 1878). The underlying, adopted probability of a test particle hitting a gravitating object (e.g. a planet) during one orbit was derived by Öpik (1951). However, this approach can only be used for the orbital collision rate. For the gas coupled collision rate, we apply the ”particle in a box” approach (Safronov 1969). A detailed description of the approaches can be found in appendix C and the underlying, assumed mass distribution of bodies as a power law with slope $\alpha$ is described in appendix B.

The resulting integrated orbital collision rates of our nominal target body with impactors larger than the indicated minimum mass (x-axis) are shown in Fig. 14(a). Collisions of the target with large impactors $m_{i}>m_{t}$ are rare for all considered eccentricity and inclination distributions and are thus negligible. Collisions with smaller bodies have to be treated with the Stokes collision rate prescription (appendix C.2) and are shown in Fig. 14(b). In both cases, the rates were integrated from the minimum impactor mass (on the x-axis) to the maximum mass of $1\text{\times}{10}^{24}\text{\,}\mathrm{g}$.777The Stokes collision rates for bodies more massive than the target agree to an order of magnitude precision with the orbital collision rates and do not contribute significantly to the integral.

For a very flat mass distribution, relatively high-energy impacts with bodies with diameters larger than $1\text{\,}\mathrm{m}$ – leading to fragmentation (Windmark et al. 2012; Blum 2018) – are frequent, i.e. are happening about once per $100\text{\,}\mathrm{yr}$, which is comparable to the simulation time of the nominal, drifting body discussed in Sect. 3. If the mass distribution is steeper, the target is less likely to encounter this kind of collisions.

In terms of energetics, the collisions do not contribute large amounts of energy compared to the thermal energy of the body or the total sublimation energy (see table 2): Integrating over all sizes, the kinetic energy is $\lesssim$7\text{\times}{10}^{6}\text{,}\mathrm{J}\text{,}{\mathrm{yr}}^{-1} (using the definition of the ”reduced mass kinetic energy” in Stewart & Leinhardt ([2009](#bib.bib101))), which is $\sim$4\text{\times}{10}^{-4}\text{\,}{\mathrm{yr}}^{-1} of the energy required to heat the body by one kelvin and $\sim$3\text{\times}{10}^{-7}\text{,}{\mathrm{yr}}^{-1}$$ of the total sublimation energy. The yearly collisional energy deposited on the target by a population of impactors in units of the heat capacity of the target is shown in Fig. 15. For a smaller, $1\text{\,}\mathrm{m}$ sized body, the relative numbers increase by almost an order of magnitude. However, the drift and sublimation timescales are also reduced for this smaller body, again resulting in negligible heating by collisions during this stage of the sublimation process.

Most of the energy input results from collisions with bodies with radii smaller than one meter (see Fig. 15). Locally on the targets, the impacts by these smaller bodies are able to erode away target material. This could be the most severe constraint on the applicability of the presented model. The mass encountered per year by the nominal target is $\sim$4\text{\times}{10}^{-4}\text{,} times its own mass. Windmark et al. ([2012](#bib.bib113)) fitted erosion efficiencies based on laboratory experiments for silicate grains. Using their velocity and mass dependent fit (Windmark et al. [2012](#bib.bib113), equation 17), the total eroded mass relative to the target mass is $\sim$8\text{\times}{10}^{-2}\text{\,}\mathrm{\char 37\relax}\text{\,}{\mathrm{yr}}^{-1}. This would imply that the assumption of a collision free sublimation is only applicable on timescales $\lesssim$10\text{,}\mathrm{yr}. Furthermore, for collisions involving impactors with sizes comparable to the target, fragmentation of both objects can happen and only a remnant with mass smaller than the masses of each object remains but Fig. [14(b)](#S4.F14.sf2) shows that this comparable-size case is rare and can be safely ignored. The erosion rates in the regime of collisions with meter-sized bodies is not well studied and applying the fit of Windmark et al. ([2012](#bib.bib113)) is therefore an extrapolation with its inherent flaws. Using the lower limit of the erosional prescription for porous icy agglomerates used in Krijt et al. ([2015](#bib.bib56)) yields smaller erosion rates $\sim$2\text{\times}{10}^{-2}\text{\,}\mathrm{\char 37\relax}\text{\,}{\mathrm{yr}}^{-1} translating to erosion of less than $7.5\text{\,}\mathrm{\char 37\relax}$ of the bodies mass over the time where sublimation was active ($T>$150\text{,}\mathrm{K}$$) in the numerical simulations shown in e.g. Fig. 4.

We conclude, that during the crucial short phase ($\sim$100\text{,}\mathrm{yr}$-$1000\text{,}\mathrm{yr}$$), where sublimation and fast radial drift take place, collisions with small bodies are happening frequently. The results presented in Sect. 3 are only strictly valid if either the surface density of solids is reduced (e.g. by not converting all solids to pebbles, by less efficient settling, or by accumulation of solids in planets), or erosion is less efficient in the relevant mass regime (large uncertainties of extrapolation of laboratory experiments). Otherwise, erosion by collisions could become an additional relevant mass loss mechanism. In terms of thermal energy, collisions do not heat the body, thereby justifying the thermal balance model we presented in Sect. 2. The fast erosion of meter sized bodies is the main argument against their presence in disks. In this work, however, we postulate their presence, which could be justified by frequent enough fragmentation of larger bodies.

We note, that the retention of a dust mantle is very hard to achieve if collisions are eroding away the uppermost layers of the body. A mantle of centimeter thickness is eroded by collisions with pebbles in less than $10\text{\,}\mathrm{yr}$.

4.2 Gas versus surface temperatures

D’Angelo & Podolak (2015) showed that for small bodies ($R<$10\text{,}\mathrm{km}$$) the bulk temperature of the body is in equilibrium with the gas after less than 500 orbits (see Fig. 20 in D’Angelo & Podolak (2015)). In their work, the entire body was heated and reached equilibrium temperature in this amount of time, whereas in our full model, we only assume equilibration of the temperature in an uppermost, thin layer. Instantaneous heat exchange from the thermal bath, i.e. the disk, to the body is thus well justified.

4.3 Frictional heating

In this work, heating due to interactions with the non-Keplerian gas is not taken into account. D’Angelo & Podolak (2015) calculate the equilibrium value of the surface temperature for their planetesimals to be (D’Angelo & Podolak 2015, equation 38)

[TABLE]

where $C_{D}$ is the drag coefficient, $\sigma$ is the Stefan-Boltzmann constant and $\varepsilon_{s}$ is the thermal emissivity (for a black body $\varepsilon_{s}=1$). To derive this equilibrium value, a fraction of $C_{D}/4$ of the total collisional energy is assumed to be transmitted as heat to the body, which corresponds to an upper limit (Podolak et al. 1988).

For a simple estimate, using the typical values $\rho_{0}=${10}^{-9}\text{,}\mathrm{g}\text{,}{\mathrm{cm}}^{-3}, $T_{g}=$140\text{\,}\mathrm{K} and $\eta=$4\text{\times}{10}^{-3}\text{,}$$, frictional heating yields a negligibly small temperature increase of $8\text{\times}{10}^{-4}\text{\,}\mathrm{K}$ for a black body. Only in very dense regions of the disk or potentially in the atmospheres of planets could a significant change occur.

4.4 Water vapor pressure

Disk water vapor can change the sublimation rate and even lead to deposition of water onto bodies if present in high enough abundance ($P_{\text{vapor}}>P^{s}(T)$). In an ideal case, where all other solid bodies in the disk do not move radially and the disk is not evolving, $P_{\text{vapor}}<P^{s}(T)$ everywhere. Hence, no water would condense onto the surface of a body drifting by. However, if fast drifting pebbles are present, the water vapor surface density can be replenished by diffusion of the freshly released vapor starwards of the snowline (Ros & Johansen 2013; Schoonenberg & Ormel 2017; Dra̧żkowska & Alibert 2017). Another source for out of thermal equilibrium water vapor could potentially be stellar outbursts which episodically heat up the disk (Hartmann & Kenyon 1996). To study constraints for deposition or suppressed sublimation in detail, a model including the evolution of all solids and water vapor in the disk would be needed. If a significant amount of vapor is transported further away from the star than the snowline, it could be deposited onto a drifting body, reducing (for Stokes numbers $s>1$) the drift speed, which allows for even more deposition, potentially leading to growth to planetesimal size.

A local source for enhanced water vapor pressure could also be the drifting body itself, due to exhibiting a coma-like region with enhanced partial pressure of water, reducing the sublimation rate. The transport of gas or vapor away from the body is not modeled here and would differ from the case of a comet due to the disk gas interacting with the released vapor and the lack of solar wind, which is absorbed in the disk. Our assumption of no increased local partial pressure due to the coma is consistent with a complete erosion of the coma by the disk gas.

If the partial pressure of vapor is smaller than the sublimation pressure $P^{s}$ but not zero, it reduces the sublimation rate compared to the nominal results without vapor. In Fig. 16 we show the influence of an artificially chosen, exponential increase of partial pressure of water vapor (motivated by results of the pebble based evaporation models of Schoonenberg & Ormel (2017)) up to one percent of the total pressure. The location to reach the percent level is set to where the temperature is $176.6\text{\,}\mathrm{K}$ to avoid deposition of water.

Under these assumptions, the water vapor moves the location of complete disintegration farther in towards the star. Hence, in the context of presence of water ice in solid bodies (the ”dynamical” snowline), the results for the case without a mantle and without water vapor in the disk is an upper boundary for the dynamical snowline location.

5 Conclusions

We presented the application of a cometary nucleus model to disk-embedded, radially drifting, spherical bodies, tracking the thermodynamic evolution of the object. The body is assumed to consist of only dust and water ice. Different properties of the disks and the drifting bodies were explored and the time evolution of the disk was taken into account. The main focus of the work was to constrain the regions that can be reached by drifting icy bodies, ultimately determining the zone where some water can be incorporated in solids and thus be accreted by growing terrestrial planets.

Here, we summarize the key findings:

Almost independently of the properties and temporal evolution of the disk, drifting bodies with radii $\geq$1\text{,}\mathrm{m}$$ can transport water ice at least ten percent closer to the star than the location of the ”classical” snowline before they completely disintegrate. 2. 2.

If surface sublimation is not impeded in any way, e.g. by the presence of a dust mantle, it is the dominant process for the evolution of the object and can be modeled in a simple, analytic way with good agreement with the results of a full numerical model. 3. 3.

These results are applicable to bodies with radii ranging from meters to $100\text{\,}\mathrm{m}$. Smaller bodies never experience fast radial drift, therefore the effect is suppressed, whereas bodies larger than $100\text{\,}\mathrm{m}$ do not drift fast enough to even cross the snowline. In the range from tens to hundreds of meters, the difference in locations of complete disintegration is small. This implies that if bodies in this size range are present, a quantifiable smearing of the water snowline results. In the absence of meter-sized bodies, the snowline is given by the local disk properties only. 4. 4.

A dust mantle covering the body suppresses surface sublimation and forces the internally released vapor to diffuse through the mantle. For the extreme case of a non-breakable mantle, this results in icy bodies drifting starwards to about one half of the classical snowline position. However, the presence and formation of a global dust mantle on a body embedded in a protoplanetary disk is hindered by collisions with pebble sized objects, because these collisions occur at relative velocities typically leading to net mass loss, i.e. erosion of the uppermost layers of the body. In particular, for bodies smaller than meter-size a dust mantle is highly unlikely to be kept due to the – relative to the total mass – large erosion rates.

Multiple processes were not included and several assumptions were made to obtain the above results. We identified two key processes that could affect our results and which should be addressed in future works:

• •

Collisions with different sized bodies, mainly stemming from the difference in radial and azimuthal velocities due to gas drag, are frequent for large solid fractions and would mostly lead to erosion. For a model including multiple bodies of different sizes, tracking the thermodynamic evolution of each is necessary to properly estimate the general outcome. A potential approach to reduce the numerical cost is to use the analytic surface sublimation expression for objects with low thermal variability over the course of an orbit (i.e. low eccentricity and inclination).

• •

Specific water vapor pressures influence sublimation rates and thus the results are sensitive to this. To get fully consistent results, it is required to take the water vapor distribution in the disk, including the contribution of the evaporation bodies, into account.

The approximations of imposing the gas temperature as surface temperature, neglecting frictional heating, and using a simple formula for the radial drift is found justified for all discussed parameters.

Of particular interest for future works is to test and potentially apply the analytic sublimation formula in complete N-body terrestrial planet formation models (as suggested by Coleman & Nelson 2016). This would also include larger than $100\text{\,}\mathrm{m}$ sized bodies because they could be moved across the snowline by N-body interactions (e.g. scattering or resonant trapping) and bodies on significantly eccentric and inclined orbits for which further research on their thermal evolution is necessary. Furthermore, we did not include different chemical species that could either be present as icy layers on the grains or as clathrates and we leave the treatment of the evolution of bodies at different, potentially observable ice lines to future works. Moreover, the influence of including amorphous water ice and the phase change to crystalline ice along with a model for dust mantle growth including cohesive strength and predicting the properties (pore size, porosity, tortuosity, thickness) of the formed mantle should also be addressed in the future for a complete model.

The presented and proposed steps will help to constrain compositions and available masses for terrestrial planet growth, which will be increasingly required to match the precisions on future observational constraints on planetary compositions.

Acknowledgements.

The authors thank the anonymous referee who’s comments helped improve the manuscript and Cléa Serpollier for her early investigative work in the subject. This work has been carried out within the frame of the National Centre for Competence in Research PlanetS funded by the Swiss National Science Foundation (SNSF). The authors acknowledge the financial support from the SNSF under grant 200020_172746.

Appendix A Radial drift formula

Treating the fastest drifting bodies in protoplanetary disks correctly, might require additional changes to the radial drift formula shown in equation (8). We show here the validity of the assumptions made to derive this form of the equation, i.e. assuming orbit averaged drift $\tau_{\text{drift}}\gg\tau_{\text{orb}}$, neglecting terms quadratic in $\eta$, assuming no radial acceleration ($dv_{r,s}/dt=0$), and setting the particle’s azimuthal speed to Keplerian ($v_{\theta,s}=v_{\text{K}}$) in the derivative term (first term in equation (18))888See Takeuchi & Lin (2002) for an instructive derivation of the simplified equations. The equations of motion in the disk plane ($v_{z}=0$) are (Takeuchi & Lin 2002)

[TABLE]

where the subscript $s$ and $g$ are for the solid body and gas, respectively, and $t_{\text{stop}}$ is given by equation (10).

For a test, we assumed $v_{r,g}=0$ and solved the equations of motion numerically. The results can be seen in Fig. 17 and are compared to the results of the analytical equation (8) with the same initial conditions. After one stopping time has passed (dashed vertical line), the initially Keplerian azimuthal ($v_{\theta,s}(t=0)=v_{\text{K}}$) speed slowed down to an equilibrium value and the analytical expression (8) reproduces the differential equation results well. The radial drift speed is slowing down because the body moves towards the star ($dr/dt\propto r\eta(r)$). The order of percent difference after equilibration of the azimuthal speed can thus be explained by this non-zero $dv_{r,s}/dt$, which is assumed to be zero to derive equation (8).This difference is small compared to the uncertainties of the other processes treated in this work.

Appendix B Mass distribution

Before assessing the collision rates, we discuss here briefly the required mass or size distributions of the bodies in the disk. The differential mass distribution $n(m_{i})$ is defined such that $n(m_{i})dm_{i}$ is the number of bodies with masses in the interval [$m_{i}$,$m_{i}+dm_{i}$]. We describe $n(m_{i})$ as a power-law with exponent $\alpha$. Data constraining the mass distribution is mainly available from solar system observations or from theoretical works treating collisional cascades or related effects (e.g. Dohnanyi 1969; Tanaka et al. 1996; Makino et al. 1998; Benz & Asphaug 1999; Jutzi et al. 2010; Pan & Schlichting 2012; Belton 2015). The observational data is either gathered by direct measurements of Jupiter family comets (e.g Fernández et al. 1999; Tancredi et al. 2006; Fernández et al. 2013), trans-Neptunian objects (e.g. Bernstein et al. 2004) or asteroids (e.g. Gladman et al. 2009) or inferred from distributions of craters on planets, satellites or other minor planets (Zahnle et al. 2003; Singer et al. 2019). The measured and predicted values of the slope $\alpha$ of the differential mass distribution – assuming a fixed density for converting size distributions – lie in the interval $[-1.5,-2.1]$. We note that multiple studies found different slopes for bodies with radii smaller than $\text{\,}\mathrm{km}$ (Zahnle et al. 2003; Fernández & Morbidelli 2006; Fernández et al. 2013; Singer et al. 2019). Nevertheless, we adopt for our order of magnitude estimates simple unbroken power laws with three fixed values for $\alpha$: the upper (-1.5 as Morbidelli & Rickman 2015) and lower (-2.1 slightly lower than the -2.05 found by Belton 2015) limits and the $\alpha$ resulting from the self-similar solution to the collisional cascade (-1.83, Dohnanyi 1969).

To avoid divergence, the distribution needs to be cut at a lower and an upper boundary. We choose an upper limit to the mass of $1\text{\times}{10}^{24}\text{\,}\mathrm{g}$, corresponding to a radius of $827\text{\,}\mathrm{km}$. The lower cut is of particular importance for the resulting collisions rates. We choose the typical pebble that can form by coagulation for the lower limit: according to laboratory experiment it has a size of \sim$$1\text{\,}\mathrm{cm} and a corresponding mass of \sim$$1\text{\,}\mathrm{g} (Blum 2018). To not underestimate the amount of solids, we assume a $100\text{\,}\mathrm{\char 37\relax}$ conversion of dust to pebbles and larger bodies.

Appendix C Collisions

C.1 Orbital collision rate

The averaged number of collisions between a target with mass $m_{t}$ and a population of bodies with mass $m_{i}$ per unit time is written as (Nakazawa et al. 1989a; Ohtsuki 1999; Inaba et al. 2001)

[TABLE]

where $\braket{P_{\text{col}}}_{ti}$ is a non-dimensional mean collision rate between bodies with masses $m_{t}$ and $m_{i}$ that is independent of the total number of bodies, but depends on the common semi-major axis, the radii and masses of the two bodies, and the mass of the central star. The brackets indicate, that the mean collision rate is an average over all eccentricities and inclinations given by a Reyleigh-type distribution function with eccentricity (inclination) dispersions $e^{*}$ ($i^{*}$), which also influence the mean collision rate (Inaba et al. 2001). $n_{s}(m_{i})dm_{i}=\Sigma_{s}/m_{i}\,n(m_{i})dm_{i}$ is the surface number density of bodies with masses between $m_{i}$ and $m_{i}+dm_{i}$ with $\Sigma_{s}$ the surface density of solids, whereas $h_{ti}$ is the reduced Hill radius of two bodies with masses $m_{t}$ and $m_{i}$ given by

[TABLE]

For the entire range of realistic eccentricity and inclination distributions, Inaba et al. (2001) found that numerical results are well reproduced if the non-dimensional mean collision rate is set to

[TABLE]

where the individual parts are

• •

[TABLE]

where $I^{*}\equiv\tilde{i}^{*}/\tilde{e}^{*}$,

[TABLE]

where $K$ and $E$ are the complete elliptic integrals of the first and second kinds,

• •

[TABLE]

• •

[TABLE]

with the reduced eccentricity and inclination dispersions

[TABLE]

and

[TABLE]

C.2 Stokes collision rate

For drifting, small particles we use the ”particle in a box” approximation, where the collision rate $\Gamma_{\text{col,ti}}$ of a gravitating target with radius $R_{t}$ and with a number of impactors with radius $R_{i}$ is (Safronov 1969)

[TABLE]

with the volume number density of impactors $n_{V}(m_{i})$, the relative velocity $\Delta v$ of the impactors with respect to the target and the mutual escape speed

[TABLE]

To estimate the collision rate, we use the squared sum of the difference in radial (equation (8)) and azimuthal drift velocities of the target and impactors according to their size. The azimuthal difference in velocity is given by $\eta v_{k}|\frac{1}{1+s_{i}^{2}}-\frac{1}{1+s_{t}^{2}}|$, where $s_{i}$ and $s_{t}$ are the Stokes numbers of the impactor and the target (Birnstiel et al. 2016). This excludes the additional velocity components due to Brownian motion and turbulence. For a more complete discussion of relative velocities we refer to Ormel & Cuzzi (2007).

The number density of a given size of impactors can be estimated given three ingredients: their mass (or size) distribution discussed above, the density of the gas and the local dust to gas ratio $f_{\text{solid}}$, which is locally enhanced due to dust settling and can for low masses be described by (Youdin & Lithwick 2007; Birnstiel et al. 2016)

[TABLE]

where $s$ is the Stokes number and $\alpha_{Z}$ is a dimensionless parameter for turbulent diffusion in the vertical direction. Here, we assumed a global dust to gas fraction of $0.01\text{\,}$ and we assume $\alpha_{Z}\simeq\alpha$, which is true on orders of magnitude level (Youdin & Lithwick 2007). For larger than meter-sized bodies the settling is no longer well described by the processes considered in Youdin & Lithwick (2007), since gravitational interactions between the particles become important. Therefore $f_{\text{solid}}$ is restricted to be $f_{\text{solid}}\leq 1$ and this maximum is reached at $R_{2}\approx$5\text{,}\mathrm{m}$$. At meter size, the inclination caused by the viscous stirring of a larger planetesimal in the vicinity of our target leads to approximately the same elevation above the midplane as the reduced scale height (Ida 1990; Thommes et al. 2003; Fortier et al. 2013).