© ESO 2021

1 Introduction

Concentric rings and gaps in the millimeter continuum emission associated with pebbles are commonly observed features in protoplanetary disks (e.g., ALMA Partnership 2015; Andrews et al. 2018). Recent evidence suggests that at least some of these radial pebble concentrations are due to trapping at pressure bumps (Dullemond et al. 2018). The origin of pressure bumps is often ascribed to planet–disk interactions, where a growing protoplanet carves a gap in the gas disk, leading to the formation of a pressure maximum at the planetary gap edge (e.g., Pinilla et al. 2012; Dipierro et al. 2015; Fedele et al. 2018; Zhang et al. 2018; Favre et al. 2019). The direct observation of protoplanets orbiting within protoplanetary disk gaps support this scenario (Pinte et al. 2019, 2020). Since the ratio of solids-to-gas in these pressure bumps can rise to values significantly higher than the global one, they are favorable spots for planetesimal formation via the streaming instability (SI; Youdin & Goodman 2005; Johansen et al. 2009; Lyra et al. 2009; Bai & Stone 2010; Carrera et al. 2015; Yang et al. 2017).

In Eriksson et al. (2020, hereafter Paper I), we performed global 1D simulations of dust evolution and planetesimal formation via the SI in a protoplanetary disk with multiple gap-opening planets. We found that planetesimals form at the gap edges for a wide range of planetary masses, particle sizes, and disk parameters. A similar study by Stammler et al. (2019) shows that planetesimal formation via the streaming instability can further explain the observed range of optical depths among dust rings in the DSHARP survey. Results from the first 3D simulations of planetesimal formation in a pressure bump were reported in Carrera et al. (2021), who concluded that planetesimal formation in pressure bumps is a robust process in the case of cm-sized particles. They did not, however, find any planetesimal formation in the case of mm-sized particles, but a search based on higher resolution and a broader parameter range is needed to assess this finding.

The above studies assert that planetesimal formation in pressure bumps is likely to be a common process. In this work, we assume that these pressure bumps are formed by growing planets and we consider what the fate of planetesimals formed at these planetary gap edges might be. To answer this question, we performed a suite of N-body simulations, taking into account the effect of gas drag on the planetesimals and mass loss via ablation. We consider two planetary systems: theSolar System with Jupiter and Saturn and an HL Tau inspired system with three planets. The main points we want to address are: (1) the extent to which gravitational interactions with the forming planets eventually redistribute the planetesimals formed at their gap edge; and (2) whether the frictional heating of planetesimals on eccentric orbits drives the production of a significant amount of dust through surface ablation; and (3) how common collisions are between pairs of planetesimals and between planetesimals and planets.

We find that planetesimals that form at the edges of planetary gaps do not remain at their birth location. Gravitational scattering by the embedded planets causes the planetesimals to spread out across the disk, resulting in high orbital eccentricities. If the planetesimals end up in the inner part of the disk, efficient ablation due to high surface temperatures results in fast mass loss; this prevents any planetesimal from entering the innermost disk region. Collisions between pairs of planetesimals are rare and should not affect the outcome of our simulations, whereas planet–planetesimal collisions do occur, and often more so in the HL Tau system than in the Solar System. In general, the closer proximity to the Sun results in much more ablation in the Solar System than in the HL Tau system. If this ablated material re-condenses to form pebbles in the disk midplane, the result is a significant flux of pebbles both interior and exterior of Jupiter’s orbit. Since the ablation efficiency depends strongly on the gas density, this pebble flux should only be present in relatively young and massive disk.

In Sect. 2, we present our disk model, our prescriptions for gas drag and mass ablation, as well as our calculations of condensation temperatures and planetsimal–planetesimal collision timescales. The numerical setup and initial conditions for our simulations are presented in Sect. 3. Results from the Solar System simulations are presented in Sect. 4 and results from the HL Tau simulations are presented in Sect. 5. In Sect. 6, we calculate the pebble flux that derives from planetesimal ablation and discuss how this fits into our understanding of the isotope dichotomy in the Solar System. The effect of planet-planetesimal collisions are discussed in Sect. 7. In Sect. 8, we summarize the key findings of the paper. Further information about the ablationmodel and some additional figures can be found in Appendix A–F.

2 Theory

We modeled the protoplanetary disk gas using a static 1D accretion disk. We added an acceleration to the equation of motion of the planetesimals in order to account for the friction force they are subject to as they move through the gas. This friction acts to heat up the planetesimals and causes a vaporization of the surface layers. We calculate the rate at which material is ablated from the planetesimals and update their mass at each timestep in the simulation. The temperature and mass loss rate are highly dependent on the composition of the planetesimals, which is set by their formation location relative to the location of the major icelines in the disk. We consider the disk to contain the following volatiles – H₂O, CO₂ and CO – and calculate their respective condensation temperatures. Finally, we estimate the frequency of planetesimal–collisions using an analytic expression.

2.1 Disk model

2.1.1 Surface density Σ

The gas surface density profile was defined to be that of a 1D steady disk, $$\Sigma =\frac{{\dot{M}}_0}{3\pi \nu }\exp \left[-\frac{r}{r_{\text{out}}}\right],$$(1)

where Σ is the gas surface density, Ṁ₀ is the disk accretion rate, ν is the kinematic viscosity of the disk, r is the semimajor axis, and r_out is the position of the outer disk edge (e.g., Pringle 1981). We used the alpha approach from Shakura & Sunyaev (1973) to approximate the kinematic viscosity: $$\nu =\alpha \Omega H^2,$$(2)

where α is a parameter that determines the efficiency of viscous transport, $\Omega ={(G{M}_{*}/r^3)}^{1/2}$ is the Keplerian angular velocity, and H is the scale height of the gas disk. We calculated the scale height as $$H=c_{\text{s}}/\Omega ,$$(3)

where c_s is the sound speed, $$c_{\text{s}}={\left(\frac{k_{\text{B}}T}{\mu m_{\text{H}}}\right)}^{1/2}.$$(4)

In the equation above, k_B is the Boltzmann constant, T is the temperature, μ is the mean molecular weight, and m_H is the mass of the hydrogen atom. The mean molecular weight was set to be 2.34, corresponding to a solar–composition mixture of hydrogen and helium (Hayashi 1981). The midplane temperature of the disk was approximated using a fixed powerlaw structure: $$T=150\text{K}\times {(r/\text{AU})}^{-3/7},$$(5)

where 150 K is the temperatureat 1 AU (Chiang & Goldreich 1997).

2.1.2 Planetary gaps

Massive planets push away material from the vicinity of their orbits, and as a result opens up a gap in the disk. We modeled these planetary gaps using a simple approach with Gaussian gap profiles. The Gaussian is described by the equation $$G(r)=\frac{{\Sigma }_{a,0}}{{\Sigma }_{a,\text{min}}}\exp \left[-\frac{{(r-a)}^2}{2H_a^2}\right],$$(6)

where Σ_a,0 is the unperturbed surface density at the location of the planet, Σ_a,min is the surface density at the bottom of the gap, and a is the semimajor axis of the planetary orbit. The depth of the gap is calculated as $$\frac{{\Sigma }_{a,\text{min}}}{{\Sigma }_{a,0}}=\frac{1}{1+0.04K},$$(7)

where K is given by $$K=q^2{\left(\frac{H_a}{a}\right)}^{-5}{\alpha }^{-1},$$(8)

and q is the planet to star mass ratio (e.g., Kanagawa et al. 2015). Each planet contributes their own Gaussian, and the final surface density profile is then obtained by dividing Eq. (1) by $1+G{(r)}_{\text{1}}+G{(r)}_{\text{2}}+\mathrm{...}$, where $G{(r)}_{\text{1}}$ is the Gaussian of Planet 1 (etc.).

2.1.3 Gas density ρ

We use the following expression to obtain the midplane gas density, $$\rho (z=0)=\frac{\Sigma }{\sqrt{2\pi }H}.$$(9)

In order to find the gas density at some height, z, away from the midplane, we assume vertical hydrostatic equilibrium for the gas in the disk. We then end up with the equation: $$\rho (z)=\rho (z=0)\exp \left[\frac{-z^2}{2H^2}\right].$$(10)

2.2 Drag force

Planetesimals that are moving through the disk at a velocity different from that of the gas either experience a headwind or a tailwind and, as a result, they are decelerated or accelerated towards the value of the gas velocity. This change in velocity occurs on a timescale called the stopping time, t_s. Assuming that the planetesimals are spherical, their stopping time is $$t_{\text{s}}={\left(\frac{\rho }{{\rho }_{•}}\frac{v_{\text{th}}}{R_{\text{pl}}}\min \left[1,\frac{3}{8}\frac{v_{\text{rel}}}{v_{\text{th}}}{C}_D(Re)\right]\right)}^{-1},$$(11)

where ρ_∙ and R_pl are the solid density and radius of the planetesimals, v_rel is the relative velocity between the planetesimals and the gas, $v_{\text{th}}=\sqrt{8/\pi }{c}_{\text{s}}$ is the gas thermal velocity, C_D is the dimensionless drag coefficient, and Re is the Reynolds number (e.g., Perets & Murray-Clay 2011; Guillot et al. 2014). The dimensionless drag coefficient is calculated as (Perets & Murray-Clay 2011): $$C_D=\frac{24}{Re}{(1+0.27Re)}^{0.43}+0.47\left(1-\exp \left[-0.04R{e}^{0.38}\right]\right)$$(12)

and the Reynolds number is $$Re=\frac{4R_{\text{pl}}{v}_{\text{rel}}}{c_{\text{s}}{l}_{\text{g}}},$$(13)

where l_g ~ 5 × 10⁻⁶ kg m⁻³∕ρ is the mean-free path of the gas (Supulver & Lin 2000).

The acceleration of the planetesimal due to gas drag is calculated as $$a_{\text{drag}}=-\frac{1}{t_{\text{s}}}(v_{\text{pl}}-v_{\text{gas}}),$$(14)

where v is a velocity vector (e.g., Whipple 1972). We assume that the gas velocity in the z-direction is zero and obtain the Cartesian components of the gas velocity by projecting the orbital velocity of the gas, $$v_{\varphi ,\text{gas}}\approx v_{\text{K}}+\frac{1}{2}\frac{c_{\text{s}}^2}{v_{\text{K}}}\frac{\partial \ln P}{\partial \ln r},$$(15)

in the x and y plane. In the above equation, $v_{\text{K}}=\sqrt{G{M}_{*}/r}$ is the Keplerian orbital velocity and ∂lnP∕∂lnr is the radial gas pressure gradient. The radial gas pressure gradient is calculated as ∂ ln P∕∂lnr = ∂lnΣ∕∂lnr + ∂lnT∕∂lnr − ∂lnH∕∂lnr, where ∂ ln Σ∕∂lnr = 15∕14.

2.3 Mass ablation

If a planetesimal becomes sufficiently heated, then the material at its surface can undergo phase transitions, resulting in mass loss. Here it is assumed that the mass loss occurs due to solid material transitioning to gas phase by frictional heating and irradiation from the surrounding gas through a process called ablation. In order to estimate the mass ablation rate from a planetesimal surface, we first need to know its composition. We assume that all planetesimals are non-differentiated and that they consist of a mixture of silicate grains, carbon grains, and volatile ices. In this paper, we consider the three volatiles H₂O, CO₂, and CO, and the volatile content of a planetesimal is set by which volatile ices were present at its formation site. The total ablation rate from a planetesimal surface is then taken to be the sum of the ablation rates for each present volatile ice and since the silicate and carbon grains are well mixed with the ices, we assume them to be released along with the ices (i.e., we assume no crust formation; see Appendices A and B for a discussion of this assumption). For a planetesimal that forms in a region of the disk where H₂O and CO₂ are in solid form, but CO is in gas form, the total ablation rate will thus be calculated as: $${\dot{m}}_{\text{abl}}={\dot{m}}_{\text{abl},{\text{H}}_2\text{O}}+{\dot{m}}_{\text{abl},{\text{CO}}_2}.$$(16)

We follow Ronnet & Johansen (2020) and use the following expression for the ablation rate of an element X: $${\dot{m}}_{\text{abl},X}=-4\pi R_{\text{pl}}^2{P}_{\text{sat},X}(T_{\text{pl}})\sqrt{\frac{{\mu }_X}{2\pi R_{\text{g}}{T}_{\text{pl}}}},$$(17)

where P_sat,X is the saturated vapor pressure of element X, T_pl is the surface temperature of the planetesimal, μ_X is the molecular weight of element X, and R_g is the ideal gas constant (e.g., D’Angelo & Podolak 2015). Expressions for the saturated vapor pressure as polynomials of the temperature are given in Fray & Schmitt (2009) for all three volatile ices under consideration. These polynomial expressions are only accurate above a certain temperature, which varies depending on the element. We therefore introduce “floor values” to the saturated vapor pressure. A plot of P_sat,X versus T_pl for all three volatile ices, with the floor values included, is presented in Fig. C.1.

The equilibrium surface temperature of the planetesimals is obtained using the following equation from Ronnet & Johansen (2020): $$T_{\text{pl}}^4=T^4+\frac{c_D\rho v_{\text{rel}}^3}{32{\sigma }_{\text{sb}}}-{\sum }_X\frac{P_{\text{sat},X}(T_{\text{pl}})}{{\sigma }_{\text{sb}}}\sqrt{\frac{{\mu }_X}{8\pi R_{\text{g}}{T}_{\text{pl}}}}{L}_X,$$(18)

where σ_sb is the Stefan-Boltzmann constant and L_X is the latent heat of vaporization of element X. We calculated the latent heat of vaporization using the Clausius-Clapeyron relation for low temperatures and pressures, and obtained the following results: $L_{{\text{H}}_2\text{O}}=2.8\times {10}^6\text{J}{\text{kg}}^{-1}$, $L_{{\text{CO}}_2}=6.1\times {10}^5\text{J}{\text{kg}}^{-1}$, and L_CO = 2.8 × 10⁵ J kg⁻¹. The second term in Eq. (18) reflects heating due to gas friction and the third term is cooling due to the ablation of volatile ices. Each volatile ice contributes its own cooling term, meaning that the planetesimal temperature will vary depending on what volatile ices it consists of. Since the cooling depends itself on the surface temperature of the planetesimal, the equation has to be solved iteratively. This is done using a bisection method. A plot of the mass ablation rate as a function of the semimajor axis is presented in Fig. D.1 for all three volatile ices and for different surface temperatures of the planetesimals.

2.4 Condensation temperature

We assume an isothermal equation of state for the disk pressure of an element X: $$P_{\text{d},X}=c_{\text{s},X}^2{\rho }_X.$$(19)

The sound speed of X is calculated using Eq. (4), and exchanging the mean molecular weight of the disk with the molecular weight of X. The density of X in the disk midplane is obtained by multiplying the total density in the midplane (Eq. (9)) with the mass fraction of X with respect to the disk. To obtain the relevant mass fractions, we use abundances from Öberg et al. (2011), and for the calculation, we assume that 25% of the disk mass is in He and the remaining 75% is in H. This results in the following mass fractions: 1.2 × 10⁻³, 9.7 × 10⁻⁴, and 3.1 × 10⁻³ for H₂O, CO₂ and CO. The mass fractions for silicate grains and carbon grains are 3.6 × 10⁻³ and 5.3 × 10⁻⁴.

The condensation temperature of element X is found by comparing Eq. (19) with the saturated vapor pressure and solving for the temperature. Since the density of the disk depends on Ṁ₀, α and r_out and as these are parameters that are varied among the simulations, the condensation temperatures will not be the same in all simulations. This means that a planetesimal which is initiated at the same semimajor axis in two separate simulations can have different compositions.

2.5 Collision timescale

We do not consider planetesimal-planetesimal collisions in our simulations; however, we still want to obtain an estimate of how common they are. If they happen on a timescale that is longer than the actual simulated time, it is a justified choice to neglect them. However, if the opposite is true, we need to consider how they would have affected our results.

Johansen & Bitsch (2019) provide an expression for the collision timescale of planetesimals in the gravitationally unfocused case: $$t_{\text{coll}}=\frac{1}{n_{\text{pl}}{\sigma }_{\text{pl}}\delta v}\approx \frac{R_{\text{pl}}{\rho }_{•}}{{\Sigma }_{\text{pl}}\Omega }=\frac{R_{\text{pl}}{\rho }_{•}}{2\delta v{\rho }_{\text{pl}}},$$(20)

where n_pl is the number density, σ_pl is the physical cross section, δv is the relative speed, and ρ_pl = Σ_pl∕(2H_pl) is the volume density of planetesimals. The relative speed between the planetesimals is approximated as (Lissauer & Stewart 1993): $$\delta v={\left(\frac{5}{4}{e}^2+i^2\right)}^{1/2}{v}_{\text{K}}.$$(21)

The scale height of the planetesimals is taken to be H_pl = i × a, where a is the semimajor axis.

3 Numerical setup

We used the N-body code REBOUND to perform our simulations and modified it to take into account the effect of gas drag and mass ablation (Rein & Liu 2012). The simulations were executed using the hybrid symplectic integrator MERCURIUS and the timestep was set to be one twentieth of the innermost planet’s dynamical timescale (Rein et al. 2019). Additional simulations with smaller timesteps were performed as well, in order to check that the outcome was not affected. We added the planets and the central star as active particles, and the planetesimals as semi-active particles with a mass. Semi-active particles only interact gravitationally with active particles. Collisions between active particles and semi-active particles were detected and recorded using a direct search method, and resultedin perfect merging. If a particle were to leave the simulation domain, it is recorded and removed from the simulation.

We considered two planetary systems that would be representative of the young Solar and HL Tau systems. In the Solar System simulations, the protoplanetary disk stretches from 0.1 to 100 au with R_out = 20 au, and it is modeled using a linear grid with 1000 grid cells. The location of a particle on this grid is found using a binary search algorithm. The simulation box is centered on the sun and stretches 100 au in x and y-direction, and 20 au in z-direction. In the HL Tau simulations, instead we use a disk with 2000 grid cells that stretches from 0.5–200 au and has R_out = 100 au. The corresponding simulation box is 500 au in x and y-direction and 100 au in z-direction.

We use two planets in the Solar System simulations (Jupiter and Saturn), and three planets in the HL Tau simulations (placed at the locations of the major gaps in the disk). We do not consider planet growth or migration. Jupiter and Saturn are initiated with their current eccentricity and inclination, and we use their current bulk density to calculate the planetary radius, given the masses in Table 1. The three planets in HL Tau are initiated with close to zero eccentricity and inclination, and their radius is calculated using the masses in Table 1 and assuming a constant density of 1000 kg m⁻³. We use a central star of solar mass and solar luminosity in the HL Tau simulations.

The planetesimals are initiated uniformly between 1 and 2 Hill radii away from the planets, in the direction away from the central star. This is roughly the region in which planetesimals form in Paper I and additional simulations show that small changes to this formation location does not affect the results. A study on how the simulation outcomes are affected by larger changes to the planetesimal formation location is presented in Appendix E. Generally, as long as the planetesimals do not form further away than about 5 Hill radii from the planets, the results do not change significantly. We initiate 50 planetesimals beyond each planetary gap, meaning that each individual Solar System simulation harbors 100 planetesimals and each individual HL Tau simulation harbors 150. In order to provide better statistics, we performed 10 simulations for each set of parameter that we study, so that the total number of planetesimals per parameter set amounts to 1000 and 1500 for the Solar System and HL Tau system simulations, respectively.

The planetesimals are given an initial radius of 100 km, consistent with constraints from Solar System observations (Bottke et al. 2005; Morbidelli et al. 2009) and streaming instability simulations (e.g., Johansen et al. 2015), and have a constant solid density of 1000 kg m⁻³. We add a property to the planetesimals that is the temperature of the disk at their formation location. This temperature, in relation to the condensation temperatures of the volatile ices, determines the composition of the planetesimals. The effect of gas drag is added as a velocity dependent force and mass ablation is added as a post-timestep modification. We keep track of how much mass is ablated from the planetesimals in each radial bin and at what time during the simulation. If a planetesimal loses more than 99% of its original mass, it is removed from the simulation and its remaining mass is considered to be ablated at the time and location where it was removed. The same procedure is applied if the amount of mass ablated in one timestep is larger than the total remaining mass of the planetesimal.

We performed a parameter study in order to explore how the fate of the planetesimals formed at planetary gap edges is affected by: (1) their composition; (2) the mass of the planets; and (3) the density of the disk (controlled by the disk parameters α and Ṁ₀). The parameter values used in the different simulations can be found in Table 1. In simulation #1, which we will hereafter refer to as the nominal Solar System simulation, the masses of Jupiter and Saturn were set, respectively, to 90 and 30 M_⊕, α = 10⁻² and Ṁ₀ = 10⁻⁷ M_⊙ yr⁻¹. In such a disk, the CO₂ iceline is located interior to the orbit of Jupiter and the CO iceline is located well beyond the orbit of Saturn, meaning that all planetesimals in this simulation contain H₂O and CO₂ ice. The nominal HL Tau simulation is labeled #7 in Table 1. In this simulation, the planetary masses where set to equal the pebble isolation mass, calculated using an analytical fitting formula from Bitsch et al. (2018). As in the nominal Solar System simulation, these planetesimals contain H₂O and CO₂ ice.

4 Simulations of the Solar System

In the SolarSystem simulations, we include two planets, Jupiter and Saturn, which are placed at their current semimajor axes. The mass ratio between Jupiter and Saturn is 3:1 in all simulations, which is roughly the current value, but the actual masses are varied. Results from the nominal model are presented in Sect. 4.1 and discussed in detail. We examine how varying the parameters affects the results in Sect. 4.2.

4.1 Nominal model

In the nominal model, the mass of Jupiter is set to be 90 M_⊕ and the mass of Saturn is set to be 30 M_⊕. All parameters used in the nominal model, also called simulation #1, can be found in Table 1. The planetesimals are initiated in a narrow region just beyond the orbits of Jupiter and Saturn, as suggested by the streaming instability simulations in Paper I. At these distances from the Sun, and with the disk parameters stated in Table 1, H₂O and CO₂ are in solid form, while CO is in gas phase. We therefore consider these planetesimals to have a volatile content of H₂O and CO₂ and we consider the ablation of both these molecules.

4.1.1 Dynamical evolution

In this section, we address the question of the dynamical redistribution of the planetesimals formed at planetary gap edges. In Figs. 1 and 2, we show the eccentricity and semimajor axis evolution for the 1000 planetesimals in the nominal simulation. The planetesimals are initiated with an eccentricity close to zero, but from Fig. 1, we see that the close proximity to the planets leads to rapid scattering, resulting in eccentricities as high as 0.4 already after 1 kyr of evolution. Most planetesimals that are not immediately scattered interior to the orbit of Jupiter end up in the scattered disk of Jupiter or Saturn. The eccentricities within these scattered disks increases with time, as the planetesimals are continuously scattered at perihelion. Such large eccentricities lead to large velocities at perihelion (v_peri) relative to the gas, $$v_{\text{peri}}=v_{\text{K}}\sqrt{\frac{1+e}{1-e}},$$(22)

which, in turn, affects the thermodynamic evolution of the planetesimals. The velocity at perihelion is higher within Jupiter’s scattered disk than it is within Saturn’s for the same orbital eccentricity since the Keplerian velocity decreases with increasing semimajor axis.

In Fig. 2, we separate the planetesimals forming at the gap edge of Jupiter and Saturn into different panels, making the diffusionof the semimajor axes clearly visible. Strong scatterings causes the semimajor axes to diffuse over several au in less than 100 yr. After a few thousand years, planetesimals from both gap edges are spread out across the giant planet region and the ones closest to the star have semimajor axes of ~ 3 au. A small number of planetesimals, 13∕1000 in this particular simulation, suffer a collision with either Jupiter or Saturn. Most of these collisions are between Jupiter and planetesimals formed at Jupiter’s gap edge. Concerning the planetesimals that obtain very large eccentricities and semimajor axes, many of them are eventually scattered outside the simulation domain, and are considered to have been ejected from the system. In this simulation, ~15% of all planetesimals are eventually ejected.

In summary, planetesimals that form at the edges of planetary gaps do not remain at their initial birth location. The close proximity to the gap-opening planets results in strong scatterings, causing the planetesimals to spread out across the protoplanetary disk. Many planetesimals initially become members of the gap-opening planet’s scattered disk, while some are scattered towards the inner disk region. The resulting high eccentricities lead to a high velocity relative to the gas, which has substantial implications for the mass evolution, as shown in Sect. 4.1.2. In Appendix E, we present simulations where the planetesimals are formed further from the planet. Placing planetesimals further from the planet does not change our results qualitatively, but it does lead to a delayed and slower scattering phase.

Fig. 1 Eccentricity, semimajor axis, and surface temperature evolution for the 1000 planetesimals in the nominal Solar System simulation. The presented orbital parameters have been averaged over 100 yr and the surface temperatures are the maximum values during the same time period, which roughly corresponds to the surface temperatures at perihelion. The temperature of the surrounding gas (Tdisk) is shown as a colorbar in the left panel. The solid black lines mark the orbits with a perihelion corresponding to Jupiter’s and Saturn’s location. The planetesimals which cluster around these lines are (at least momentarily) members of Jupiter’s or Saturn’s scattered disk. Planetesimals that are scattered interior of Jupiter’s orbit obtain high surface temperatures at perihelion and experience efficient ablation.

Fig. 2 Semimajor axis and surface temperature evolution for the 1000 planetesimals in the nominal Solar System simulation (same data as in Fig. 1). For simulation times longer than 1 kyr, we average the semimajor axes over 100 yr and show the maximum surface temperatures during the same time period; whereas for shorter simulations times, we show non-averaged values. Upper panel: planetesimals formed at the gap edge of Jupiter and lower panel: planetesimals formed at the gap edge of Saturn. The close proximity to the giant planets result in continuous scatterings, which causes the fast semimajor axis diffusion that is displayed in the plot.

4.1.2 Mass evolution

In this section, we look into the mass evolution of planetesimals; that is, the loss of mass due to ablation. The ablation rate depends strongly on the planetesimal surface temperature, which increases towards the central star as the disk temperature and gas density becomes higher. The surface temperature is also highly dependent on the velocity relative to the gas, which increases with increasing orbital eccentricity according to the process described in Sect. 4.1.1. Based on this, the highest surface temperature, and thus the ablation rates, should be obtained at the perihelion passage of eccentric orbits close to the star. This can be seen in Figs. 1 and 2, where we also show the surface temperature evolution of the planetesimals at perihelion. The planetesimals that are scattered interior to Jupiter’s orbit obtain surface temperatures around 100 K at perihelion. At these temperatures, CO₂ ablation quickly sublimates the planetesimal, which is why there are no planetesimals in the innermost part of the disk. We note here that the ablation process itself severely decreases the surface temperatures through cooling due to the latent heat of vaporization. A similar simulation without taking into account ablation would result in surface temperatures that are many times higher (see Appendix A). We also note that a planetesimal with the same surface temperatureand orbital parameters in a gas-free disk would not experience any ablation due to the lack of frictional heating.

In Fig. 3, we track the mass loss of six selected planetesimals from the nominal simulation, as they are scattered around by the planets. The time evolution of the semimajor axes, perihelia, aphelia, and eccentricity for the same planetesimals can be found in Fig. F.1. The planetesimals that form at Jupiter’s gap edge (the first three legend entries) begin to experience mass loss already after 100 yr of evolution. The “red” planetesimal is a member of Jupiter’s scattered disk for a few thousand years, until it is kicked towards Saturn’s scattered disk, where it remains until the end of the simulation. Since its orbit never enters the inner disk region, ablation is slow and the planetesimal only loses about 2% of its mass. The “yellow” planetesimal obtains multiple kicks by Jupiter during the first few thousand years, after which it ends up interior to Jupiter’s orbit with a perihelion of just above 2 au and becomes completely ablated within 10 kyr. The “green”planetesimal is a member of Jupiter’s scattered disk for 100 kyr, with a relatively small eccentricity, until a strong planetary encounter leaves it on an orbit interior to Jupiter, where it very quickly becomes ablated.

The planetesimals formed at Saturn’s gap edge (the last three legend entries), begin to lose mass much later. The “light-blue” planetesimal sits on a low-eccentric orbit in between Jupiter and Saturn for 150 kyr, after which it obtains a strong kick and becomes a high-eccentric member of Jupiter’s scattered disk. Its eccentricity continues to increase until it leaves the simulation domain, having lost 85% of its mass. The “dark-blue” planetesimal is scattered onto a 4 au orbit by Saturn, where it slowly becomes circularized. Once circularized, the relative velocity between the planetesimal and the gas turns to zero, and it remains in the same orbit until the end of the simulation. Finally, the “purple” planetesimal never enters the region interior to Saturn, but remains on a low-eccentric orbit in the outer part of the disk for the entire simulation.

The mass loss tracks presented in Fig. 3 are examples of what can happen to a planetesimal in the simulation. In Fig. 4, we present the mass loss and semimajor axis evolution for all planetesimals in the nominal simulation. The planetesimals formed at Jupiter’s gap edge start to lose mass much earlier than those formed at Saturn’s gap edge and they generally lose mass at a faster rate. This is very much expected, mostly because the planetesimals form in a warmer part of the disk, but also because the mass of Saturn is lower than that of Jupiter and, therefore, the planetesimals are not scattered as much. Towards the end of the simulation, ~ 50% of all planetesimals have become completely ablated. In Appendix E we show that placing the planetesimals further from the planet results in a delayed and slower ablation phase.

From the plots presented here, it is evident that mass loss due to ablation plays a major role in the evolution of planetesimals formed at planetary gap edges, at least for the parameters used in the nominal model. In Sect. 4.2, we will investigate exactly how much mass is lost and where and compare that to simulations with varying planetary masses and disk parameters.

Fig. 3 Mass loss versus semimajor axis for 6 selected planetesimals from the nominal Solar System simulation. The mass loss is calculated as 1 − M(t)∕M(t = 0). Filled circlesmark 100 kyr of evolution and the formation location of the planetesimals (shown at the bottom of the plot) has been marked “ST”. The time before the planetesimals first experience mass loss, that is, the time before they appear on the plot, has been included in the legend. The first three legend entries are for planetesimals formed at the gap edge of Jupiter and the following three are for planetesimals formed at the gap edge of Saturn. The dotted lines mark the semimajor axis of Jupiter and Saturn. Planetesimals which are scattered interior of Jupiter’s orbit lose mass at a high rate, while those which are scattered exterior of Saturn’s orbit experience little mass loss.

4.1.3 Planetesimal–planetesimal collisions

In this section, we use simple calculations to estimate how common planetesimal–planetesimal collisions are. If they occur on the same timescale as dynamical scattering, then they could affect our results. Furthermore, if these collisions are strong enough to disrupt the planetesimals that are involved, it would constitute another mechanism to replenish dust and pebbles in the disk.

We splitthe disk into multiple semimajor axis rings and find the planetesimals which are located within each ring. We then use the average eccentricity and inclination of those planetesimals to get an estimate of the velocity dispersion (Eq. (21)) and the scale height.To calculate the collision timescale, we further need to know the volume density of planetesimals within each ring. If all solids between the semimajor axis of Jupiter and Saturn are converted into planetesimals at Jupiter’s gap edge, then there would be 16 M_⊕ of planetesimals at that location, based on the disk parameters of the nominal model and a dust-to-gas ratio of 1%. We assume that the amount of planetesimals forming at Saturn’s gap edge is the same. The total planetesimal mass in a ring is then simply taken to be the sum of the mass of all planetesimals in that ring.

When using a ring-width of 2 au, this resulted in collision timescales above 1 Myr, at all locations in the disk and at all times during the simulation (see Fig. 5). Decreasing the ring-width to 1 au leads to slightly smaller collision timescales. These results indicate that an individual planetesimal suffer a very low risk of colliding with another planetesimal. Since, in reality, there are millions of planetesimals forming at the gap edges, collisions will still be occurring, but not often enough to affect any of our results.

4.2 Parameter study and mass loss profiles

In this section, we vary certain disk and planet parameters in order to find out how the simulation results are affected. The parameters used in each simulation can be found in Table 1. In Fig. 6, we show how the total amount of ablated mass varies across the disk (left panel) and as a function of time (right panel). The mass loss evolution for each simulation is shown in Fig. 7, where the number of planetesimals that have become completely ablated or scattered or that have suffered planetary collisions can be noted and compared as well. The collision timescale for all simulations in the parameter study are very similar to those in the nominal simulation, meaning that planetesimal-planetesimal collisions can be safely ignored as a means of producing dust.

In the nominal simulation, 70% of the initial planetesimal mass has been ablated at the end of the simulation. The majority of this mass is released either just interior of Jupiter’s orbit, between 2 and 5 au, or just outside of it, between 6 and 8 au. The mass loss outside of Jupiter’s orbit is possible because of the large orbital eccentricities that result in high surface temperatures. The planetesimals in the nominal simulation all form beyond the CO₂ iceline, and thus both H₂O and CO₂ ablation is considered. In simulation #2 we remove the CO₂ ablation and consider the volatile content of the planetesimals as only made up of H₂O. As a result, the amount of ablated mass at the end of the simulation drops from 70% to 30%. Furthermore, since the ablation rate of H₂O is lower than that of CO₂ for the sametemperature (see Fig. D.1), the planetesimals need to be scattered further towards the Sun in order for efficient ablation to occur. This can be seen in the left panel of Fig. 6, where the peak of the ablation curve interior to Jupiter has been shifted closer to the Sun. Additionally, there is no longer any mass released beyond the orbit of Jupiter, telling us that all mass released in this region in the nominal simulation is attributed to the presence of CO₂.

In simulation #3, we decrease the planetary masses by a factor of 3, mimicking an earlier stage of the planet formation process. The mass of Jupiter is now equal to the pebble isolation mass, while the mass of Saturn is only 20% of the pebble isolation mass. In Paper I, the amount of planetesimals that form at the gap edges drops by 80% when the mass is decreased from one pebble isolation mass to half a pebble isolation mass. The amount of planetesimals forming beyond Saturn in this simulation would thus have been very small, if any at all, and so we only include planetesimals at Jupiter’s gap edge. The number of planetesimals in the simulation is kept the same as in all other simulations, but since they represent only half the mass, we divide the amount of ablation by a factor of 2, with the results provided in Fig. 6.

Lowering the planetary masses and removing the planetesimals at Saturn’s gap edge result in about 40% less mass ablation than in the nominal simulation. The removal of planetesimals at Saturn’s gap edge results in fewer planetesimals far out in the disk. Lowering the planetary masses results in weaker planetary scatterings, which, in turn, has several effects: (1) only 1/1000 planetesimals are ejected beyond the simulation domain; (2) the orbits are not as excited, resulting in lower surface temperatures and slower ablation; (3) the planetesimals end up on orbits closer to their birth locations, and thus the mass loss is more concentrated towards Jupiter’s location.

In simulation #4, the planetary masses are increased, mimicking a later stage in the planet formation process. The stronger planetary scatterings initially result in higher mass loss rates as the planetesimals orbits quickly become excited. However, the strong planetary scatterings also lead to many more planetesimals being ejected beyond the simulation domain, which leaves fewer planetesimals in the system to be ablated during later times.

In the final two simulations of the parameter study, simulations #5 and #6, the viscosity parameter and disk accretion rate are decreased by a factor of 10. Decreasing the viscosity parameter slightly changes the gap profiles, but mostly it results in a surface density that is ten times larger, thus mimicking an earlier stage of disk evolution. The ice lines are also shifted closer to the Sun, but in this case, that does not change the composition of the planetesimals compared to the nominal model. On the other hand, decreasing the disk accretion rate results in a surface density that is ten times lower, causing the ice lines to move further out in the disk. As a result, the CO₂ ice line ends up beyond the gap edge of Jupiter and, thus, the planetesimals forming there do not contain any CO₂ ice.

Looking at Fig. 6, we see that lowering the viscosity parameter results in more ablation, while lowering the disk accretion rate results in less ablation, just as expected. The higher surface density in simulation #5 results in higher surface temperatures, allowing for planetesimals to lose mass further out in the disk and at a higher rate. This effect in a combination with larger gas friction lead to a result whereby there are almost no planetesimals being ejected from the system. The opposite is the case in simulation #6 and since the planetesimals formed at Jupiter’s gap edge do not contain CO₂ ice, there is little mass loss exterior of Jupiter’s orbit, just as in simulation #2.

The results from the parameter study show that the surface density of the disk is a key parameter in determining how much mass is ablated from the planetesimals. Since the surface density in disks decreases with time, ablation is expected to be much more efficient in young disks than in old ones. Thus, planetesimals forming at the gap edge late during the disk lifetime might not suffer any significant mass loss due to ablation, but could, for example, be implanted into the asteroid belt or aid in water delivery to Earth (Raymond & Izidoro 2017). The masses of the planets determine how excited the planetesimal orbits become and how far they are scattered within the disk. To some extent, the planetary mass also affects how many planetesimals form at the gap edge (see Paper I). Finally, the composition of the planetesimals plays a major role in determining the efficiency of the ablation process. This is set by the location of formation relative to the location of the major ice lines, which, in a real disk, shifts inwards over time as the temperature of the disk decreases.

Fig. 4 Mass loss versus semimajor axis evolution for the 1000 planetesimals in the nominal Solar System simulation (same data as in Figs. 1 and 2). The thin grey lines mark the perihelion and aphelion of the planetesimal orbits. The number of planetesimals formed at Jupiter’s (red dots) respective Saturn’s (blue dots) gap edge which do not appear on the plot, because they have been either: completely ablated; collided with a planet; or ejected beyond the simulations domain, is written in each panel. Planetesimals formed at the gap edge of Jupiter generally experience more ablation than planetesimals formed at the gap edge of Saturn. About 50% of all planetesimals have become completely ablated after 500 kyr.

Fig. 5 Time evolution of the timescale for planetesimal-planetesimal collisions for the nominal Solar System simulation. The collision timescale was calculated using Eq. (20) and assuming an initial planetesimal mass of 16 M⊕ per gap edge. The dotted line marks where the collision timescale equals the time of the simulation. Since the resulting collision timescale for one planetesimal is several orders of magnitude larger than the actual simulated time, planetesimal–planetesimal collisions can be safely ignored in our simulations and will not constitute a significant path to replenishing the dust component in the disk.

Fig. 6 Distribution of ablated mass as a function of semimajor axis (left) and time (right) for all Solar System simulations. Left plot: Distribution of ablated mass across the disk after 500 kyr. The values on the y-axis represent the amount of mass that has been ablated in a 0.2 au semimajor axis bin. The dotted lines mark the semimajor axes of Jupiter and Saturn. Most mass loss occur in the region just interior to and exterior to Jupiter’s orbit. Right plot: Total amount of mass that has been ablated as a function of time, for the same data as in the left plot. The colored lines show the average over the ten simulations, and the colored region shows the one standard deviation away from this value. The dotted black lines show the best fit to the curve f(x); and the parameters to the fits can be found in Table 2. The mass ablation rate is highly dependent on the gas surface density, where a high surface density results in a high degree of ablation (#5) and a low surface density results in low degree of ablation (#6).

5 Simulations of HL Tau

We include three planets in the HL Tau simulations and place them at the locations of the major gaps in the disk, that is, at 11.8, 32.3, and 82 au (Kanagawa et al. 2016). The masses of the planets are set to be 59.6, 141.2, and 313.7 M_⊕, which is exactly equal to the pebble isolation mass using the parameters in simulation #7. For reference, this is the same set-up as in simulation #4 of Paper I. Each simulation contains 150 planetesimals, with 50 formed at each planetary gap edge, and we run ten simulations per parameter set. Results from the nominal simulation are presented in Sect. 5.1, while in Sect. 5.2, we study how these results change when we vary the planetesimal compositions and the disk parameters.

5.1 Nominal model

The three planets in the HL Tau disk are located well beyond the CO₂ iceline, which in the nominal model (simulation #7) is at 4.7 au. The CO iceline sits much further out in the disk at 99.3 au, placing it just beyond the formation location of the outermost planetesimal. Since all planetesimals in the nominal model form in between the CO₂ and CO iceline, we thus consider them to be made up of H₂O and CO₂ ice and we consider the ablation of these two molecules.

The dynamical evolution of the 1500 planetesimals from the nominal simulation is presented in Figs. 8 and 9. Just as in the Solar System simulations, most planetesimals end up in the scattered disks of the planets, with eccentricities increasing over time. The two outermost planets are very massive and deliver strong kicks to the planetesimals formed in their vicinity, causing them to quickly spread out over the entire disk. Many of these planetesimals are eventually ejected beyond the simulation domain. The planetesimals that form at the gap edge of the innermost planet instead suffer weaker kicks and tend to remain in the inner disk region.

The timescale for planetesimal-planetesimal collisions is calculated in the same way as for the Solar System simulations and we use the results from Paper I to infer the total mass of planetesimals formed at each gap edge. In simulation #4 of Paper I, after 1 Myr, 281 M_⊕ of planetesimals have formed. We make the assumption that the amount of planetesimals forming at each gap edge is the same and based on the results of Paper I, this is a reasonable approximation. This yields a collision timescale of ~ 10⁴ yr in the beginning of the simulation, which then increases roughly linearly with time as the planetesimals spread out across the disk (see Fig. 10). Since the collision timescale remains an order of magnitude larger than the actual simulated time throughout the simulation, for HL Tau we also conclude that planetesimal-planetesimal collisions are not a major source of dust production

The amount of mass that has been ablated as a function of semimajor axis and time for the nominal model is presented in Fig. 11. From the left plot, we learn that all mass loss occurs interior of the innermost planet, with a peak at ~ 7 au. In total, 11% of the initial planetesimal mass has been ablated after 1 Myr. This is a much smaller number than what we found in the Solar System simulations and the main reason for this is that the planetesimals form further out in the disk. In order for a planetesimal forming at 90 au to experience efficient ablation, it needs to be scattered inwards by approximately 80 au. Since the planetesimals need to travel a further distance before mass loss can begin to occur, the initial mass loss rate is much lower than in the Solar System simulations, where planetesimals already start to ablate efficiently after a few thousand years. This can also be seen in the top panel of Fig. 12, where the mass loss evolution for all planetesimals in the ten nominal simulations has been plotted. After 10 kyr, the planetesimal with the smallest semimajor axis has still only lost about 30% of its mass.

Fig. 7 Mass loss versus semimajor axis evolution for all Solar System simulations. The masses of the planets affect how far the planetesimals are scattered, with more planetesimals leaving the simulation domain when the planetary mass is increased (#4). A large surface density (#5) generally results in efficient ablation and few scatterings beyond the simulation domain since the planetesimals are ablated before they become scattered. The opposite is the case for a low surface density (#6) and ablation also occurs closer to the star.

Fig. 8 Eccentricity, semimajor axis, and surface temperature evolution for the 1500 planetesimals in the nominal HL Tau simulation, produced in a similar manner as Fig. 1. The solid black lines mark the perihelia of the planets. Many planetesimals end up in the scattered disks of the planets, where they obtain high eccentricities. Only planetesimals scattered tointerior of the innermost planet’s orbit obtain perihelion temperatures large enough for ablation to take place.

Fig. 9 Semimajor axis and surface temperature evolution for the 1500 planetesimals in the nominal HL Tau simulation, produced in a similar manner as in Fig. 2 (same data as in Fig. 8). Planetesimals formed at the gap edge of the innermost planet is shown in the upper panel. For the middle and outermost planet, these are shown in the middle and bottom panels, respectively. Since the planetary masses are set to equal the pebble isolation mass, and the pebble isolation mass increases with semimajor axis, the planetesimals formed at the gap edge of the outermost planet experience stronger scatterings and diffuse faster than those formed at the innermost planet’s gap edge.

Fig. 10 Time evolution of the timescale for planetesimal–planetesimal collisions for the nominal HL Tau simulation. The collision timescale was calculated using Eq. (20) and using an initial planetesimal mass of 281 M⊕, following the results from simulation #4 of Paper I. The dotted line marks where the collision timescale equals the time of the simulation. Even though the collision timescale for one planetesimal can be as low as 104 yr, it is still roughly an order of magnitude larger than the actual simulated time, meaning that planetesimal-planetesimal collisions can be safely ignored in our simulations as a source of dust production.

5.2 Parameter study and mass loss profiles

In this parameter study, we chose to vary the composition of the planetesimals, the α value, and the initial mass accretion rate Ṁ₀. Changing these disk parameters does not result in a significant change to the collision timescale compared to the nominal simulation. As previously explained, lowering α results in an increased surface density, while lowering Ṁ₀ results in a decreased surface density. This change in the surface density also affects the positions of the icelines. Lowering Ṁ₀ by a factor of 10 (simulation #11) does not result in any change of the planetesimal compositions compared to the nominal simulation; however, when α is lowered by a factor of 10, then the planetesimals forming at the gap edge of the outermost planet end up interior of the CO iceline. This gives us an opportunity to study CO ablation since the planetesimals will consist of both H₂O, CO₂, and CO ice.

In order to see the effect of CO ablation, we performed one simulation that includes it (simulation #8) and one simulation without (simulation #9). We note, however, that the only planetesimals which will experience CO ablation in simulation #8 are those that form at the outermost planetary gap edge. Furthermore, we perform one simulation where we only consider the ablation of water ice (simulation #10). The results from these simulations are presented in Fig. 11 and 12.

When only the ablation of water ice is considered, all mass loss occurs interior to the innermost planet. The total amount of mass that has been ablated after 1 Myr is only about 5%. When CO₂ ablation is added, this value increases to 25%, and there is significant mass loss also in the region around the gap edge of the innermost planet. Compared to the nominal simulation, the initial mass loss rate is much higher, which is because of the larger surface density of the disk. Finally, when CO ablation is added, the total amount of ablated mass increases by another few percent. Another key observation is that at this point, a small amount of mass loss is occurring as far out in the disk as 100 au. This is clearly visible in Fig. 12, where planetesimals forming at the outermost planetary gap edge are losing mass at a slow rate all over the disk.

The effect of decreasing Ṁ₀ is very much the same as in the Solar System simulations: the mass loss rate is lower and the mass loss occurs closer to the central star. In such a disk, less than 5% of the planetesimal mass is lost due to ablation. In summary, planetesimals forming at the edges of planetary gaps in the HL Tau system are not as affected by ablation as those forming at the gap edges of Jupiter and Saturn. This is expected since the planets in the HL Tau system are located much further away from the central star. Nevertheless, if the surface density of the disk is large, the amount of mass loss can still be significant. Most of this mass is released closer to the star than the innermost planet, but some small amounts are also released around the gap edge of the innermost planet. If the surface density of the disk is high enough for CO to be in solid phase where some planetesimals form, these planetesimals can ablate at a slow rate far out in the disk.

6 Pebble flux due to ablation

The ablated material consist of a mixture of silicate grains, carbon grains, and vaporized ices. For simplicity, we assume that all of this material lands on the midplane, which is not a bad approximation given that ablation should be most efficient in the dense midplane of the disk. Depending on where in the disk ablation occurs, the vaporized ices will either re-condense to form solid ice, or remain in the gas phase. For example, if the ablation occurs in between the H₂O and CO₂ iceline, then the vaporized H₂O ice will re-condense to form solid ice, but CO₂ and CO will remain in the gas phase. We assume that re-condensation happens instantaneously. We further assume that all solid material (silicate grains, carbon grains, and re-condensed ices) grow to ~mm-sized pebbles with Stokes number 10⁻³ directly after the ablation occurs. From Paper I, we find that μm-sized grains grow to mm-sizes within ~10 kyr in the inner region of the disk, so given the timescales under consideration, it is a valid approximation.

The solid material from ablation gives rise to a flux of pebbles, which can be calculated given the initial planetesimal formation rate at the gap edges, the rate at which these planetesimals are ablating (right panel of Fig. 6 and 11), and the distribution of the ablated material (left panel of Fig. 6 and 11). We assume that all planetesimals ablate at the same rate and that their ablated material has the same distribution, which gives the correct behavior for the population as a whole. The planetesimal formation rate is taken directly from Paper I in the case of HL Tau and calculated using some assumptions in the case of the Solar System (description in Sect. 6.1). This information is combined togive the total mass ablation rate (gas + solids) in each semimajor axis bin and at each timestep of the simulation Ṁ_abl(r, t). We then remove the mass fraction from ices that do not re-condense (using abundances from Öberg et al. (2011) and assuming a chemical composition of Mg₂ SiO₄ for the silicate grains, and C for the carbon grains), thus, we are left with only the solid (pebble) component.

The flux of pebbles in an annulus of with Δr is: $$F(r,t)=\frac{{\dot{M}}_{\text{abl}}(r,t)}{2\pi r\Delta r},$$(23)

where F is in units of kg s⁻¹ m⁻². The total mass flux F_tot(r, t) as a function of semimajor axis and time is then obtained by numerically integrating 2πrF(r, t) from some r_in to some r_out. Here, we make another assumption that the planetary gaps act as hard barriers, which completely block the flow of pebbles past the gaps. In other words: when calculating the total mass flux interior of the innermost planet, we take the integral from the inner disk edge to the semimajor axis of the innermost planet; when calculating the total mass flux in a region between two planets, we take the integral from the inner planet to the outer planet; and when calculating the mass flux exterior to the outermost planet, we take the integral from the outermost planet to the edge of our ablation array (20 au for SS, 100 au for HL Tau). Finally, we calculate the corresponding pebble-to-gas surface density ratio using Eqs. (8)–(13) from Johansen et al. (2019).

Fig. 11 Left plot: distribution of ablated mass across the disk after 1 Myr for all HL Tau simulations. The values on the y-axis represent the amount of mass that has been ablated in a 0.25 au semimajor axis bin. The dotted lines mark the semimajor axes of the two innermost planets. In the nominal simulation, all mass lossoccurs just interior of the innermost planet. When the surface density is increased (#8), there is also some mass loss occurring further out in the disk. Right plot: total amount of mass that has been ablated as a function of time, for the same data as in the left plot. The colored lines show the average over the 10 simulations and the colored region shows the one standard deviation away from this value. The dotted black lines show the best fit to the curve f(x), and the parameters to the fits can be found in Table 3. Comparison of simulation #8–#10 shows that most mass loss occurs due to the presence of CO2.

6.1 Solar System

In order to calculate how much mass is ablated in the Solar System simulations, we first need an estimate for the initial planetesimal formation rate at Jupiter’s and Saturn’s gap edges. For Jupiter, we simply assume that all the pebbles initially located between the semimajor axis of Jupiter and Saturn are turned into planetesimals at Jupiter’s gap edge (this is the same assumption we used for the collision timescale). Using the parameters of the nominal model and a dust-to-gas ratio of 1%, this gives a result of 16 M_⊕. We calculate that it takes55 kyr for all these pebbles to reach Jupiter’s gap edge, using equations for the pebble drift timescale from Johansen et al. (2019). Assuming a constant drift velocity for all pebbles, this yields an initial planetesimal formation rate of 3 × 10⁻⁴ M_⊕ yr⁻¹ at Jupiter’s gap edge. This is the rate during the first 55 kyr of the simulation, after which it is zero. In order to keep things simple, we assume the exact same formation rate at Saturn’s gap edge.

As a first approximation, we neglect any later generations of planetesimals forming at the gap edges from the ablated material. The resulting pebble flux and pebble-to-gas surface density ratio for the nominal model is presented in Fig. 13 as dotted lines. The highest pebble flux is obtained at 50 kyr, after which it decreases with time. The drop in pebble flux at the H₂O iceline is clearly visible. The CO₂ iceline is located in the gap just interior of Jupiter and, therefore, there is no corresponding drop in the pebble flux.

Since we assumed a hard pebble barrier at the planetary gaps, all pebbles flowing into the gap edge will become trapped. If we use the same assumption as for the initial planetesimal formation rate, all these pebbles will collapse into newer generations of planetesimals. We follow the ablation of these planetesimals and even later generations, assuming the same ablation rate and distribution of ablated material as for the initial population. This results in an increased pebble flux and surface density ratio as compared to our first approximation (see solid lines in Fig. 13). In total, about ~ 15 M_⊕ of pebbles are delivered to 1 au over the course of the simulation. Inside of the water ice line, this amount drops to ~ 12 M_⊕. The presented fluxes are obtained in a young and massive disk, and would not be present in older or low-density disks, thus the formation of transition disks, which are generally old, would not be prevented in such a way.

Fig. 12 Mass loss versus semimajor axis evolution for all HL Tau simulations. The red dots are for planetesimals forming at the gap edge of the innermost planet; the blue dots are for planetesimals forming at the gap edge of the middle planet; and the yellow dots are for planetesimals forming at the gap edge of the outermost planet. Ablation is most efficient for planetesimals forming at the gap edge of the innermost planet. Ablation of CO (#8) causes planetesimals to lose mass far out in the disk. The massive planets deliver strong kicks to the planetesimals, which causes many of them to leave the simulation domain.

Fig. 13 Material ablated from the planetesimal surfaces re-condense and grows to mm-sized pebbles in the disk’s midplane. Top: Mass flux of pebbles as a function of time and semimajor axis for the nominal Solar System simulation. The solid lines take into account the formation and ablation of newer generations of planetesimals forming at the gap edges from this flux of pebbles, while the dotted lines do not. Bottom: Corresponding pebble-to-gas surface density ratios.

Fig. 14 Fluxof pebbles (top panel) and pebble-to-gas surface density ratio (bottom panel) as a function of semimajor axis and time for the nominal HL Tau simulation (similar to Fig. 13). Since the large majority of mass ablation occurs interior of the innermost planet, we do not consider the formation of newer generations of planetesimals at the planetary gap edges.

6.2 HL Tau

In simulation #4 of Paper I (which has the same parameters as in the nominal simulation) we have 281 M_⊕ of planetesimals forming at the gap edges within 1 Myr. By using the planetesimal formation rate from this simulation and combining it with the ablation rate and distribution of ablated mass from this paper, we obtain the amount of mass which is ablated in each semimajor axis bin and at each timestep. The resulting pebble flux and pebble-to-gas surface density ratio is shown in Fig. 14. Since the amount of ablated material being deposited beyond the innermost planet is very small, we do not consider the formation of newer generation of planetesimals.

Due to the relatively low ablation rate in the nominal simulation (only 11% of the initial planetesimal mass is ablated during 1 Myr), the pebble-to-gas surface density ratio only reaches about 0.1%. The maximum pebble flux is about half that of the Solar System, 40 M_⊕ Myr⁻¹ and the total integrated mass flowing past 1 au is 25 M_⊕. In the Solar System simulations, most ablation occurs within the first 50 kyr, after which the ablation rate drops quickly. In the HL Tau simulations, ablation happens on a longer timescale, and the rate does not drop as quickly after the peak at 200 kyr. Because of this, the flux of pebbles remains high for a much longer time; however, we emphasize that these results are for a non-evolving disk and in an evolving disk, the pebble flux would decrease with time along with the gas surface density.

6.3 Asteroid dichotomy

In this work, we use a simple model of the Solar System with only two planets, Jupiter and Saturn, where neither planet nor disk properties evolve with time. Due to these simplifications, comparisons to Solar System data and observations remain conceptual. We nevertheless comment on how pebble drift due to ablation could fit into our understanding of the accretion history of the inner Solar System, specifically the isotopic dichotomy among the parent bodies of iron meteorites (irons) and chondritic meteorites (chondrites).

Meteorites collected on Earth can be broadly classified as either carbonaceous (CC) or noncarbonaceous (NC) based on their distinct isotope compositions (Warren 2011; Kruijer et al. 2017). This isotopic dichotomy, which could reflect either an infall of material with two compositions (Jacquet et al. 2019) or a temperature-dependent destruction of presolar grains (Trinquier et al. 2009; Schiller et al. 2018) hints towards the formation of meteorite parent bodies from two distinct and spatially separated reservoirs (Kruijer et al. 2017). Since NC chondrites come from dry bodies, whereas CC chondrites show hydration features, it is generally envisioned that the NC population represents inner Solar System bodies while members of the CC group would have formed further out.

Iron meteorites are fragments of cores from melted and differentiated planetesimals that accreted within 1 Myr after the formation of calcium–aluminum-rich inclusions (CAIs; Kruijer et al. 2017). Chondrites are fragments of non-molten and non-differentiated planetesimals that accreted around 2–4 Myr after CAI formation (Kita & Ushikubo 2012). The parent bodies of NC and CC irons accreted ~ 0.4 Myr and ~ 0.9 Myr after CAI formation, respectively. Since the two populations have distinct compositions, they must have remained separated even after their formation. Similarly, the NC and CC reservoirs must have been separated at the time of formation of the NC chondrite parent bodies and CC chondrite parent bodies at ~ 2 Myr and 3–4 Myr after CAIs, respectively. In the context of pebbles drift, such constraints imply that pebbles bearing a CC isotopic signature either never entered the inner Solar System region where the parent bodies of NC irons and chondrites formed or, if they had, they were not efficiently accreted by the NC population of objects. Below, we discuss these issues in the context of our results.

The formation of Jupiter is likely to have begun early and resulted in gap-opening in less than 1 Myr (Kruijer et al. 2017). At this point, planetesimals should have started to form at its gap edge. As shown in our simulations, a large fraction of these planetesimals would be scattered into the inner Solar System. While the disk is young and massive, a large fraction of these planetesimals become ablated, resulting in a flux of pebbles with (possibly) CC composition both interior and exterior of Jupiter’s gap. Given the expected decrease in disk density with time, this process should have stopped being relevant at the time of formation of the NC chondrite parent bodies at ~ 2 Myr after CAIs. However, since the NC iron parent bodies formed as early as ~ 0.4 Myr after CAIs, they should have been in place at the time Jupiter reached the pebble isolation mass and started to scatter planetesimals inside of its orbit. Thus, there must have been something preventing the NC iron parent bodies from accreting the pebbles produced by the early ablation of planetesimals originating from beyond Jupiter’s orbit. This process was likely the early formation of planetary embryos interior of Jupiter’s orbit, which caused the inclinations and eccentricities of the planetesimals to be excited, thereby, disconnecting them from the pebbles drifting through the midplane (Schiller et al. 2018). As long as the planetesimals remain excited, they would be unable to efficiently accrete pebbles (e.g., Levison et al. 2015; Liu et al. 2019; Johansen et al. 2015) and should thus have preserved their birth compositions. In general, the described pebble flux was present at a time in the Solar System so early that there is likely no record of it.

7 Planet–planetesimal collisions

In the Solar System simulations, about 1–2% of all planetesimals suffer a collision with a planet. The majority of these collisions are between Jupiter and planetesimals formed at Jupiter’s gap edge. Only about 20% of the collisions are between a planet and a planetesimal not formed at that planet’s gap edge. Following the assumption made in Sect. 6.1, and assuming that the collisions occur before any mass loss has taken place, the total planetesimal mass colliding with Jupiter is about half an Earth mass. Depending on how high up in the atmosphere this material is deposited, this could be relevant for the composition of Jupiter’s atmosphere (e.g., Shibata & Ikoma 2019).

Collisions between planets and planetesimals occur much more frequently in the HL Tau system, where about 15–20% of all planetesimals suffer such a collision. Due to the larger separation between the planetary orbits, nearly all collisions take place between a planet and a planetesimal formed at that planet’s gap edge. If a large fraction of the colliding planetesimal mass is accreted high up in the planet’s atmospheres, this would certainly be relevant for the atmosphere’s compositions.

The collision algorithm employed in this work is rather simple and by making it more advanced, the frequency of collisions could certainly be altered. We use a direct collision search whereby the radius of the planet is determined by assuming a constant density, however, in reality, planets that are accreting gas have an inflated envelope with a significantly larger radius. Thus, the capture radius used in this paper is significantly smaller than it would be in real life, resulting in a smaller number of collisions. The initial formation location of the planetesimals relative to the planet (discussed in Appendix E) could also have an effect on the collision frequency. Another process which has been ignored in this paper but would likely lead to more collisions is planetary migration (Pirani et al. 2019; Carter & Stewart 2020).

8 Conclusions and future studies

In this work, we study the evolution of planetesimals formed at planetary gap edges. To this end, we performed a suite of N-body simulations and considered two planetary systems: (1) the Solar System with Jupiter and Saturn and (2) a system inspired by HL Tau with three planets. We then followed the evolution of planetesimals initiated at the gap edges, where the planetesimals are further subjected to gas drag and mass loss via ablation. We assumed that the mass which is ablated from the planetesimals will re-condense to form pebbles in the disk midplane and we calculated the corresponding pebble flux. We reached the following conclusions regarding the main questions posed in the introduction section:

• 1.

  With regard to the extent to which gravitational interactions with the forming planets redistribute the planetesimals formed at their gap edge, we find that the close proximity to the gap-opening planets results in large orbital excitations, causing the planetesimals to leave their birth location soon after formation. Within ten orbital periods, the semimajor axes of the planetesimals have diffused over several au and after a few hundred orbital periods, the planetesimals are spread out across the entire disk. If the gap-opening planets are massive, as in the HL Tau system, ~40% of the planetesimals become ejected from the system within a million years. These planetesimals could potentially be circularized by galactic tides and end up in the Oort cloud (Brasser & Morbidelli 2013). The ejection efficiency generally increases with time as the disk density decreases and the planetary masses increases. If the planetary gaps were to be much wider than assumed in this paper, resulting in formation locations that are further from the planets, the scattering phase would be slower and delayed.

• 2.

  Regarding the possibility that frictional heating of planetesimals on eccentric orbits may drive the production of a significant amount of dust through surface ablation, we find that planetesimals on eccentric orbits within ~10 au experience efficient ablation and because of this there are no planetesimals entering the innermost disk region. In the nominal Solar System simulation with Ṁ₀ = 10⁻⁷ M_⊙ yr⁻¹, α = 10⁻² and Jupiter and Saturn at 30% of their current masses, 70% of the initial planetesimal mass has been ablated after 500 kyr. Planetesimals formed at Jupiter’s gap edge lose about two times more mass than planetesimals formed at Saturn’s gap edge. The ablation rate is significantly lower in the HL Tau system due to the larger planetary orbits, and only 11% of the initial planetesimal mass has been ablated after 1 Myr for the same disk parameters. Planetesimals that contain CO₂ ice ablate much more efficiently than planetesimals that contain only H₂O ice, and planetesimals that additionally contain CO ice lose a small amount of mass in the furthest regions of the disk. Given the high dependency on the disk density, ablation is only expected to be important in relatively young and massive disks.

• 3.

  Furthermore, we consider how common the occurrence is of collisions between pairs of planetesimals and between planetesimals and planets. We estimated the timescale for planetesimal–planetesimal collisions in our simulations and we found that it is at least an order of magnitude longer than the actual simulated time. In other words, each planetesimal suffers a low risk of colliding with another planetesimal. About 1–2% of all planetesimals in the Solar System simulations will collide with either Jupiter or Saturn, while for the HL Tau system, this value is about ten times higher. Depending on where the material from these planetesimals is deposited in the atmosphere, this could have an impact on the atmospheric compositions.

The material that is ablated from the planetesimal surfaces consists of a mixture of silicate grains, carbon grains, and vaporized ices. A large fraction of these vaporized ices re-condense to form solid ice. The solid material quickly grows to millimeter-sized pebbles in the disk midplane, drift towards the star, and gives rise to a flux of pebbles. In the Solar System, there is a flux of pebbles produced by planetesimal ablation both within and outside of Jupiter’s orbit. Beyond Jupiter, this pebble flux could give rise to newer generations of planetesimals. The total integrated mass that reaches 1 au is ~ 15 M_⊕ in the nominal Solar System simulation and ~25 M_⊕ in the nominal HL Tau simulation. The pebble flux is expected to drop with time as the density in the disk decreases and ablation becomes less important.

In this work (as well as in our previous paper Eriksson et al. 2020), we focus on studying the formation and fate of planetesimals formed at stationary planetary gap edges; however, the reality is that planets migrate. For the formation of planetesimals, this is expected to lead to planetesimals forming in a larger region of the disk (Shibaike & Alibert 2020). The size of this region depends on how far the planet migrates between the onset of planetesimal formation and the halt of migration when the gap becomes very deep. If planetesimals are forming while the migration timescale is smaller than or comparable to the dynamical timescale, it could result in larger excitations of the planetesimal orbits and more planetesimal–planetesimal collisions (Carter & Stewart 2020). The effects of planetary migration will be investigated in a follow-up study.

Acknowledgements

L.E. and A.J. are supported by the Swedish Research Council (Project Grant 2018-04867). T.R. and A.J. are supported by the Knut and Alice Wallenberg Foundation (Wallenberg Academy Fellow Grant 2017.0287). A.J. further thanks the European Research Council (ERC Consolidator Grant 724 687-PLANETESYS), the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine, and the Wallenberg Foundation (Wallenberg Scholar KAW 2019.0442) for research support. The computations were performed on resources provided by the Swedish Infrastructure for Computing (SNIC) at the LUNARC-Centre in Lund, and are partially funded by the Royal Physiographic Society of Lund through grants. Simulations in this paper made use of the REBOUND code which is freely available at http://github.com/hannorein/rebound.

Appendix A No ablation simulation

In Fig. A.1, we show a comparison of simulations with and without surface ablation. When surface ablation is neglected, there is no cooling due to latent heat of vaporization (last term in Eq. (18)). In regions where surface ablation would still have been inefficient, that is, in regions where the surface temperatures are low, this cooling term is very small and thus there is no significant difference between the surface temperatures in the two simulations. This is the case for planetesimals located in Saturn’s scattered disk and beyond. For planetesimals located in Jupiter’s scattered disk and closer to the star, the surfaces can reach several hundred degrees higher temperatures than they do in simulations including ablation. If an ablating planetesimal in this region were to suddenly stop ablating for some reason, its surface temperature would rapidly increase.

Fig. A.1 Comparison of nominal Solar System simulations without surface ablation (top panels), and with surface ablation (bottom panels). The temperature range in the plots is set to 0–300 K to allow for easy comparison; the hottest planetesimals in the top panels have surface temperatures of 443 K, 362 K, and 281 K, respectively, from left to right. The high temperatures obtained in the non-ablation simulations are due to the lack of cooling due to the release of latent heat of vaporization. The population of planetesimals with low eccentricities and short periods in the non-ablation simulations have been circularized by gas drag, these populations do not exist in the ablation simulations since they are removed by ablation before they become circularized.

In runs without ablation, the effect of eccentricity damping by gas drag in the innermost disk becomes visible. Over time, the eccentricities of planetesimals close to the star begin to decrease and after 500 kyr, a population of planetesimals with small and close-to circular orbits emerges. This population does not appear in the ablation simulations since planetesimals become completely ablated before they are circularized.

The disk parameters of the nominal simulation correspond to a young and relatively massive disk. The surface density of disks decreases with time, meaning that both ablation and circularization by gas drag will become less efficient. The difference between the two simulations will thus be less prominent in older disks and, just like ablation, it will occur closer to the star. This means that as the disk evolves with time and ablation becomes inefficient, planetesimals should eventually survive in the asteroid belt region. Finally, we would like to highlight that ablation is a process which occurs in the presence of a major gas component and it could not occur for planetesimals with similar surface temperatures in a system without gas (e.g., comets).

Appendix B Comments on crust formation

In this work, we assume that solid grains embedded in the ice are released along with the vaporized ices. In a scenario where the outflow of vapor is too small to carry along the solid grains, we might instead expect that a dry dust mantle would build up. Once this mantle has become thick, it would severely limit further mass loss. Here, we present a few arguments for why this should not be the case for the efficiently ablating planetesimals in our simulations.

In Appendix A, we show that if an efficiently ablating planetesimal suddenly stops ablating, then the lack of cooling due to latent heat of vaporization causes the surface temperatures to increase by up to several hundred degrees. In other words, if a dust mantle starts to build up and begins to limit the amount of mass loss via ablation, the surface temperaturequickly increases. A similar effect has been studied in the case of comets, although the heating source is different (Orosei et al. 1995; Coradini et al. 1997; Prialnik 1997). The main differences are that: (1) planetesimals can reach much higher surface temperatures than comets if the gas densities are high, which goes against the sustainability of a dust mantle; and (2) the strong headwind around the planetesimals should blow of the dust that is lifted up (whereas dust can fall back onto comets). The formation of dust mantles in the case of planetesimals is thus more difficult than in the case of comets, and even for comets it is not a process which always occurs (Prialnik 1997).

Fig. B.1 Initial vertical velocity for mm-sized pebbles in the outflowing CO2 gas caused by surface ablation. The calculation is done for planetesimals around perihelion in the nominal Solar System simulation. Pebbles are released along with the gas in the case of high eccentric planetesimals interior of Jupiter’s orbit. Lighter gas molecules, smaller pebble sizes, and increased surface temperatures due to a lack of cooling from release of latent heat of vaporization would result in more pebbles being released in this plot.

Fig. B.2 Comparison between the saturated vapor pressure and the gravitational pressure at 0.1 m and 1 m below the surface of a 100 km-sized planetesimal at different temperatures and for different molecules. The temperature down to 1 m below the surface is assumed to be the same as on the surface. The temperature needs to be about twice as high when considering H2O molecules than when considering CO2 molecules inorder for the vapor pressure to win over the gravitational pressure.

The process described above relies on the presence of a significant gas component, which is the case for planetesimals on short orbits in young disks. Further away from the star, where the gas density is lower and ablation less efficient, this process is much less dramatic and a dust mantle might be allowed to form; however, since planetesimals in this region do not ablate very efficiently to begin with, the effect should not modify our results significantly. The production of a dust mantle will also get easier with time, as the gas density decreases and ablation becomes less and less efficient.

In order to further validate our method we have performed a calculation to check whether the escaping gas from planetesimals is strong enough to carry mm-sized pebbles along with it. This is assumed to be the case when the drag force exerted on a pebble by the outflowing gas is stronger than the gravitational pull of the planetesimal (we ignore the centrifugal correction on the gravity of the planetesimal). Because of the strong headwind in the disk, any pebble that is lifted up from the surface will be immediatelycarried away by the headwind. As previously mentioned, this is in contrast with the activity of airless bodies such as comets for which dust grains should fall back if their velocity is smaller than the escape velocity of the body.

The initial vertical velocity of a pebble at the surface is calculated as: $$v_0=v_{\text{th}}-\frac{G{M}_{\text{pl}}{t}_{\text{s}}}{R_{\text{pl}}^2}.$$(B.1)

Here, v_th is the outflow speed of the gas molecules and we perform the calculation for CO₂ molecules. The density of the outflowing gas is calculated as: $$\rho =\frac{{\dot{m}}_{\text{abl}}}{A_{\text{pl}}{v}_{\text{th}}},$$(B.2)

where A_pl is the surface area of the planetesimal. The relative velocity between the pebble and the gas, which is required to calculate the stopping time, is then simply $$v_{\text{rel}}=v_0-v_{\text{th}}=\frac{G{M}_{\text{pl}}{t}_{\text{s}}}{R_{\text{pl}}^2}.$$(B.3)

Since t_s depends also on v_rel this equation has to be solved using a bisection method.

The results from this calculation is presented in Fig. B.1 for planetesimals close to perihelion in the nominal Solar System simulation. In this plot, a pebble is being ejected from the planetesimal when its initial vertical velocity is larger than 0. The results show that pebbles are lifted up by the escaping gas in the case of the most efficiently ablating planetesimals. As explained earlier, ablating planetesimals that start to build up a dust mantle would obtain significantly increased surface temperatures. The corresponding initial vertical velocity of pebbles on those planetesimals would then increase, leading to many more pebbles released than what is suggested by the plot. The above calculation is performed for CO₂ molecules. Using H₂O or CO molecules instead results in higher pebble velocities, provided the ablation rate is the same. The pebble size also affects the calculation, with larger pebbles being more difficult to achieve lift off the surface (e.g., Orosei et al. 1995).

Finally, we ask whether a potential crust would have survived the pressure from the vapor sublimated underneath its surface. To answer this, we performed a simple calculation where we compare the saturated vapor pressure with the gravitational pressure at 0.1 m and 1 m below the surface of a 100 km-sized planetesimal, assuming constant density, hydrostatic equilibrium and that the temperature down to 1 m below the surface is the same as on the surface. The result of this calculation is presented in Fig. B.2. When considering CO₂ vapor, a temperature of about 120–130 K is required to blow off the crust, for H₂O a temperature of 220–240 K is required. Comparison with the top panel of Fig. A.1 shows that these temperatures are obtained for planetesimals sitting inside and interior of Jupiter’s scattered disk in the case of CO₂ vapor, and for highly eccentric planetesimals interior of Jupiter’s orbit in the case of H₂O vapor. Any crust that forms in these regions of the disk is thus likely to be blown off. If a crust were to form further out in the disk on the other hand, it would likely remain and limit any further ablation from the planetesimal surface. However, we note that the formed crust is likely to be porous with a density smaller than what is assumed for this calculation, so the temperatures required to blow off the crust would be lower than what is suggested in Fig. B.2.

Appendix C Saturated vapor pressure

Figure C.1 shows the saturated vapor pressure as a function of the planetesimal surface temperature. The saturated vapor pressure is calculated using polynomial expressions from Fray & Schmitt (2009). These expressions are valid in the following temperature range: 0–273.16 K, 40–216.58 K, and 14–68.1 K for H₂O, CO₂, and CO ice,respectively. Looking at Fig. C.1, there is a smooth extrapolation to temperatures above the stated temperature range for all ices. This is not the case below the temperature range. For CO₂ and CO ice the polynomial turns upwards at low temperatures; therefore, we introduce floor-values at 40 K and 14 K, respectively. For H₂O the polynomial heads to minus infinity as the temperature decreases; therefore, we introduce a floor-value at 60 K to prevent numerical issues. Additional tests show that lowering this floor-value does not affect the results. Since there are no issues at the upper temperature range, we do not introduce any similar constraints there.

Fig. C.1 Saturated vapor pressure as a function of planetesimal surface temperature for H2O, CO2 and CO ice. The dotted lines mark our floor values, which are used for numerical reasons.

Appendix D Mass ablation rate

Figure D.1 shows the mass ablation rate as a function of the semimajor axis for H₂O, CO₂, and CO ice (calculated using Eq. (17)). The ablation of H₂O ice is only efficient in the very innermost part of the disk and does not generate very high ablation rates. The ablation rates for CO ice are on the contrary high enough to completely sublimate a planetesimal on a ~ 10 au circular orbit within a year. Planetesimals formed beyond the CO iceline should thus lose mass much further out in the disk than planetesimals formed interior of the CO iceline. The ablation rates for CO₂ ice lie in between the ablation rates for H₂O and CO ice.

Fig. D.1 Mass ablation rate as a function of semimajor axis for H2O, CO2, and CO ice.The planetesimal surface temperature is set to equal integer multiples of the disk temperature. The mass of one planetesimal is marked as a purple line on the plot; above this line, a planetesimal becomes completely ablated within one year.

Appendix E Varying the gap width and the planetesimal formation location

The planetesimals in our simulations are initiated between 1–2 Hill radii away from the planets, which is roughly the region in which planetesimals form in Paper I. The formation location in Paper I is strongly affected by (1) the width and depth of the planetary gap; and (2) our assumption that the planetesimals form at the location where the SI is triggered. Regarding (1): There are numerous gap prescriptions in the literature and they often disagree on the resultant gap shape. Wider gaps result in planetesimal formation further from the planet and deeper gaps result in planetesimals forming in a narrower region. Furthermore, 1D simulations generally produce deeper and narrower gaps than their 2D or 3D analogs. Regarding (2), Carrera et al. (2021) performed 3D simulations of the streaming instability in the presence of a pressure bump and found that the planetesimals form inward of the pressure bump that initiated their growth, so not at the exact location where the SI is triggered. Taken together, this means that we do not know exactly where the planetesimals form relative to the planet location. The formation location might be both closer or further away from the planet than it is in our simulations.

In order to study how different formation locations affect our results, we performed three nominal Solar System simulations where we placed the planetesimals at either 1–3, 3–5, or 5–7 Hill radii away from the planets and extended the simulation time to 1 Myr. The gas surface density profile was kept unchanged. The results of these simulations are presented in Figs. E.1 and E.2. Generally, the excitation of the orbital eccentricities is slower when the planetesimals are formed further away from the planets. Because of this, some planetesimals are still located close to their birth locations at the end of the simulation. The slower orbital excitation also results in lower ablation rates, especially in the beginning of the simulation; however, the ablation continues for a much longer period of time. This would result in the flux of pebbles through the midplane being lower, but persisting for a longer time. The presented results indicate that ablation is important also for planetesimals forming further away from the planets, but it happens on a longer timescale.

Fig. E.1 Eccentricity, semimajor axis, and surface temperature evolution for three nominal Solar System simulations with 500 planetesimals per simulation and varying planetesimal formation locations. In the top panel, the planetesimals are initiated between 1 and 3 Hill radii away from the planets, in the middle panel between 3 and 5 Hill radii, and in the bottom panel between 5 and 7 Hill radii. The solid black lines mark the perihelia of the planets. Eccentricity excitation occurs much earlier and at a much faster speed for planetesimals initiated close to the planets. When the planetesimals are initiated between 5 and 7 Hill radii away from the planets, a population of planetesimals still remain around their birth location at the end of the simulation.

Fig. E.2 Distribution of ablated mass as a function of semimajor axis (left) and time (right) for the same three simulations as in Fig. E.1. Left: Distribution of ablated mass across the disk after 1 Myr. The values on the y-axis represent the amount of mass that has been ablated in a 0.2 au semimajor axis bin. The dotted lines mark the semimajor axes of the planets. The distribution of ablated mass is very similar in all three simulations, with the majority of mass being lost either just within or just outside of Jupiter’s orbit. Right: Total amount of mass that has been ablated as a function of time for the same data as in the left plot. The colored lines show the average over the five simulations and the colored region shows the one standard deviation away from this value. The ablation rate decreases as the planetesimal formation location with regard to the planet increases; however, the ablation continues for a longer period of time.

In the above simulations, we only varied the formation location of the planetesimals, not the shape of the planetary gas-gap. Narrower gaps will result in mass loss closer to the planet location, and shallower gaps will result in that more mass is lost around the planet location. We confirmed this by performing a simulation with half the gap-width and a much shallower gap. We found that the total amount of mass loss did not change, only the distribution of the ablated material around the location of the planets did, in fact, change.

Appendix F Dynamical evolution of six planetesimals

Figure F.1 shows the dynamical evolution for the same 6 planetesimals as in Fig. 3. The “red” planetesimal forms at the gap edge of Jupiter, and is a member of Jupiter’s scattered disk until it obtains a strong kick and becomes a member of Saturn’s scattered disk, where it remains until the end of the simulation. The “yellow” planetesimal is a member of Jupiter’s scattered disk for a few thousand years, until a kick by Jupiter places it on an orbit interior of the scattered disk, where it quickly becomes ablated. The “green” planetesimal is a member of Jupiter’s scattered disk during most of its lifetime until it is kicked towards the sun and becomes ablated.

Fig. F.1 Semimajor axis and eccentricity evolution for the same data as in Fig. 3, with the same color-coding. Planetesimals formed at Jupiter’s gap edge are plotted in the top panels, and planetesimals formed at Saturn’s gap edge are plotted in the bottom panels. Left: evolution of the planetesimals semimajor axis (thick colored lines), perihelion, and aphelion (thin colored lines). The semimajor axes, perihelia and aphelia of Jupiter and Saturn (thick and thin black lines) are included as well. Right: eccentricity and semimajor axis tracks, with dots to mark 100 kyr of evolution.The black lines mark the perihelia of Jupiter and Saturn.

The “light-blue” planetesimal forms at the gap edge of Saturn, and sits on an orbit just within Saturn’s scattered disk for about 150 kyr, when it is then scattered towards Jupiter’s scattered disk and obtains a large eccentricity, which increases with time until the planetesimal leaves the simulation domain. The “dark-blue” planetesimal is a member of Saturn’s scattered disk for 50 kyr, until a strong encounter with Saturn places it on an orbit interior of Jupiter, where it becomes circularized. The “purple” planetesimal is a member of Saturn’s scattered disk for most of its lifetime, except during a brief period where it has an orbit that is slightly outside the scattered disk.

References

All Tables

All Figures

Fig. 1 Eccentricity, semimajor axis, and surface temperature evolution for the 1000 planetesimals in the nominal Solar System simulation. The presented orbital parameters have been averaged over 100 yr and the surface temperatures are the maximum values during the same time period, which roughly corresponds to the surface temperatures at perihelion. The temperature of the surrounding gas (Tdisk) is shown as a colorbar in the left panel. The solid black lines mark the orbits with a perihelion corresponding to Jupiter’s and Saturn’s location. The planetesimals which cluster around these lines are (at least momentarily) members of Jupiter’s or Saturn’s scattered disk. Planetesimals that are scattered interior of Jupiter’s orbit obtain high surface temperatures at perihelion and experience efficient ablation.

In the text

Fig. 2 Semimajor axis and surface temperature evolution for the 1000 planetesimals in the nominal Solar System simulation (same data as in Fig. 1). For simulation times longer than 1 kyr, we average the semimajor axes over 100 yr and show the maximum surface temperatures during the same time period; whereas for shorter simulations times, we show non-averaged values. Upper panel: planetesimals formed at the gap edge of Jupiter and lower panel: planetesimals formed at the gap edge of Saturn. The close proximity to the giant planets result in continuous scatterings, which causes the fast semimajor axis diffusion that is displayed in the plot.

In the text

Fig. 3 Mass loss versus semimajor axis for 6 selected planetesimals from the nominal Solar System simulation. The mass loss is calculated as 1 − M(t)∕M(t = 0). Filled circlesmark 100 kyr of evolution and the formation location of the planetesimals (shown at the bottom of the plot) has been marked “ST”. The time before the planetesimals first experience mass loss, that is, the time before they appear on the plot, has been included in the legend. The first three legend entries are for planetesimals formed at the gap edge of Jupiter and the following three are for planetesimals formed at the gap edge of Saturn. The dotted lines mark the semimajor axis of Jupiter and Saturn. Planetesimals which are scattered interior of Jupiter’s orbit lose mass at a high rate, while those which are scattered exterior of Saturn’s orbit experience little mass loss.

In the text

Fig. 4 Mass loss versus semimajor axis evolution for the 1000 planetesimals in the nominal Solar System simulation (same data as in Figs. 1 and 2). The thin grey lines mark the perihelion and aphelion of the planetesimal orbits. The number of planetesimals formed at Jupiter’s (red dots) respective Saturn’s (blue dots) gap edge which do not appear on the plot, because they have been either: completely ablated; collided with a planet; or ejected beyond the simulations domain, is written in each panel. Planetesimals formed at the gap edge of Jupiter generally experience more ablation than planetesimals formed at the gap edge of Saturn. About 50% of all planetesimals have become completely ablated after 500 kyr.

In the text

Fig. 5 Time evolution of the timescale for planetesimal-planetesimal collisions for the nominal Solar System simulation. The collision timescale was calculated using Eq. (20) and assuming an initial planetesimal mass of 16 M⊕ per gap edge. The dotted line marks where the collision timescale equals the time of the simulation. Since the resulting collision timescale for one planetesimal is several orders of magnitude larger than the actual simulated time, planetesimal–planetesimal collisions can be safely ignored in our simulations and will not constitute a significant path to replenishing the dust component in the disk.

In the text

Fig. 6 Distribution of ablated mass as a function of semimajor axis (left) and time (right) for all Solar System simulations. Left plot: Distribution of ablated mass across the disk after 500 kyr. The values on the y-axis represent the amount of mass that has been ablated in a 0.2 au semimajor axis bin. The dotted lines mark the semimajor axes of Jupiter and Saturn. Most mass loss occur in the region just interior to and exterior to Jupiter’s orbit. Right plot: Total amount of mass that has been ablated as a function of time, for the same data as in the left plot. The colored lines show the average over the ten simulations, and the colored region shows the one standard deviation away from this value. The dotted black lines show the best fit to the curve f(x); and the parameters to the fits can be found in Table 2. The mass ablation rate is highly dependent on the gas surface density, where a high surface density results in a high degree of ablation (#5) and a low surface density results in low degree of ablation (#6).

In the text

Fig. 7 Mass loss versus semimajor axis evolution for all Solar System simulations. The masses of the planets affect how far the planetesimals are scattered, with more planetesimals leaving the simulation domain when the planetary mass is increased (#4). A large surface density (#5) generally results in efficient ablation and few scatterings beyond the simulation domain since the planetesimals are ablated before they become scattered. The opposite is the case for a low surface density (#6) and ablation also occurs closer to the star.

In the text

Fig. 8 Eccentricity, semimajor axis, and surface temperature evolution for the 1500 planetesimals in the nominal HL Tau simulation, produced in a similar manner as Fig. 1. The solid black lines mark the perihelia of the planets. Many planetesimals end up in the scattered disks of the planets, where they obtain high eccentricities. Only planetesimals scattered tointerior of the innermost planet’s orbit obtain perihelion temperatures large enough for ablation to take place.

In the text

Fig. 9 Semimajor axis and surface temperature evolution for the 1500 planetesimals in the nominal HL Tau simulation, produced in a similar manner as in Fig. 2 (same data as in Fig. 8). Planetesimals formed at the gap edge of the innermost planet is shown in the upper panel. For the middle and outermost planet, these are shown in the middle and bottom panels, respectively. Since the planetary masses are set to equal the pebble isolation mass, and the pebble isolation mass increases with semimajor axis, the planetesimals formed at the gap edge of the outermost planet experience stronger scatterings and diffuse faster than those formed at the innermost planet’s gap edge.

In the text

Fig. 10 Time evolution of the timescale for planetesimal–planetesimal collisions for the nominal HL Tau simulation. The collision timescale was calculated using Eq. (20) and using an initial planetesimal mass of 281 M⊕, following the results from simulation #4 of Paper I. The dotted line marks where the collision timescale equals the time of the simulation. Even though the collision timescale for one planetesimal can be as low as 104 yr, it is still roughly an order of magnitude larger than the actual simulated time, meaning that planetesimal-planetesimal collisions can be safely ignored in our simulations as a source of dust production.

In the text

Fig. 11 Left plot: distribution of ablated mass across the disk after 1 Myr for all HL Tau simulations. The values on the y-axis represent the amount of mass that has been ablated in a 0.25 au semimajor axis bin. The dotted lines mark the semimajor axes of the two innermost planets. In the nominal simulation, all mass lossoccurs just interior of the innermost planet. When the surface density is increased (#8), there is also some mass loss occurring further out in the disk. Right plot: total amount of mass that has been ablated as a function of time, for the same data as in the left plot. The colored lines show the average over the 10 simulations and the colored region shows the one standard deviation away from this value. The dotted black lines show the best fit to the curve f(x), and the parameters to the fits can be found in Table 3. Comparison of simulation #8–#10 shows that most mass loss occurs due to the presence of CO2.

In the text

Fig. 12 Mass loss versus semimajor axis evolution for all HL Tau simulations. The red dots are for planetesimals forming at the gap edge of the innermost planet; the blue dots are for planetesimals forming at the gap edge of the middle planet; and the yellow dots are for planetesimals forming at the gap edge of the outermost planet. Ablation is most efficient for planetesimals forming at the gap edge of the innermost planet. Ablation of CO (#8) causes planetesimals to lose mass far out in the disk. The massive planets deliver strong kicks to the planetesimals, which causes many of them to leave the simulation domain.

In the text

Fig. 13 Material ablated from the planetesimal surfaces re-condense and grows to mm-sized pebbles in the disk’s midplane. Top: Mass flux of pebbles as a function of time and semimajor axis for the nominal Solar System simulation. The solid lines take into account the formation and ablation of newer generations of planetesimals forming at the gap edges from this flux of pebbles, while the dotted lines do not. Bottom: Corresponding pebble-to-gas surface density ratios.

In the text

	Fig. 14 Fluxof pebbles (top panel) and pebble-to-gas surface density ratio (bottom panel) as a function of semimajor axis and time for the nominal HL Tau simulation (similar to Fig. 13). Since the large majority of mass ablation occurs interior of the innermost planet, we do not consider the formation of newer generations of planetesimals at the planetary gap edges.
In the text

Fig. A.1 Comparison of nominal Solar System simulations without surface ablation (top panels), and with surface ablation (bottom panels). The temperature range in the plots is set to 0–300 K to allow for easy comparison; the hottest planetesimals in the top panels have surface temperatures of 443 K, 362 K, and 281 K, respectively, from left to right. The high temperatures obtained in the non-ablation simulations are due to the lack of cooling due to the release of latent heat of vaporization. The population of planetesimals with low eccentricities and short periods in the non-ablation simulations have been circularized by gas drag, these populations do not exist in the ablation simulations since they are removed by ablation before they become circularized.

In the text

Fig. B.1 Initial vertical velocity for mm-sized pebbles in the outflowing CO2 gas caused by surface ablation. The calculation is done for planetesimals around perihelion in the nominal Solar System simulation. Pebbles are released along with the gas in the case of high eccentric planetesimals interior of Jupiter’s orbit. Lighter gas molecules, smaller pebble sizes, and increased surface temperatures due to a lack of cooling from release of latent heat of vaporization would result in more pebbles being released in this plot.

In the text

Fig. B.2 Comparison between the saturated vapor pressure and the gravitational pressure at 0.1 m and 1 m below the surface of a 100 km-sized planetesimal at different temperatures and for different molecules. The temperature down to 1 m below the surface is assumed to be the same as on the surface. The temperature needs to be about twice as high when considering H2O molecules than when considering CO2 molecules inorder for the vapor pressure to win over the gravitational pressure.

In the text

	Fig. C.1 Saturated vapor pressure as a function of planetesimal surface temperature for H2O, CO2 and CO ice. The dotted lines mark our floor values, which are used for numerical reasons.
In the text

	Fig. D.1 Mass ablation rate as a function of semimajor axis for H2O, CO2, and CO ice.The planetesimal surface temperature is set to equal integer multiples of the disk temperature. The mass of one planetesimal is marked as a purple line on the plot; above this line, a planetesimal becomes completely ablated within one year.
In the text

Fig. E.1 Eccentricity, semimajor axis, and surface temperature evolution for three nominal Solar System simulations with 500 planetesimals per simulation and varying planetesimal formation locations. In the top panel, the planetesimals are initiated between 1 and 3 Hill radii away from the planets, in the middle panel between 3 and 5 Hill radii, and in the bottom panel between 5 and 7 Hill radii. The solid black lines mark the perihelia of the planets. Eccentricity excitation occurs much earlier and at a much faster speed for planetesimals initiated close to the planets. When the planetesimals are initiated between 5 and 7 Hill radii away from the planets, a population of planetesimals still remain around their birth location at the end of the simulation.

In the text

Fig. E.2 Distribution of ablated mass as a function of semimajor axis (left) and time (right) for the same three simulations as in Fig. E.1. Left: Distribution of ablated mass across the disk after 1 Myr. The values on the y-axis represent the amount of mass that has been ablated in a 0.2 au semimajor axis bin. The dotted lines mark the semimajor axes of the planets. The distribution of ablated mass is very similar in all three simulations, with the majority of mass being lost either just within or just outside of Jupiter’s orbit. Right: Total amount of mass that has been ablated as a function of time for the same data as in the left plot. The colored lines show the average over the five simulations and the colored region shows the one standard deviation away from this value. The ablation rate decreases as the planetesimal formation location with regard to the planet increases; however, the ablation continues for a longer period of time.

In the text

Fig. F.1 Semimajor axis and eccentricity evolution for the same data as in Fig. 3, with the same color-coding. Planetesimals formed at Jupiter’s gap edge are plotted in the top panels, and planetesimals formed at Saturn’s gap edge are plotted in the bottom panels. Left: evolution of the planetesimals semimajor axis (thick colored lines), perihelion, and aphelion (thin colored lines). The semimajor axes, perihelia and aphelia of Jupiter and Saturn (thick and thin black lines) are included as well. Right: eccentricity and semimajor axis tracks, with dots to mark 100 kyr of evolution.The black lines mark the perihelia of Jupiter and Saturn.

In the text