© ESO 2019

1 Introduction

Recent years of observations and theoretical work on planet formation have highlighted the importance of the physics at the various snow- or icelines, that is, 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¹.

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 snowline most accessible to observation because of the low freezing point (30–40 K Fray & Schmitt 2009), implying a large distance to the star, and the high abundance of CO in the disk gas. Unfortunately, the H₂O 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 and 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), and Schoonenberg & Ormel (2017) show increased planetesimal formation rates by streaming instability at the H₂O 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 physical processes occurring 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 present here. Schorghofer (2008) found that ice can persist on asteroids of kilometer size up to temperatures of at least 145 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), and Schoonenberg & Ormel (2017). Here, we investigate the potential impact of boulder-sized bodies (sizes from cm to 100 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 toward the star. In general, this fast drift leads to fast removal of these bodies, which is usually used as an argument to not 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 and analytic) of the model. The results are presented inSect. 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 150 K), 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 α-disk model assumes axis-symmetry, vertical hydrostatic equilibrium, flatness (z ≪ r), and no self-gravity (M_Disk ≪ M_*). The surface density $\Sigma \equiv {\displaystyle {\int }_{-\infty }^{\infty }\rho }(z)\text{d}z={\rho }_{0}H\sqrt{2\pi }$, where ρ₀ 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 used². In their work, the mid-plane temperature is approximated as a sum of terms accounting for optically thick, optically thin, and stellar contributions: $$\sigma T^4=\frac{9}{8}\left(\frac{3}{8}{\tau }_{\text{R}}+\frac{1}{2{\tau }_{\text{p}}}\right)\Sigma \nu {\Omega }_{\text{K}}^2+\sigma T_l^4,$$(1)

where σ is the Stephan–Boltzmann constant, ν is the viscosity described in Sect. 2.1.2, and τ_P is the Planck mean optical depth, which is assumed to be τ_P = 2.4τ_R as in Nakamoto & Nakagawa (1994). τ_R is the Rosseland mean optical depth, which is defined in terms of the Rosseland mean opacity κ_R as $${\tau }_{\text{R}}={\kappa }_{\text{R}}\Sigma .$$(2)

The Rosseland mean opacity is calculated using the modified Alexander/Cox/Stewart opacities by Bell & Lin (1994), $${\kappa }_{\text{R}}={\kappa }_i{\rho }_0^a{T}^b,$$(3)

where the exponents a and b, and the factor κ_i differ in different temperature regimes, that is, depend on the gas state. The values for a, b, and κ_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 $$T_l=T_{*}{\left[\frac{2}{3\pi }{\left(\frac{R_{*}}{r}\right)}^3+\frac{1}{2}{\left(\frac{R_{*}}{r}\right)}^2\frac{H}{r}\left(\frac{d\ln H}{d\ln r}-1\right)\right]}^{1/4},$$(4)

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∕dlnr = 9∕7 as in Hueso & Guillot (2005). In contrast to their work, however, we did not include a molecular cloud temperature in Eq. (4) and instead use a fixed floor value of 10 K for the temperature.

2.1.2 Disk evolution

For the disk, we assume an α-viscosity ν = αc_sH, where α is a numerical factor on the order of 10⁻⁴ –10⁻² (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) $$\frac{\text{d}\Sigma }{\text{d}t}=\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.$$(5)

We have added an external photo-evaporation term $$\{\begin{array}{ll}{\dot{\Sigma }}_w=0,\hfill & \text{if }r\le r_g\hfill \\ {\dot{\Sigma }}_w\propto r^{-1},\hfill & \text{otherwise},\hfill \end{array}$$(6)

where the gravitational radius r_g is taken to be equal to 5 AU (Veras & Armitage 2004). The mass loss parameter Ṁ_wind, corresponds to the mass that a disk extending to 1000 AU would lose due to photo-evaporation. The actual mass loss due to photo-evaporation is approximately 1% of this value because the typical disk only extends to ~100 AU. Ṁ_wind can be chosen to reproduce reasonable lifetimes (Ribas et al. 2014), which is the case for values ~ 10⁻⁷ M_⊙ yr⁻¹. In our disk model, the disk evolution Eq. (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, $$\eta \equiv 1-v_g/v_{\text{K}}\approx -\frac{r}{2v_{\text{K}}^2{\rho }_0}\frac{\text{d}P}{\text{d}r}$$(7)

is defined, where the density ρ₀ 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_{{\text{H}}_2}^2{\rho }_0)$. Here, $d_{{\text{H}}_2}$ and $m_{{\text{H}}_2}$ are the kinetic diameter and mass, respectively, of the hydrogen molecule, which is the dominant species in the disk. With that, the drag regime is determined by the following conditions:

• If R < 3λ∕2, we use the Epstein drag law

• Else:

  • if R_e < 27, where R_e = 3(v_K − v_g)R∕(λv_therm) is the microscopic Reynolds number and v_therm the midplane mean thermal velocity, the Stokes drag is used (Rafikov 2004)³,

  • 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) and Takeuchi & Lin (2002): $$\frac{\text{d}r}{\text{d}t}\approx \{\begin{array}{cc}-\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}$$(8)

where $$s=t_{\text{stop}}{\Omega }_{\text{K}}$$(9)

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 ~ 5000 for radii ~ 100 m for which the drag regime changes).

The stopping times differ in the three drag regimes and are given by (Chambers 2008; Takeuchi & Lin 2002; Whipple 1972) $$t_{\text{stop}}=\{\begin{array}{cc}\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_{\text{K}}-v_g)}& \text{(Quadratic)},\end{array}$$(10)

where R and ρ_s are the radius and the density of the solid body⁴. 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). This 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 describehere 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 considerable 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 m, radioactive heating is negligible: The simple estimates and values from Merk & Prialnik (2003)⁵ for heating by ²⁶Al yield a net cooling for radii below 1 km, due to conduction and radiation from the surface. Because the heat produced by radioactive decay scales ~ R³, whereas the cooling scales ~R², radioactive heating becomes more important for larger objects with radii > 10 km.

The crystallization of amorphous water ice depends on the temperature and becomes efficient above 100 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ⁿ > 1) or viscous flow (Kⁿ ≪ 1), depending on the Knudsen number Kⁿ = λ∕(2r_p), where λ 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⁻² < Kⁿ < 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, where tortuosity isthe 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 Eq. (12).

The gas pressure of the disk in the vicinity of the body is neglected, meaning 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 does not include disk gas. We discuss the impact of disk material on the dominating surface sublimation rate in Sect. 4.4.

Fig. 1 Schematic view of structure model. Adapted from Marboeuf & Schmitt (2014).

Fig. 2 Tortuosity of a path through a porous structure. In (A) a path through the material is shown, in (B) the length of the pore L and the distance between the endpoints X is indicated. Tortuosity is defined as L∕X. Image adapted from O’Connell et al. (2010) under a creative common licence.

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, Eq. (7)) $$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}},$$(11)

where C_D,Kn ~ 2 is the drag coefficient in the free molecular (Knudsen) flow regime (Prialnik et al. 2004), R_nucleus and M_nucleus are the radius and the mass of the whole nucleus, $m_{{\text{H}}_2\text{O}}$ the molar mass, ${\varphi }_{{\text{H}}_2\text{O}}$ the molar flow, and $v_{{\text{H}}_2\text{O}}$ the velocity of water vapor⁶. 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, for example in 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 also the gas (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 Sect. 3.3.3 three cases: (a) the nominal case for which no initial dust mantle is present nor is it allowed to form subsequently, that is, 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 fixedthickness over the full evolution of the body. In case (c) most of the ejected dust is still lost, but a fraction iskept 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) are 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, called ablation, is outlined here. This can be tested against similar models from the literature, for instance 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) $$\phi (T)=\frac{P^s(T)}{\sqrt{2\pi m_{{\text{H}}_2\text{O}}{R}_{\text{g}}T}}\left(\text{mole} {\text{m}}^{-2} {\text{s}}^{-1}\right),$$(12)

where P^s is the water vapor sublimation pressure (Pa), $m_{{\text{H}}_2\text{O}}$ is the molar weight of water and R_g is the ideal gas constant (J mole⁻¹ K⁻¹). 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 δ ≪ R at the surface, the total water mass loss is $$\frac{\text{d}m}{\text{d}t}{|}_{{\text{H}}_2\text{O}}=\phi (T)m_{{\text{H}}_2\text{O}}4\pi R^2.$$(13)

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 $$\frac{\text{d}R}{\text{d}t}=\phi (T)\frac{m_{{\text{H}}_2\text{O}}}{{\rho }_{{\text{H}}_2\text{O,matrix}}},$$(14)

where ${\rho }_{{\text{H}}_2\text{O,matrix}}$ is the macroscopic water density in the matrix (taking porosity into account). The initial conditions shown in Table 2 yield a ${\rho }_{{\text{H}}_2\text{O,matrix}}=276 \text{kg} {\text{m}}^{-3}$. At fixed porosity, increasing the amount of refractory components reduces the available water, hence reducing ${\rho }_{{\text{H}}_2\text{O,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, Eq. (14) is independentof the body’s total radius, leading to a constant decrease in radius over time. Figure 3 shows the result of a cometary nucleus model run of an initially ten-meter-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 ≫ δ) 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, meaning 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.

Fig. 3 Almost linear decrease in radius over time (green line, left axis) for a fixed surface temperature of 169 K using the cometary nucleus model. The derivative d R∕d t is plotted inorange (right axis).

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 Σ₀, corresponding to the surface density at approximately 5.2 AU, which determines the total disk mass. $$\Sigma (r)={\Sigma }_0{\left(\frac{r}{5.2 \text{AU}}\right)}^{-\beta }\exp \left[-{\left(\frac{r}{R_{\text{out}}}\right)}^{(2-\beta )}\right],$$(15)

where R_out is a constant exponential cut-off radius, β is the power law exponent, determining the slope of the surface density profile. The disk evolution is then given by the Shakura–Sunyaev α parameter and the photo-evaporation (Sect. 2.1.2). We leave α 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 M_⊙, 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 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 ~ 0.42 g cm⁻³ which is of the same order of magnitude as the recently found bulk density of (0.533 ± 0.006) g cm⁻³ (Pätzold et al. 2016) and the previous value of (0.470 ± 0.045) g cm⁻³ 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 (~ 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% 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

Figure 4 shows the results for a test case, in which we placed a body with an initial radius of 10 m and the composition shown in Table 2 into the nominal disk (see Table 1). The initial semi-major axis is set to 6 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 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 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 K, which is the local gas temperature), and one without pre-heating (i.e., an initial bulk temperature of 20 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 (see Eq. (12)) and the cometary nucleus model do agree well in the tested size range of bodies (i.e., meters to 100 m). The agreement worsens with increasing size even in the absence of a dust mantle (see Sect. 3.3.3).

Fig. 4 Comparison of the comet model solving the internal structure and the analytical solution (Eq. (12)) for ten-meter-sized bodies. The solid lines show the distance to the star (left axis) with dots representing the locations where the bodies shrank to a size of 10 cm while the dash-dotted lines show the remaining mass fraction (right axis). The lines of the two different model solutions are essentially indistinguishable. The initial position is 10% above the snowline location at time zero in the nominal disk. The barely visible kink in the mass fractions at 600 yr is due to reaching the threshold temperature of 150 K, where the sublimation models are started.

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, also called the snowline, was determined (see Fig. 6). Due to external photo-evaporation, the disk vanishes almost completely after 2.8 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 Eq. (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, that is, depending solely on α, and the snowlinemotion 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⁻² m to 10⁵ m over time, as well as the change of the snowline position. For objects bigger than 10⁵ 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 the top panel of Fig. 7. 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 the top panel of Fig. 7 the color code is chosen such that this transition lies in the white region. We found that planetesimals with R ≳ 100 m will no longer drift toward 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 the bottom panel of Fig. 7. A size of roughly 100 m happens to coincide with the transition from Stokes to quadratic drag regime emphasizing the need to take into account the different drag regimes.

Fig. 5 As Fig. 4, but for a hundred-meter-sized body. The initial position is chosen starwards of the snowline, which is at 4.3 AU. One of the cometary nucleus model runs is started with a low initial bulk temperature of 20 K (green line), whereas the other is pre-heated to the local gas temperature at the starting location (170.16 K), which is theimplicit assumption of the Analytical Sublimation. The lines of the analytical sublimation model and the preheated comet model are barely distinguishable. The local gas temperatures at the end of the calculation are 174.95, 175.07, 174.94 K for the preheated, analytical and the cold model respectively.

Fig. 6 Surface density evolution for the nominal disk. The dashed, blue line shows the snowline position.

Fig. 7 Drift velocity compared to water snow line velocity (top panel) and drag regime (bottom panel) in an irradiated disk with photo-evaporation. All the bodies with sizes in the red area in the top panel cross the snowline, since they drift faster than it moves toward the central star. The snowline is determined using tabulated values for the temperature and pressure. After approximately 2.4 Myr, the snowline position starts to move away from the star, due to the disks dispersal. Thus, the ratio of the body’s drift speed to the snowline speed is negative and the log in the top panel is no longer defined and the ratio is set to a value of 12 to indicate that all the bodies will cross the snowline in that phase. In order to smooth out numerical artifacts, we applied a Gaussian filter in horizontal direction.

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 Sects. 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

The top andbottom panel of Fig. 8 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, meaning 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 cm the location is marked as a dot in the top panel of Fig. 8. This location is considered to be 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). Twenty bodies are modeled starting at different times over the disk lifetime for each evaluated size. Initially, the bodies are located 10% further away from the star than the snowline location at the specific starting time.

The lowest included initial radius is 0.5 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 the top panel of Fig. 7 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 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 body with a size of 0.5 m 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 thanan 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 Fig. 8. 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 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, the bottom panel of Fig. 8 shows the same results as the top panel, but instead of measuring the distance from the central star in units of AU, it is measured in units of the snowline position r_snowline (1 corresponds to the snowline location, 0 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 kg to 1 × 10⁹ kg corresponding to 8 to 76 m). For simplicity, the density was fixed to the nominal value of 0.422 g cm⁻³ 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% of solids remain in the disk is moved from the snowline to two percent starwards of the snowline.

Fig. 8 Comparison of locations of complete disintegration of different sized bodies without a dust mantle. In the top panel, the distance to the star is measured in AU and the dots represent the locations where the body shrank to a size of 10 cm. The dashed, cyan line indicates the evolving position of the P–T tabulated snowline, and the regions where only icy solid bodies, only water depleted solid bodies, and the region that is injected with drifting icy bodies are colored and labeled. In the bottom panel, the same data is shown but measured in units of the evolving, tabulated snowline position (1 corresponds to the snowline position, 0 to the central star).

Fig. 9 Remaining mass fraction in the overall population of bodies crossing the snowline (1 kg ≤ m ≤ 1 × 109 kg) with shaded bands indicating the standard deviation due to the evolving disk. The mass shown is an integral over a distribution of masses with the indicated power-law slope and a mean over time in the disk. More details can be found at the end of Sect. 3.3.1.

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 m embedded in disks for which we modified one parameter compared to the nominal case: a light disk (M_disk = 0.01 M_⊙, i.e., Σ₀ = 53.704 g cm⁻²), a massive disk (M_disk = 0.1 M_⊙, i.e., Σ₀ = 537.046 g cm⁻²), a long-lived disk (Ṁ_wind = 3 × 10⁻⁹M_⊙ yr⁻¹), and a disk without heating by irradiation of the central star (T_l = 0 in Eq. (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% 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 labeled “no mantle” in Fig. 10 show the temporal mean of all these calculations for the different disk cases with indicated standard deviations (1σ error-bars). As in the bottom panel of Fig. 8, 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 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, that is, it is only 3% 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.

Fig. 10 Mean relative locations of complete disintegration in different disks for initially ten-meter-sized bodies with and without constant dust mantles. As in the bottom panel of Fig. 8, the distance from the central star is measured in units of snowline positions (snowline location at 1, star at 0). The mean over the calculations at different times is taken and the 1σ error is indicated. Refer to the text for the different disk properties.

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 considerhere thin dust mantles to initially exist on bodies in the gas disk with radii of 10 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 AU closer to the star for a 10 cm thick mantle. In Figs. 10 and 12, the less extreme case of a 5 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 permittivityprobe 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 mantleproperties 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 nolonger 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 temperatureis visible. This behavior is found for all small body cases with radii of up to 100 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 ~ 0.9, whereas the one with a dust mantle goes down to 0.5 of 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.

Fig. 11 Sublimation comparison of a ten-meter-sized bodies with different dust mantle thicknesses and removal processes. The legend is ordered in increasing sublimation time. The smooth line marks the location of the body in time (left axis), while the dots at the end of the line indicate shrinking to a radius of 10 cm as in Figs. 8 and 12. The dash-dotted lines indicate the mass fraction compared to the initial mass of the same colored case (right axis). The kink that is visible in the mass fraction of the unstable (initially 5 cm thick) mantle stems from the mantle breaking up at that point in time.

Fig. 12 Locations of complete disintegration in the nominal disk of initially ten-meter-sized bodies with and without a dust mantle. The classical water ice line (snowline) in the disk is indicated for reference.

Fig. 13 Interior temperature of an initially ten-meter-sized body covered by a 10 cm thick dust mantle. The number of layers is reduced to 15 compared to nominal runs for better visibility and 60 timesteps are merged into one block. The uppermost, dark framed layer shows the dust mantle. Outside the body, the local disk temperature is plotted. A radial temperaturegradient is only barely visible.

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 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. For instance, for a ten-meter-sized body, the thermal timescale τ_T ≃ R²ρc∕K, where ρc is the density times the heat capacity and K is the heat conductivity (see Table 2), is approximately 0.3 yr. Thus, the internal temperature of a meter-sized body spiraling toward the star on timescales of thousands of years is expected to be in thermal equilibrium with the disk.

In the case of a hundred-meter-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 ten-meter-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, that is, 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 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 km, differentiation due to heating by ²⁶Al (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, for instance, 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 ~1 × 10³ yr.

Both, the target and the population of impactors undergo radial drift. Even though radial drift timescales can be as short as 1000 orbitalperiods (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), that is, to neglect the influence of the central star, was used to estimate collision rates until Nakazawa et al. (1989a,b), and 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 α 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 the top panel of Fig. 14. 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 the bottom panel of Fig. 14. In both cases, the rates were integrated from the minimum impactor mass (on the x-axis) to the maximum mass of 1 × 10²⁴ g⁷.

For a very flat mass distribution, relatively high-energy impacts with bodies with diameters larger than 1 m – leading to fragmentation (Windmark et al. 2012; Blum 2018) – are frequent, that is, are happening about once per 100 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 ≲ 7 × 10⁶ J yr⁻¹ (using the definition of the “reduced mass kinetic energy” in Stewart & Leinhardt 2009), which is ~ 4 × 10⁴ yr⁻¹ of the energy required to heat the body by one kelvin and ~ 3 × 10⁷ yr⁻¹ 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, meter-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 ~ 4 × 10⁴ times its own mass. Windmark et al. (2012) fit erosion efficiencies based on laboratory experiments for silicate grains. Using their velocity andmass dependent fit (Windmark et al. 2012, Eq. (17)), the total eroded mass relative to the target mass is ~ 8 × 10²% yr⁻¹. This would imply that the assumption of a collision free sublimation is only applicable on timescales ≲10 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 the bottom panel of Fig. 14 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) 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) yields smaller erosion rates ~ 2 × 10²% yr⁻¹ translating to erosion of less than 7.5% of the bodies mass over the time where sublimation was active (T > 150 K) in the numerical simulations shown in Fig. 4.

We conclude, that during the crucial short phase (~ 100–1000 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 yr.

Fig. 14 Collision rates of the nominal target body (Table 2) with a radius of 10 m integrated over impactor masses larger than the indicated minimum mass. Top panel: orbital collision rate and bottom panel: stokes collision rate. Results are shown for three different slopes of the impactor mass distribution and in the bottom panel additionally for three different eccentricity and inclination values. $e_{\text{eq}}^{M-m}$ is the equilibrium eccentricity in the vicinity of a mars mass perturber (Ida 1990; Thommes et al. 2003). The red shaded region depicts collisions more frequent than once every thousand years and the target mass mt is indicated.

Fig. 15 Collisional energy calculated with the Stokes collision rate, integrated over impactor masses larger than the indicated minimum mass. The energy is measured in units of the energy required to heat the body by one Kelvin. The target properties and impactor mass distributions are the same as in Fig. 14 and the target mass mt is indicated.

4.2 Gas versus surface temperatures

D’Angelo & Podolak (2015) showed that for small bodies (R < 10 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, that is, 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, Eq. (38)) $${\left(T_s^{\text{eq}}\right)}^4\approx T_g^4+\frac{C_{\text{D}}{\rho }_0}{32\sigma {\epsilon }_s}|v_g-v_{\text{K}}{|}^3,$$(16)

where C_D is the drag coefficient, σ is the Stefan–Boltzmann constant, and ε_s is the thermal emissivity (for a black body ε_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 ρ₀ = 10⁻⁹g cm⁻³, T_g = 140 K and η = 4 × 10³, frictional heating yields a negligibly small temperature increase of 8 × 10⁻⁴ 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_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_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 K to avoid deposition of water.

Under these assumptions, the water vapor moves the location of complete disintegration farther in toward 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:

• 1.

  Almost independently of the properties and temporal evolution of the disk, drifting bodies with radii ≥1 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.

  If surface sublimation is not impeded in any way, for instance 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.

  These results are applicable to bodies with radii ranging from meters to hundred meters. Smaller bodies never experience fast radial drift, therefore the effect is suppressed, whereas bodies larger than 100 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.

  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, that is, 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 hundred-meter-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.

Fig. 16 Evolution of a ten-meter-sized body in the nominal disk calculated with the analytical surface ablation model, with and without water vapor pressure. The water vapor increases exponentially depending on the local disk temperature up to a maximum of one percent of the local pressure.

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 thefastest drifting bodies in protoplanetary disks correctly, might require additional changes to the radial drift formula shown in Eq. (8). We show here the validity of the assumptions made to derive this form of the equation, that is, assuming orbit averaged drift τ_drift ≫ τ_orb, neglecting terms quadratic in η, assuming no radial acceleration (d v_r,s∕dt = 0), and setting the particle’s azimuthal speed to Keplerian (v_θ,s = v_K) in the derivative term (first term in Eq. (A.2))⁸. The equations of motion in the disk plane (v_z = 0) are (Takeuchi & Lin 2002) $$\begin{array}{l}\frac{\text{d}{v}_{r,s}}{\text{d}t}=\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{\text{d}}{\text{d}t}\left(r{v}_{\theta ,s}\right)=-\frac{v_{\text{K}}}{t_{\text{stop}}}\left(v_{\theta ,s}-v_{\theta ,g}\right),\end{array}$$

where the subscript s and g are for the solid body and gas, respectively, and t_stop is given by Eq. (10). For a test, we assumed v_r,g = 0 and solved the equations of motion numerically. The results can be seen in Fig. A.1 and are compared to the results of the analytical Eq. (8) with the same initial conditions. After one stopping time has passed (dashed vertical line), the initially Keplerian azimuthal (v_θ,s(t = 0) = v_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 toward the star (d r∕dt ∝ rη(r)). The order of percent difference after equilibration of the azimuthal speed can thus be explained by this non-zero d v_r,s∕dt, which is assumed to be zero to derive Eq. (8). This difference is small compared to the uncertainties of the other processes treated in this work.

Fig. A.1 Radial drift speed comparision between the approximate Eq. (8) and the numerical solution to the differential equations. The radius is 1 m, the initial location is at 2 AU, temperature and midplane density are constant over the disk and set to 140 K and 2.5 × 1010 g cm−3. Initially thebody moves with Keplerian speed in azimuthal direction and no radial velocity. The dashed, vertical line marks the stopping time.

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 α. 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, 2013; Tancredi et al. 2006), 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 α 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 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 α: 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 α 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 × 10²⁴ g, corresponding to a radius of 827 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 ~1 cm and a corresponding mass of ~1 g (Blum 2018). To not underestimate the amount of solids, we assume a 100% 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) $${\langle {\Gamma }_{\text{col}}\rangle }_{ti}=h_{ti}^2{a}^2{\Omega }_{\text{K}}{n}_s(m_i)\text{d}{m}_i\langle P_{\text{col}}{}_{ti},$$(C.1)

where ${\langle P_{\text{col}}\rangle }_{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 = Σ_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 Σ_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 $$h_{ti}={\left(\frac{m_t+m_i}{3M_{*}}\right)}^{1/3}.$$(C.2)

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 $$\langle P_{\text{col}}\rangle =\min \left({\langle P_{\text{col}}\rangle }_{\text{med}},{\left({\langle P_{\text{col}}\rangle }_{\text{high}}^{-2}+{\langle P_{\text{col}}\rangle }_{\text{low}}^{-2}\right)}^{-1/2}\right),$$(C.3)

where the individual parts are

• $${\langle P_{\text{col}}\rangle }_{\text{high}}=\frac{{\tilde{r}}_p^2}{2\pi }\left(\mathcal{F}(I^{*})+\frac{6}{{\tilde{r}}_p}\frac{\mathcal{G}(I^{*})}{{({\tilde{e}}^{*})}^2}\right),$$(C.4)

  where I^*≡ĩ^*∕ẽ^*, $$\mathcal{F}(I^{*})\equiv 8{\displaystyle {\int }_0^1\frac{I^{*}{}^{2}E[\text{}\sqrt{3(1-{\lambda }^2)}/2]}{{[I^{*}{}^2+(1-I^{*}{}^2){\lambda }^2]}^2}}\text{d}\lambda$$(C.5)

  and $$\mathcal{G}(I^{*})\equiv 8{\displaystyle {\int }_0^1\frac{K[\text{}\sqrt{3(1-{\lambda }^2)}/2]}{[I^{*}{}^2+(1-I^{*}{}^2){\lambda }^2]}}\text{d}\lambda ,$$(C.6)

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

• $${\langle P_{\text{col}}\rangle }_{\text{med}}=\frac{{\tilde{r}}_p^2}{4\pi {\tilde{i}}^{*}}\left(17.3+\frac{232}{{\tilde{r}}_p}\right)$$(C.7)

• $${\langle P_{\text{col}}\rangle }_{\text{low}}=11.3\sqrt{{\tilde{r}}_p},$$(C.8)

with the reduced eccentricity and inclination dispersions $${\tilde{e}}^{*}\equiv e^{*}/h_{ti},{\tilde{i}}^{*}\equiv i^{*}/h_{ti}$$(C.9)

and $${\tilde{r}}_p\equiv \frac{R_t+R_i}{h_{ti}a}.$$(C.10)

C.2 Stokes collision rate

For drifting, small particles we use the particle in a box approximation, where the collision rate Γ_col,ti of a gravitating target with radius R_t and with a number of impactors with radius R_i is (Safronov 1969) $${\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),$$(C.11)

with the volume number density of impactors n_V(m_i), the relative velocity Δv of the impactors with respect to the target and the mutual escape speed $$v_{\text{esc}}^2=2G\frac{m_t+m_i}{R_t+R_i}.$$(C.12)

To estimate the collision rate, we use the squared sum of the difference in radial (Eq. (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_solid, which is locally enhanced due to dust settling and can for low masses be described by (Youdin & Lithwick 2007; Birnstiel et al. 2016) $$f_{\text{solid}}=\frac{0.01}{\sqrt{\frac{{\alpha }_Z}{s+{\alpha }_Z}}},$$(C.13)

where s is the Stokes number and α_Z is a dimensionless parameter for turbulent diffusion in the vertical direction. Here, we assumed a global dust to gas fraction of 0.01 and we assume α_Z ≃ α, 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_solid is restricted to be f_solid ≤ 1 and this maximum is reached at R₂ ≈ 5 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).

References

All Tables

All Figures

	Fig. 1 Schematic view of structure model. Adapted from Marboeuf & Schmitt (2014).
In the text

	Fig. 2 Tortuosity of a path through a porous structure. In (A) a path through the material is shown, in (B) the length of the pore L and the distance between the endpoints X is indicated. Tortuosity is defined as L∕X. Image adapted from O’Connell et al. (2010) under a creative common licence.
In the text

	Fig. 3 Almost linear decrease in radius over time (green line, left axis) for a fixed surface temperature of 169 K using the cometary nucleus model. The derivative d R∕d t is plotted inorange (right axis).
In the text

Fig. 4 Comparison of the comet model solving the internal structure and the analytical solution (Eq. (12)) for ten-meter-sized bodies. The solid lines show the distance to the star (left axis) with dots representing the locations where the bodies shrank to a size of 10 cm while the dash-dotted lines show the remaining mass fraction (right axis). The lines of the two different model solutions are essentially indistinguishable. The initial position is 10% above the snowline location at time zero in the nominal disk. The barely visible kink in the mass fractions at 600 yr is due to reaching the threshold temperature of 150 K, where the sublimation models are started.

In the text

Fig. 5 As Fig. 4, but for a hundred-meter-sized body. The initial position is chosen starwards of the snowline, which is at 4.3 AU. One of the cometary nucleus model runs is started with a low initial bulk temperature of 20 K (green line), whereas the other is pre-heated to the local gas temperature at the starting location (170.16 K), which is theimplicit assumption of the Analytical Sublimation. The lines of the analytical sublimation model and the preheated comet model are barely distinguishable. The local gas temperatures at the end of the calculation are 174.95, 175.07, 174.94 K for the preheated, analytical and the cold model respectively.

In the text

	Fig. 6 Surface density evolution for the nominal disk. The dashed, blue line shows the snowline position.
In the text

Fig. 7 Drift velocity compared to water snow line velocity (top panel) and drag regime (bottom panel) in an irradiated disk with photo-evaporation. All the bodies with sizes in the red area in the top panel cross the snowline, since they drift faster than it moves toward the central star. The snowline is determined using tabulated values for the temperature and pressure. After approximately 2.4 Myr, the snowline position starts to move away from the star, due to the disks dispersal. Thus, the ratio of the body’s drift speed to the snowline speed is negative and the log in the top panel is no longer defined and the ratio is set to a value of 12 to indicate that all the bodies will cross the snowline in that phase. In order to smooth out numerical artifacts, we applied a Gaussian filter in horizontal direction.

In the text

Fig. 8 Comparison of locations of complete disintegration of different sized bodies without a dust mantle. In the top panel, the distance to the star is measured in AU and the dots represent the locations where the body shrank to a size of 10 cm. The dashed, cyan line indicates the evolving position of the P–T tabulated snowline, and the regions where only icy solid bodies, only water depleted solid bodies, and the region that is injected with drifting icy bodies are colored and labeled. In the bottom panel, the same data is shown but measured in units of the evolving, tabulated snowline position (1 corresponds to the snowline position, 0 to the central star).

In the text

	Fig. 9 Remaining mass fraction in the overall population of bodies crossing the snowline (1 kg ≤ m ≤ 1 × 109 kg) with shaded bands indicating the standard deviation due to the evolving disk. The mass shown is an integral over a distribution of masses with the indicated power-law slope and a mean over time in the disk. More details can be found at the end of Sect. 3.3.1.
In the text

Fig. 10 Mean relative locations of complete disintegration in different disks for initially ten-meter-sized bodies with and without constant dust mantles. As in the bottom panel of Fig. 8, the distance from the central star is measured in units of snowline positions (snowline location at 1, star at 0). The mean over the calculations at different times is taken and the 1σ error is indicated. Refer to the text for the different disk properties.

In the text

Fig. 11 Sublimation comparison of a ten-meter-sized bodies with different dust mantle thicknesses and removal processes. The legend is ordered in increasing sublimation time. The smooth line marks the location of the body in time (left axis), while the dots at the end of the line indicate shrinking to a radius of 10 cm as in Figs. 8 and 12. The dash-dotted lines indicate the mass fraction compared to the initial mass of the same colored case (right axis). The kink that is visible in the mass fraction of the unstable (initially 5 cm thick) mantle stems from the mantle breaking up at that point in time.

In the text

	Fig. 12 Locations of complete disintegration in the nominal disk of initially ten-meter-sized bodies with and without a dust mantle. The classical water ice line (snowline) in the disk is indicated for reference.
In the text

	Fig. 13 Interior temperature of an initially ten-meter-sized body covered by a 10 cm thick dust mantle. The number of layers is reduced to 15 compared to nominal runs for better visibility and 60 timesteps are merged into one block. The uppermost, dark framed layer shows the dust mantle. Outside the body, the local disk temperature is plotted. A radial temperaturegradient is only barely visible.
In the text

Fig. 14 Collision rates of the nominal target body (Table 2) with a radius of 10 m integrated over impactor masses larger than the indicated minimum mass. Top panel: orbital collision rate and bottom panel: stokes collision rate. Results are shown for three different slopes of the impactor mass distribution and in the bottom panel additionally for three different eccentricity and inclination values. $e_{\text{eq}}^{M-m}$ is the equilibrium eccentricity in the vicinity of a mars mass perturber (Ida 1990; Thommes et al. 2003). The red shaded region depicts collisions more frequent than once every thousand years and the target mass mt is indicated.

In the text

	Fig. 15 Collisional energy calculated with the Stokes collision rate, integrated over impactor masses larger than the indicated minimum mass. The energy is measured in units of the energy required to heat the body by one Kelvin. The target properties and impactor mass distributions are the same as in Fig. 14 and the target mass mt is indicated.
In the text

	Fig. 16 Evolution of a ten-meter-sized body in the nominal disk calculated with the analytical surface ablation model, with and without water vapor pressure. The water vapor increases exponentially depending on the local disk temperature up to a maximum of one percent of the local pressure.
In the text

Fig. A.1 Radial drift speed comparision between the approximate Eq. (8) and the numerical solution to the differential equations. The radius is 1 m, the initial location is at 2 AU, temperature and midplane density are constant over the disk and set to 140 K and 2.5 × 1010 g cm−3. Initially thebody moves with Keplerian speed in azimuthal direction and no radial velocity. The dashed, vertical line marks the stopping time.

In the text