© ESO 2019

1 Introduction

The standard model of giant planet formation is the core-accretion mechanism (Mizuno 1980; Stevenson 1982; Pollack et al. 1996; Guilera et al. 2010; Helled et al. 2014). According to this model, the planet forms from the accretion of planetesimals onto a solid core until it has enough mass to start accreting gas from the protoplanetary disk. When the cross-over mass is reached, the mass of the envelope equals the core mass (at ~10 M_⊕), and the planet starts to accrete large quantities of gas in a short period of time; this process is known as runaway gas accretion. At some point, due to mechanisms that are not yet well understood, the accretion of gas onto the planet is limited. Finally, the planet evolves, cooling and contracting at constant mass.

An alternative model has been developed in the past few years, also based on the accretion of solid material, for the formation of giant planets. This model is based on the accretion of small-sized particles (called pebbles) for the formation of the planet core. Unlike planetesimals, pebbles can be accreted by the full Hill sphere making their accretion rates significantly larger than planetesimal accretion rates (Ormel & Klahr 2010; Lambrechts & Johansen 2012).

In the standard core accretion mechanism an important issue is the formation of a massive core before the dissipation of the gaseous component of the protoplanetary disk. In this framework, the process of accretion of planetesimals plays a fundamental role. As the planetary embryo grows, the surrounding planetesimals increase their relative velocities due to the gravitational stirring produced by the embryo, leading to disruptive collisions among them. This effect produces a cascade of collision fragments. Inaba et al. (2003), Kobayashi et al. (2011), Ormel & Kobayashi (2012), and Guilera et al. (2014) found that excessive fragmentation may cause the smaller fragments produced by these collisions to drift inwards through gas drag and a significant amount of mass could be lost stalling the oligarchic growth. Moreover, when the core reaches several Earth masses the relative planetesimal velocities become high enough to produce super-catastrophic collisions among planetesimals (Guilera et al. 2014; Chambers 2014). In this case, Guilera et al. (2014) considered that the planetesimals involved in such a collision are pulverized and the mass of such bodies is lost. On the other hand, Chambers (2014) found that small planetesimals collide frequently and produce rapid embryo growth. Unlike Guilera et al. (2014), Chambers (2014) considered that the mass that would go into fragments smaller than the minimum size considered in the model (referred to as pebbles of second generation) is assigned to become an equivalent mass of the minimum size adopted in their model, since they assume that tiny fragments quickly coagulate into new pebbles.

Despite the different planetesimal fragmentation models, during the process of planetary formation, fragmentation of planetesimals is an important effect that may inhibit or favor the formation of giant planet cores. In a collision between two planetesimals, the catastrophic impact energy threshold required to fragment and disperse 50% of the target mass (also known as specific impact energy) is an important function that has to be defined in models of planetesimal fragmentation. This function depends on many factors of the collision, such as the impact velocity, the planetesimal type of material and its porosity, the sizes of the planetesimals, and other factors that may affect the outcome of the collision (see Jutzi et al. 2015). For example, Benz & Asphaug (1999) found that a target of a fixed size composed of ice material has a catastrophicimpact energy threshold lower than if it were composed of basalt. Benz (2000) also found that low-velocity collisions are more efficient in destroying and dispersing basaltic bodies than collisions at higher velocities. From smoothed particle hydrodynamics (SPH) simulations, Jutzi et al. (2010) studied the porosity of the planetesimals in a collision and found that, in the strength regime, porous targets are more difficult to disrupt than non-porous ones, while in the gravity regime the outcome of the collision for porous targets depends on gravity and porosity, unlike non-porous ones that are also controlled by the strength. Recently, Beitz et al. (2011) and Bukhari Syed et al. (2017) showed in low-velocity impact experiments that the strength of compacted dust aggregates is much weaker than that ofporous rocks. This outcome is analogous to Jutzi et al. (2010). It is important to remark here that most models of planet formation that incorporate the planetesimal fragmentation process use the prescription of Benz &Asphaug (1999) for the catastrophic impact energy threshold for non-porous basalts at impact velocities of 3 km s⁻¹ (Kobayashi et al. 2011; Ormel & Kobayashi 2012; Guilera et al. 2014; Chambers 2014). However, Guilera et al. (2014) showed thatduring the formation of massive cores the relative planetesimal velocities are incremented from centimeters per second and meters per second to kilometers per second as planet grows.

In our previous work (Guilera et al. 2014), we incorporated a model of planetesimal fragmentation into our global model of giant planet formation (Guilera et al. 2010). We found that the formation of massive cores is only possible from an initial population of big planetesimals and massive disks. In other situations, planetesimal fragmentation tends to inhibit the formation of massive cores before the dissipation of the disk. In this new work, we incorporate into our planetesimal fragmentation model the dependence of the catastrophic impact energy threshold with material composition and planetesimal relative velocities. We also incorporate different velocity regimes for planetesimal relative velocities (Keplerian shear and dispersion-dominated regimes) into the planetesimal fragmentation model, and pebble accretion rates for the small particles produced by the collision process into the global model of giant planet formation. The aim of this work is to analyze the impact on the formation of massive cores of the newly incorporated phenomena, especially planetesimals’ relative velocities and planetesimal compositions.

This paper is organized as follows: in Sect. 2, we describe the improvements of our model of giant planet formation and the improvements of the model of planetesimal fragmentation; in Sect. 3, we present the results of our simulations for the formation of a giant planet, and finally, in Sect. 5 we present the summary and conclusions of our results.

2 Our model of giant planet formation

In a series of previous works (Guilera et al. 2010, 2014) we developed a numerical model that describes the formation of giant planets immersed in a protoplanetary disk that evolves in time. In our model, the protoplanetary disk is represented by a gaseous and a solid components. Planets grow by the simultaneous accretion of solids and gas. The solid component of the disk evolves by planet accretion, radial drift due to nebular drag, and collisional evolution, while the gaseous component evolves by an exponential decay (see Appendix A.1 for technical details on the model).

In this work we improve our planetesimal fragmentation model developed in Guilera et al. (2014), incorporating a dependence of $Q_{\text{D}}^{*}$, the catastrophic impact energy threshold, with the planetesimal relative velocities and the planetesimal composition considering a mix of non-porous ices and, as a proxy for the rocky material of planetesimals, non-porous basalts. We also incorporate pebble accretion rates for the small fragments that are products of the planetesimal fragmentation process (hereafter, pebbles of second generation), and different velocity regime models for the calculation of low and high relative planetesimal velocities. The improvements in our model are presented in the next subsections.

2.1 Improvements in the solid accretion rates

In the past few years, a new model of solid particles accretion has been proposed as an alternative to the core accretion model for the formation of solid massive cores. Such a model proposes that massive cores, precursors of giant planet cores, can be quickly formed by 100–1000 km embryos accreting centimeter-sized particles known as pebbles. Ormel & Klahr (2010) and Lambrechts & Johansen (2012) found that pebbles, particles with Stokes numbers below unity (St ≲ 1), are strongly coupled to the gas and can be very efficiently accreted by hundred-thousand kilometer-sized embryos. While planetesimals can be accreted by a fraction of the Hill radius of the growing planet, pebbles can be accreted by the full Hill radius, making pebble accretion rates significantly larger than planetesimal accretion rates.

However, the dominant initial size of the accreted solids (centimeter-sized particles named pebbles or 1–100 km-sized bodies named planetesimals) is unknown (Helled et al. 2014; Johansen et al. 2014). At the beginning of our simulations, the solid component of the protoplanetary disk is composed of planetesimals of a radius of 100 km with the distribution given by the framework of the oligarchic growth (Kokubo & Ida 1998, 2000, 2002; Ida & Makino 1993). Thus, we assume that there are no initial pebbles in our model. However, pebbles appear as a result of the planetesimal collisional evolution, that is, second generation pebbles.

Following Guilera (2016) and Guilera & Sándor (2017), we improve the solid accretion rates by incorporating the pebble accretion rate for pebbles of second generation in addition to the planetesimal accretion. For planetesimals we adopt the accretion rates given by Inaba et al. (2001) and for pebbles of second generation we use the pebble accretion rates given by Lambrechts et al. (2014) with a reduction factor that takes into account the scale height of the pebbles in comparison with the Hill radius of the planet (see Appendix A.2).

2.2 Improvements to our planetesimal fragmentation model

Our planetesimal fragmentation model, which is based on the Boulder code (Morbidelli et al. 2009 and Supplementary Material), describes the collisional evolution of the planetesimal population. A brief description of the planetesimal fragmentation model and technical details are presented in Appendix B.

2.2.1 Velocities and probabilities of collision regimes

Following Morbidelli et al. (2009), the impact rate between targets and projectiles estimated in our model is given by $$\text{Impact rate}{|}_{\text{A}}=\frac{\alpha V_{\text{rel}}}{4Ha(\delta a+2a{e}_{\text{P}})}E\left(1+b\frac{V_{\text{esc}}^2}{V_{\text{rel}}^2}\right){(R_{\text{T}}+R_{\text{P}})}^2,$$(1)

where V_rel is the relative velocity between targets and projectiles, α is a coefficient that depends on V_rel, H is the mutual scale height, a and δa are the mean semi-major axis and the width of the annulus that contains the targets and projectiles, respectively, e_P is the eccentricity of the projectiles, E is a coefficient that considers the deviation of the gravitational focusing corresponding to the two-body problem at low relative velocities, b is also a function of V_rel, and V_esc is the mutual escape velocity. Finally, R_T and R_P are the radius of the targets and projectiles, respectively (see Morbidelli et al. 2009, for details).

It is important to note, that the relative velocity among planetesimals used in the Boulder code is the dispersion velocity (the velocity corresponding to the dispersion regime). However, when the planetesimal relative velocity tends to zero, the impact rate between them is undetermined. Weidenschilling (2011) noticed this problem and incorporated the relative velocity corresponding to the Keplerian shear. Following Weidenschilling (2011), we incorporate in our planetesimal fragmentation model different velocity regimes and probabilities of collision for lower relative velocities given by Greenberg et al. (1991) to calculate collision rates more accurately. We consider three different regimes and their transitions according to Greenberg et al. (1991): regime A: dominance by random motion (impact rate given by Eq. (1)); regime B: dominance by Keplerian shear motion; regime C: Keplerian shear dominance in a very thin disk (see Appendix B.1 for details).

2.2.2 Catastrophic impact energy threshold

The catastrophic impact energy threshold per unit target mass $Q_{\text{D}}^{*}$ is the energy needed to fragment and disperse half of the target mass in an impact, that is, the threshold for catastrophic disruption. This quantity plays an important role in the collisional evolution of the planetesimal population. As we mentioned in Sect. 1, $Q_{\text{D}}^{*}$ depends on many factors of the collision; particularly important are the relative velocity among planetesimals, which determines the impact velocity of the collision, and the planetesimals’ composition.

Figure 1 represents $Q_{\text{D}}^{*}$ as a function of the radius of a non-porous basalt target for three different values of the impact velocity (~ 25 m s⁻¹, 3 and 5 km s⁻¹), while Fig. 2 represents this threshold for a non-porous icy target for two values of the impact velocity (0.5 and 3 km s⁻¹) (see Appendix B.2). For basalt targets, we can see that for a fixed target’s radius, the smaller the impact velocity, the smaller the $Q_{\text{D}}^{*}$. This phenomenon can have an important effect during the formation of massive cores due to the fact that, initially, planetesimal relative velocities are low and are then incremented as the planet grows. We remark again that most works of planetary formation that include planetesimal fragmentation use a fixed value of $Q_{\text{D}}^{*}$ for non-porous basalt targets at an impactvelocity of 3 km s⁻¹ (Ormel & Kobayashi 2012; Chambers 2014). Given the functional dependency of $Q_{\text{D}}^{*}$ with target’s radius for different impact velocities given by Benz & Asphaug (1999) and Benz (2000), for a fixed target’s radius, we implemented an interpolation between the impact velocity dependencies to get an improved value of $Q_{\text{D}}^{*}$ as a function of the impact velocity. We also include an interpolation using the different curves of $Q_{\text{D}}^{*}$ of Benz & Asphaug (1999) for different velocities for ices. We note that if the impact velocities are greater or lower than the velocities corresponding to the upper and lower curves of $Q_{\text{D}}^{*}$ adopted, we do not extrapolate $Q_{\text{D}}^{*}$. In that case, we adopt $Q_{\text{D}}^{*}$ corresponding to the maximum or minimum velocity used in Benz & Asphaug (1999) and Benz (2000).

On the other hand, according to Lodders (2003), beyond the iceline the amount of solid mass in the protoplanetary disk increases by a factor of two meaning that 50% of the material behind the iceline underwent condensation in the primitive solar system. The ice-to-rock ratio derived from trans-Neptunian objects, comets, and irregular satellites of giant planets (McDonnell et al. 1987; Stern et al. 1997; Johnson & Lunine 2005) confirms this. As we are interested in the formation of giant planets behind the iceline, we assume that planetesimals are composed of ices and basalts following Lodders (2003). Finally, to implement the dependence of $Q_{\text{D}}^{*}$ with the impact velocity and the target’s composition, we first interpolate between the curves of $Q_{\text{D}}^{*}$ as a function of the target’s radius for the different impact velocities, obtaining the values of $Q_{\text{D}}^{*}$ for each pure material (basalts and ices). Then, we make a linear combination between them, depending on the percentage of basalts and ices that we define for the solids (in this work we considered three cases, planetesimals composed purely by basalts, composed purely by ices, and composed 50% by basalts and 50% by ices). We remark here again that in all the cases $Q_{\text{D}}^{*}$ is calculated by Eq. (B.8), but using the effective radius (r_eff = 3(M_T + M_P)∕4πρ) instead of the target’s radius.

It should be noted that Leinhardt & Stewart (2012) also obtain a derivation of a general catastrophic disruption law. However, their work is focused on the gravitational regime, and for bodies of different porosities, while in our model we adopt compatible laws valid in the gravitational as well as in the strength regime for non-porous bodies.

Fig. 1 Catastrophic disruption thresholds for basalt targets at impact velocities of 20–30 m s−1 (Benz 2000), 3 and 5 km s−1 (Benz & Asphaug 1999). The blue points correspond to the discrete data extracted from Benz (2000).

Fig. 2 Catastrophic disruption thresholds for icy targets at impact velocities of 0.5 and 3 km s−1 (Benz & Asphaug 1999).

3 Results

We aim to study the impact of the improvements of our planetesimal fragmentation model discussed above on the formation of a giant planet. We carried out a set of different simulations including one phenomenon at a time. Our simulations start at the beginning of the oligarchic growth with a moon-sized embryo located at 5 au from the central star. Initially, the embryo is immersed in an homogeneous single-sized population of planetesimals of 100 km radius and the disk is ten times more massive than the minimum mass solar nebula (MMSN, Hayashi 1981). The simulations stop when the mass of the envelope equals the core mass, that is, when the cross-over mass is achieved (in this case we consider that the planet ends its formation in a short period of time after the gaseous runaway growth starts), or at 6 Myr when we consider that the gas of the disk is dissipated.

3.1 Baseline model

In Guilera et al. (2014) we included a first approach to the modeling of planetesimal fragmentation in our model of giant planet formation (Guilera et al. 2010), and we studied, for different initial planetesimal sizes and disk masses, how the collisional evolution of the planetesimal population modified the planetary formation process. We remark that in Guilera et al. (2014) the accretion rates of particles with Stokes number less than the unity were calculated using the prescription derived by Inaba et al. (2001) in the low-velocity regime. We found that only for initial large planetesimals (r_P = 100 km) and massive disks, and only if most of the mass loss in collisions was distributed in larger fragments (see Eq. (B.3)) did planetesimal fragmentation favor the relatively rapid formation of a massive core (larger than 10 M_⊕). If smaller planetesimals are considered, the planetesimal fragmentation process inhibits the formation of massive cores. We remark that these results are in agreement with previous works that adopt a similar hypothesis for the model of planetesimal fragmentation (Inaba et al. 2003; Kobayashi et al. 2011; Ormel & Kobayashi 2012). Thus, we choose as our baseline case for this work the most favorable simulation in Guilera et al. (2014), which corresponds to the case of initial planetesimals of 100 km radius and a 10 MMSN initial disk mass.

In Fig. 3, we show the time evolution of the core mass and envelope mass of the planet, for the baseline model and for the formation of the planet without considering the planetesimal fragmentation process. We can see that the inclusion of planetesimal fragmentation reduces by ~7% the time to reach the cross-over mass. We can also observe that the planet reaches the cross-over mass at a lower core mass (18% lower). However, the fragmentation model incorporated in Guilera et al. (2014) considered the catastrophic disruption threshold given by Benz & Asphaug (1999) for basalts at impact velocities of 3 km s⁻¹. This is clearly a simplification. On one hand, as we are studying the formation of a giant planet behind the iceline, it should be taken into account that planetesimals are composed of rocks and ices (Lodders 2003). Also, during the process of planet formation, the planetesimal relative velocities are increased due to the gravitational perturbations produced by the growing planet. In the left panel of Fig. 4, we show the increase of eccentricities and inclinations for planetesimals of 100 km radius, and in the right panel their relative velocities, near the planet while it grows. We can see that eccentricities and inclinations are quickly increased near the planet. This leads to an increase of the planetesimal relative velocities from velocities of some meters per second to velocities of about 5 km s⁻¹ (right panel of Fig. 4). These results show that clearly the assumption of a $Q_{\text{D}}^{*}$ given for a constant impact velocity is a simplification in the model of planetesimal fragmentation, and thus, it motivated us to incorporate the following dependencies in the catastrophic disruption threshold.

Fig. 3 Core masses (solid lines) and envelope masses (dashed lines) as function of time. Here, as in Figs. 5, 8, 9, 12, and 13, the red thick lines correspond to the baseline model, while the gray thin lines correspond to the formation of the planet without planetesimal fragmentation.

3.2 Dependencies of Q${}_{\text{D}}^{*}$ with the planetesimal impact velocities and compositions

As we mentioned above, in general models of planet formation that include planetesimal fragmentation consider the catastrophic impact energy threshold for a fixed velocity of basaltic planetesimals. In this section, we explore how the dependencies of the catastrophicdisruption threshold with the planetesimal compositions and with the planetesimal impact relative velocities impact affect the formation of the planet. To do this, we first calculate the formation of the planet considering the same initial parameters as in the baseline model (an embryo located at 5 au in a disk ten times more massive than the MMSN and with initial planetesimals of 100 km radius) but now we implement the dependence of $Q_{\text{D}}^{*}$ with different impact velocities following the approaches described in Sect. 2.2.2. The solid accretion rates are calculated following Eq. (A.6). In Fig. 5, we show the time evolution of the core mass and envelope mass of the planet for the model wherein planetesimal fragmentation is not considered, for the baseline model, and for the case where we consider basaltic planetesimals with $Q_{\text{D}}^{*}$ as a function of the relative velocities among planetesimals. We can see that the cross-over mass, for the last case, is achieved in a longer time than for the baseline model. Moreover, the formation of the giant planet takes ~ 25% more time compared to the case wherein planetesimal fragmentation is not considered. In Fig. 6, we plot the solid accretion rates for the different sizes of planetesimals and small particles for both models, the baseline model, and the model that considers $Q_{\text{D}}^{*}$ as a function of the relative velocities among planetesimals. We can see that the accretion of planetesimals smaller than 100 km and small particles generated by the collisional evolution becomes effective at a longer time for the model where $Q_{\text{D}}^{*}$ is a function of the relative velocities. We can also see that in both models, the total solid accretion rate is always dominated by plantesimalsof 100 km. However, while in the baseline model planetesimals between 1 and 25 km are the products of the fragmentation process that most contribute to the total solid accretion rate (being the planetesimals of ~ 25 km the most important), in the modelwhere $Q_{\text{D}}^{*}$ is a functionof the relative velocities, the most important contribution from the fragments is given by planetesimals between 10 and ~ 16 km. Moreover, in this last case, planetesimals of 1 km make a negligible contribution to the total solid accretion rate. This is due to the fact that these bodies are near the minimum of $Q_{\text{D}}^{*}$ for low impact velocities, which is for this size one order of magnitude lower with respect to the $Q_{\text{D}}^{*}$ at impact velocities of 3 km s⁻¹. We can see from Fig. 7 that relative velocities for planetesimals of 1 m reach 20–30 m s⁻¹ at 1 Myr (when the mass of the core is ~1 M_⊕), and 1 km-sized planetesimals exceed thesevalues rapidly; at 10⁻² Myr (when the mass of the core is only ~0.2 M_⊕) the relative velocity is over 200 m s⁻¹. It is important to remark that, in both cases, their velocities are always below 3 km s⁻¹.

On the other hand, in Fig. 8 we plot the time evolution of the core mass and envelope mass of the planet for a model that considers planetesimals purely composed of ices and adopting $Q_{\text{D}}^{*}$ for ices at impact velocities of 3 km s⁻¹, and for a model wherein planetesimals are composed of 50% of basalts and 50% of ices, and $Q_{\text{D}}^{*}$ is calculated by a linear combinationbetween $Q_{\text{D}}^{*}$ for basalts at 3 km s⁻¹ and $Q_{\text{D}}^{*}$ for ices at 3 km s⁻¹. We also show, for completeness and comparison, the baseline model and the case where planetesimal fragmentation is not considered. We can see that when we considered planetesimals purely composed of ices, the cross-over mass is never achieved. This is due to the fact that for the same fixed impact velocity (in this case 3 km s⁻¹), $Q_{\text{D}}^{*}$ is lower for ices thanfor basalts. Thus, the collisional evolution of the planetesimal population is different, small fragments are disrupted more efficiently, and the accretion rates of such fragments do not compensate for the fragmentation of big bodies. Besides, if we compare $Q_{\text{D}}^{*}$ from Figs. 1 and 2 we can see that $Q_{\text{D}}^{*}$ for ices at 3 km s⁻¹ is lower than $Q_{\text{D}}^{*}$ for basalts at 20–30 m s⁻¹ for planetesimals of ~100 km. Therefore, this difference can explain why in the model that considers $Q_{\text{D}}^{*}$ as a function of the relative velocities for basalts the giant planet core could form during the disk lifetime and, for the model that adopts $Q_{\text{D}}^{*}$ for ices at 3 km s⁻¹, the cross-over mass was not achieved in that time. When planetesimals are composed of 50% of basalts and 50% of ices, and $Q_{\text{D}}^{*}$ is calculated as a linear combination of the corresponding $Q_{\text{D}}^{*}$ for basalts and ices at impact velocities of 3 km s⁻¹, respectively, the cross-over mass is achieved in less than 6 Myr with a core mass of 16 M_⊕. However, the planet achieves the cross-over mass after a longer time, by ~12%, compared to the case without planetesimal fragmentation.

All these results show that the collisional evolution of the planetesimal population and the accretion of the fragments produced in such collisional evolution play an important role in the formation of a giant planet. We point out that in this subsection, the accretion rates of particles with Stokes number less than unity was also computed as in the baseline model, using the prescription derived by Inaba et al. (2001) in the low-velocity regime.

Fig. 4 Left panel: time evolution of the radial profiles of the eccentricity (solid lines) and inclination (dashed lines) for 100 km-sized planetesimals. Right panel: time evolution of the relative velocity between planetesimals with 100 km radius. In brackets we show the total mass of the planet at the given times. The black curve corresponds to the time (and mass) at which the planet reaches the cross-over mass.

Fig. 5 Coremasses (solid lines) and envelope masses (dashed lines) as a function of the time. Gray lines: model without planetesimal fragmentation. Red lines: baseline model. Pink lines: model in which $Q_{\text{D}}^{*}$ is a functionof relative velocities.

Fig. 6 Top panel: solid accretion rates for basaltic planetesimals and small particles of various radii, products of the collisional evolution of planetesimals, for the baseline model. Bottom panel: same as top panel but for the model that considers $Q_{\text{D}}^{*}$ as a function of the planetesimal relative velocities.

Fig. 7 Time evolution of the relative velocity between basaltic planetesimals of the same size for the model that considers $Q_{\text{D}}^{*}$ as a function of the planetesimal relative velocities. Top panel: 1 m-sized planetesimals. Bottom panel: 1 km-sized planetesimals. In brackets we show the total mass of the planet at the given times.

Fig. 8 Core masses (solid lines) and envelope masses (dashed lines) as a function of time. Gray lines: model with no planetesimal fragmentation. Red lines: baseline model. Pink lines: model with $Q_{\text{D}}^{*}$ for ices at impact velocities of 3 km s−1. Blue lines: model where $Q_{\text{D}}^{*}$ is calculated by a linear combination between $Q_{\text{D}}^{*}$ for basalts at 3 km s−1 and $Q_{\text{D}}^{*}$ for ices at 3 km s−1.

3.3 Accretion of small fragments: pebble accretion of second generation

As we mentioned before, in out model we improved the solid accretion rates, incorporating in this work the accretion rates of pebbles (Eq. (A.7)) derived by Lambrechts et al. (2014) for particles with Stokes number less than or equal to unity. In Guilera et al. (2014), the accretion rates of such particles were calculated using the prescription derived by Inaba et al. (2001) in the low-velocity regime (Eq. (A.6)). In that work, we showed that the accretion rates for particles smaller than ~ 1 m could become greater than the pebble accretion rates (see Guilera et al. 2014, for the details). Thus, in order to calculate in a more accurate way the accretion of small particles, we implemented the pebble accretion rates mentioned above.

In Fig. 9, we show the time evolution of the core mass and the envelope mass for the case where planetesimals are purely composed of basalts adopting the corresponding $Q_{\text{D}}^{*}$ at impact velocities of 3 km s⁻¹, and where small particles with Stokes number less or equal to unity are accreted as pebbles using Eq. (A.7), while in the baseline model the accretion rate for planetesimals (Eq. (A.6)) is applied also for small particles. We also plot for comparison the baseline model and the case wherein planetesimal fragmentation is not considered. We can see that a more accurate treatment for the accretion of small particles delays the formation of the giant planet by about 30%, despite the fact that pebbles do not mostly contribute to the total solid accretion rate (Eq. (A.7)).

This effect is discussed in Guilera et al. (2014) where the rates of small fragments using the prescriptions of Inaba et al. (2001) and Lambrechts & Johansen (2012) are compared. Following that analysis we here obtain the accretion rates of Inaba et al. (2001) in terms of the pebble accretion rates of Lambrechts et al. (2014): $$\dot{M}=\{\begin{array}{l}5.65\sqrt{{\tilde{R}}_{\text{C}}/R_{\text{H}}}{\dot{M}}_{\text{H}},\text{\hspace{1em}if } 0.1\le St<1,\hfill \\ 5.65\sqrt{{\tilde{R}}_{\text{C}}/R_{\text{H}}}{\left(\frac{St}{0.1}\right)}^{-2/3}{\dot{M}}_{\text{H}2},\text{\hspace{1em}if } St<0.1,\hfill \end{array}$$(2)

where ${\tilde{R}}_{\text{C}}$ is the enhanced radius due to the planet’s gaseous envelope, ${\dot{M}}_{\text{H}}=2R_{\text{H}}^2{\Sigma }_p{\Omega }_P$ and ${\dot{M}}_{\text{H2}}={\left(\frac{St}{0.1}\right)}^{2/3}{\dot{M}}_{\text{H}}$, where Ṁ_H for 0.1 ≤ St < 1 and Ṁ_H2 for St < 0.1 are the pebble accretion rate dM_C∕dt of Eq. (A.7). We can see from Fig. 15 of Guilera et al. (2014) that when the planet’s core mass becomes larger than ~0.2 M_⊕ planetesimal accretion rates are higher than pebble accretion rates for fragments with r_P ≲ 1 m. In our work this is the case for particles with 0.1 ≤ St < 1, but for fragments with St < 0.1, the difference between the planetesimal and pebble accretion rates is even higher due to the factor ${\left(\frac{St}{0.1}\right)}^{-2/3}$. Moreover, we introduced a factor β = min(1, R_H∕H_P) to include a reduction in the pebble accretion rates if the scale height of the pebbles becomes greater than the Hill radius of the planet; this is not taken into account in the planetesimal accretion rate for small fragments in the baseline model. In Fig. 10 we show R_H∕H_P as a function of time, where we can observe that objects between 1 and 10 cm have values of R_H ∕H_P below 1. The quantified differences are shown in Fig. 11 where we present the comparison of the solid accretion rates for small particles with St ≲ 1 of the baseline case and the simulation including pebble accretion. We can see from Fig. 11 that the accretion rates in the baseline model are higher than in the simulation that takes into account the accretion of pebbles as explained before.

Fig. 9 Core masses (solid lines) and envelope masses (dashed lines) as a function of time. Gray lines: model with no planetesimal fragmentation. Red lines: baseline model. Pink lines: model with $Q_{\text{D}}^{*}$ for basalts at impact velocities of 3 km s−1 and where small particles with Stokes numbers less than or equal to unity are accreted as pebbles of second generation.

Fig. 10 Ratio between the Hill radius and the scale height of small particles, with Stokes number less than or equal to unity, for the baseline model (solid lines) and for the simulation including pebble accretion (dashed lines) as a function of time.

Fig. 11 Solid accretion rates for small particles, product of the collisional evolution of planetesimals, for the baseline model (solidlines) and for the simulation including pebble accretion (dashed lines) as a function of time. Small particles of these sizes have Stokes number less than or equal to unity.

3.4 Global model

In this section, we compare the results for the formation of the giant planet obtained with the baseline model and with the global model. The global model includes all the improvements in the calculation of $Q_{\text{D}}^{*}$, that is, the dependencies of $Q_{\text{D}}^{*}$ with the planetesimal impact velocities and their compositions described in Sect. 3.2, and the accretion of second generation pebbles described in Sect. 3.3.

In the global model where we include all the new effects on $Q_{\text{D}}^{*}$ and a more accurate and realistic treatment for the accretion of the small particles, the core does not achieve the cross-over mass within the 6 Myr as shown in Fig. 12. We point out that this model also includes the different regimes for the calculation of the relative velocities (and the impact rates) discussed in Sect. 2.2.1. However, as the planet quickly excites the relative velocities of the near planetesimals, they are in general in the dispersion regime, and thus the Keplerian shear regime does not play a relevant role. Finally, we remark that despite the fact that the planet does not achieve the cross-over mass, a solid core of a few Earth masses, up to five Earth masses, is formed more quickly compared to the case where planetesimal fragmentation is not considered. Ikoma et al. (2000) and Hubickyj et al. (2005) showed that a reduction in the grain opacity in the planet envelope, as well as the pollution of the envelope (due to evaporated materials of icy planetesimals in the envelope, see Hori & Ikoma 2011), could reduce the mass of the core at which the planet reaches the cross-over mass, and thus reduce the formation time. We will explore this possibility in the next section.

Fig. 12 Coremasses (solid lines) and envelope masses (dashed lines) as a function of time. Gray lines: model without planetesimal fragmentation. Red lines: baseline model. Pink lines: model with all the improvements in the calculation of $Q_{\text{D}}^{*}$ (the dependencies of $Q_{\text{D}}^{*}$ with the planetesimal impact velocities and compositions) and the pebble accretion rates.

3.5 Reduced grain opacities

In the previous sections we studied the process of giant planet formation considering fragmentation of planetesimals including, for the catastrophic disruption threshold, the dependencies with relative velocity and composition of the colliding bodies. We also include the accretion rates for small fragments called pebbles. We analyzed every improvement separately and finally all together and we found that, in all the cases, the formation of the giant planet core is slower.

Motivated by this result we incorporate a physical phenomenon that could act in an opposite way, accelerating the formation of a giant planet, in this case a reduction in the grain opacities of the planet envelope. The grain opacities in the planet envelope play an important role in the formation of a giant planet. Ikoma et al. (2000) found that if the grain opacities are small enough, a small-sized giant planet core can capture significant amounts of gas. Later, Hubickyj et al. (2005) found that if the grain opacities are assumed to be 2% of the interstellar medium value, the planet can reach the cross-over mass with a core mass between 5 and 10 M_⊕. More recently, Hori & Ikoma (2010) showed that a core of only ~2 M_⊕ is able to capture enough gas to form a giant planet if the accreting envelope is grain-free. Thus, following these works, especially the work of Hubickyj et al. (2005), we compute a set of simulations with and without planetesimal fragmentation, reducing the grain opacities in the envelope to up to 2% of the interstellar medium value. For the case where we consider the planetesimal fragmentation, we adopt the global model described above.

In Fig. 13, we present the time evolution of the core mass and the envelope mass of the planet for the models that consider reduced grain opacities with and without planetesimal fragmentation and the baseline case. When we calculate the formation of the giant planet without considering the planetesimal fragmentation process, the time at which the planet achieves the cross-over mass is reduced by more than 50% when the grain opacities in the planet envelope are reduced. However, the cross-over mass remains practically similar. These results are in very good agreement with the previous results of Hubickyj et al. (2005). A reduction in the envelope opacity allows the planet to release more efficiently the heat generated by the accretion of the planetesimals, and as a consequence, the gas accretion becomes more efficient. We can see this effect in Fig. 14 where we plot the envelope mass as a function of the core mass for the two models without planetesimal fragmentation. At a fixed core mass, the model with the reduced grain opacities has a greater value for the envelope mass.

Finally, we note that when planetesimal fragmentation is considered, the planet achieves the cross-over mass at a time (1.66 Myr) and with a core mass (16.44 M_⊕) that are lower than for the case where planetesimal fragmentation is not considered (~7 and 24% lower in time and core mass, respectively). We associate this result to the fact that the amount of envelope mass in the planet at low-mass cores plays an important role, significantly enhancing the capture radius of the planet.

Fig. 13 Core masses (solid lines) and envelope masses (dashed lines) as a function of time. Gray lines: model without planetesimal fragmentation. Black lines: model without planetesimal fragmentation but considering grain opacities reduction (G.O.R). Red lines: baseline model. Pink lines: model taking into account all the improvements in the calculation of $Q_{\text{D}}^{*}$ (the dependencies of $Q_{\text{D}}^{*}$ with the planetesimal impact velocities and compositions), pebble accretion rates, and considering reduced grain opacities.

Fig. 14 Envelope mass as a function of the core mass for the two cases where planetesimal fragmentation is not considered. We plot the core mass until 10 M⊕ in order to show more clearly the differences at low-mass core values. The black line represents the case where grain opacities are not reduced, while the red line represents the case where grain opacities are reduced in the planet envelope.

4 Comparison with previous work

A similar model of giant planet formation was developed by Chambers (2014). They carried out numerical simulations of oligarchic growth including pebble accretion and planetesimal fragmentation. One of the major differences between their model and ours is that they consider that all the collisions involving targets smaller than their minimum bin size are assumed to coagulate and grow rapidly when they approach such minimum size, while we assume that the mass distributed under $r_{\text{P}}^{\text{min}}=1$ cm is lost. In our fragmentation model most of the mass is in the largest fragments. As discussed in Guilera et al. (2014), the results are analogous when the minimum bin size $r_{\text{P}}^{\text{min}}$ is reduced to lower values. Also, Chambers (2014) uses a fixed value of $Q_{\text{D}}^{*}$ (for basalts at 3 km s⁻¹) and includes an additional factor to reduce 10 and 100 times the strength for all the sizes of the bodies to see how this could change the growth of the embryo. He found the most significant differences in the case where $Q_{\text{D}}^{*}$ is reduced by a factor of 100, changing the growth track of the embryo (shown in Fig. 20 of his paper). Even though we cannot compare directly this result with ours since our models are different, we agree that the relative importance of pebble accretion for giant planet formation will depend strongly on the rate at which pebbles are generated from planetesimal-planetesimal collisions, depending in particular on the impact parameters.

Anothermodel of giant planet formation that includes pebble accretion is that of Alibert et al. (2018). They developed a model of the formation of Jupiter that could explain the constraints that the cosmochemical evidence gives on the existence of two main reservoirs of small bodies that remain separated during the early solar system (Kruijer et al. 2017). They conclude that Jupiter formed in a three-step process. First, Jupiter’s core grew by the accretion of pebbles. Second, the pebble accretion stopped and the core continued growing by the accretion of small planetesimals. Third, the core was massive enough for the runaway gasaccretion to start. It is important to highlight that our initial conditions correspond to the beginning of the oligarchic growth where almost all the solid material is contained in planetesimals, in contrast to Alibert et al. (2018) where they start their simulations with a first generation of pebbles. Moreover, they do not include a fragmentation model since they estimate, for the start of their second stage, the amount of fragments that could have been produced during their first stage of pebble accretion. For this estimation, they used a fixed value of $Q_{\text{D}}^{*}$ for 100 km-sized planetesimals. Different initial conditions for giant planet formation including a first generation of pebbles and a variation of the initial size of planetesimals will be studied in future works.

5 Summary and conclusions

We studied the formation of a giant planet considering the collisional evolution of the planetesimal population, including the dependence of the catastrophic impact energy threshold with composition and relative velocities of the colliding planetesimals. We also included the pebble accretion rates for small particles produced by the collisional cascade, which are pebbles of second generation strongly coupled to the gas, to calculate in a more accurate way than in our baseline case. We showed that as the planet grows, the relative velocities among the planetesimals near the planet are increased from velocities of the order of meters per second to velocities up to ~ 4.5 km s⁻¹ before the planet reaches the cross-over mass.

We analyzed the improvements incorporated in the calculation of $Q_{\text{D}}^{*}$ and the inclusion of the pebble accretion rates one at each time, and finally all the improvements together. We point out that, for each improvement in $Q_{\text{D}}^{*}$, as the collisional cascade increases, the timescale to reach the cross-over mass is delayed. We also found that in all the cases, the formation of the core is less efficient compared with the case where planetesimal fragmentation is not considered, and with our baseline model.

When we include all the improvements at once there is a sum of effects that inhibit the formation of the giant planet core. On one hand, pebble accretion rates are lower than the planetesimal accretion rates for the small particles adopted in the baseline model. This effect contributes to slow down the growth of the giant planet core. On the other hand, when we include the dependence of $Q_{\text{D}}^{*}$ with impact velocity and composition we first interpolate between the curves of $Q_{\text{D}}^{*}$ for a given velocity obtaining $Q_{\text{D}}^{*}$ for each pure material. Then, we make a linear combination between these values where we assume that planetesimals are half basalts and half ices. For low impact velocities basalts are weaker than ices (except for ~100 km-sized bodies) and for higher velocities, ice is weaker than basalts. These effects are averaged since the composition of planetesimals is 50% basalts and 50% ices. Moreover, it is important to remark that the catastrophic disruption thresholds for ices at impact velocities of 0.5 and 3 km s⁻¹ are lower for ~100 km-sized planetesimals than $Q_{\text{D}}^{*}$ for basalts at any of the three impact velocities (20–30 m s⁻¹, 3 and 5 km s⁻¹) for that size of bodies. Regardless of whether the planet arrives at the cross-over mass or not within the disk lifetime, in all the cases, planetesimal fragmentation favors the relative rapid growth of the giant planet core up to a few Earth masses compared to the case where planetesimal fragmentation is not considered.

This last result motivated us to explore alternatives for the formation of a giant planet with less massive cores. In this way, we followed the work of Hubickyj et al. (2005). These authors showed that if the grain opacity of the planet envelope is reduced to 2% of the interstellar medium value, the planet can reach the cross-over mass with a core mass between 5 and 10 M_⊕. Thus we computed two new simulations, reducing the grain opacity of the envelope up to 2% of the interstellar medium value, considering planetesimal fragmentation (including all the improvements in our model) and without planetesimal fragmentation. For this last case, we found that the planet reached the cross-over mass at a time significantly lower (by a factor two) with respect to the case without planetesimal fragmentation but not considering reduced grain opacities. However, the mass of the core at which the cross-over mass is achieved remained practically similar. For the case where planetesimal fragmentation and reduced grain opacities are included, the planet reached the cross-over time at a lower time and with a less massive core compared to the case where reduced grain opacities are included but without planetesimal fragmentation. We associated this result to the fact that if the opacity of the planet envelope is reduced, the planet releases more efficiently the luminosity through the envelope, and then, the gas accretion becomes more efficient. Thus the planet accretes significantly more gas for a low-mass core, enhancing the radius of capture of the planet and increasing the solid accretion rate.

Hori & Ikoma (2011) studied the formation of a giant planet where the gas envelope is polluted by icy planetesimals that are dissolved in the envelope. They found that the increase in the molecular weight and the reduction in the adiabatic temperature gradient produced by the pollution of planetesimals significantly reduce the mass of the core at which the planet achieves the cross-over mass. More recently, Venturini et al. (2015, 2016) found similar results, computing for the first time self-consistent models of giant planet formation that include the effect of envelope enrichment by the pollution of planetesimals dissolved as they enter the planet envelope. The increment in the envelope molecular weight and the reduction in adiabatic temperature gradient mean that the envelope layers are compressed more efficiently and the gas accretion is significantly increased. We will incorporate this phenomenon in a future work.

Acknowledgements

We would like to thank John Chambers and the anonymous referee for the comments and suggestions that helped us to improve this paper. This work was supported by the PIP 112-201501-00699CO grant from Consejo Nacional de Investigaciones Científicas y Técnicas, and by PIDT 11/G144 grant from Universidad Nacional de La Plata, Argentina. O.M.G is also supported by the PICT 2016-0023 from ANPCyT, Argentina. O.M.G. also acknowledges being hosted as an invited researcher by IA-PUC.

Appendix A Accretion model

A.1 Evolution of the envelope and planetesimal radial migration

The evolution of the envelope is calculated solving the standard equations of the stellar evolution theory. The following equations correspond to the conservation of mass, the hydrostatic equilibrium, the energetic balance, and the energy transport: $$\frac{\partial r}{\partial m_r}=\frac{1}{4\pi r^2\rho },$$(A.1) $$\frac{\partial P}{\partial m_r}=-\frac{G{m}_r}{4\pi r^4},$$(A.2) $$\frac{\partial L_r}{\partial m_r}={\varepsilon }_{\text{pl}}-T\frac{\partial S}{\partial t},$$(A.3) $$\frac{\partial T}{\partial m_r}=-\frac{G{m}_{r}T}{4\pi r^{4}P}\nabla ,$$(A.4)

where G is the universal gravitational constant, ρ is the density of the envelope, S is the entropy per unit mass, ɛ_pl is the energy release rate due to the accretion of planetesimals, and $\nabla \equiv \frac{\text{dln}T}{\text{dln}P}$ is the dimensionless temperature gradient, which depends on whether the energy is carried by radiation or by convection (see Fortier et al. 2009; Guilera et al. 2010, for details).

In our simulations we consider 2500 radial bins logarithmically equally spaced in a protoplanetary disk defined between 0.4 and 20 AU. The particles migrate inward in the protoplanetary disk due to gas drag, where the radial migration velocity for each drag regime is given by $$v_{\text{mig}}=\{\begin{array}{l}-\frac{2a\eta }{t_{\text{stop}}}\left[\frac{S{t}^2}{1+S{t}^2}\right],\text{\hspace{1em} Epstein regime}\hfill \\ \hfill \\ -\frac{2a\eta }{t_{\text{stop}}}\left[\frac{S{t}^2}{1+S{t}^2}\right],\text{\hspace{1em} Stokes regime}\hfill \\ \hfill \\ -\frac{2a\eta }{t_{\text{stop}}},\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{1em} quadratic regime},\hfill \end{array}$$(A.5)

being a the semi-major axis, η the ratio of the gas velocity to the local Keplerian velocity, St = t_stopω_k the Stokes number, with t_stop the stopping time depending on the drag regime (Epstein, Stokes or quadratic regime) and ω_k the Keplerian frequency. The rest of the model is explained in detail in Guilera et al. (2010, 2014).

A.2 Solid accretion rates

For planetesimals, we use the accretion rates given by Inaba et al. (2001), $$\frac{\text{d}{M}_{\text{C}}}{\text{d}t}=R_{\text{H}}^2{\Sigma }_{\text{P}}{\Omega }_{\text{P}}{P}_{\text{coll}},\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}when\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}St\ge 1,$$(A.6)

where M_C is the core’s mass, R_H is the planet’s Hill radius, Σ_P is the surface density of solids at the location of the planet, Ω_P is the Keplerian frequency atthe planet location, and P_coll is the collision probability, which is a function of the core radius R_C, the Hill radius of the planet, and the relative velocity between the planetesimals and the planet v_rel, thus P_coll = P_coll(R_C, R_H, v_rel). In fact, as we also consider the drag force that planetesimals experience on entering the planetary envelope (following Inaba & Ikoma 2003), the collision probability is a function of the enhanced radius ${\tilde{R}}_{\text{C}}$ instead of R_C (Guilera et al. 2014).

For the pebbles of second generation, we use the pebble accretion rates given by Lambrechts et al. (2014), $$\frac{\text{d}{M}_{\text{C}}}{\text{d}t}=\{\begin{array}{l}2\beta R_{\text{H}}^2{\Sigma }_{\text{P}}{\Omega }_{\text{P}},\text{\hspace{1em}if } 0.1\le St<1,\hfill \\ 2\beta {\left(\frac{St}{0.1}\right)}^{2/3}{R}_{\text{H}}^2{\Sigma }_{\text{P}}{\Omega }_{\text{P}},\text{\hspace{1em}if } St<\mathrm{0.1.}\hfill \end{array}$$(A.7)

In Eq. (A.7), we introduce the factor β = min(1, R_H∕H_p) in order to take into account a reduction in the pebble accretion rates if the scale height of the pebbles, H_p, becomes greater than the Hill radius of the planet. This can happen if the vertical turbulent dispersion of small particles becomes significant. The scale height of the solids at a given distance from the central star is given by (Youdin & Lithwick 2007) $$H_{\text{p}}=H_{\text{g}}\sqrt{\frac{\alpha }{\alpha +St}},$$(A.8)

being H_g the scale height of thedisk, and α the dimensionless Shakura & Sunyaev viscosity-parameter (Shakura & Sunyaev 1973); we adopt α =10⁻³.

Appendix B Fragmentation model

According to our model, a collision between a target of mass M_T and a projectile of mass M_P results in a remnant of mass M_R given by $$M_{\text{R}}=\{\begin{array}{ll}\left[-\frac{1}{2}\left(\frac{Q}{Q_{\text{D}}^{*}}-1\right)+\frac{1}{2}\right](M_{\text{T}}+M_{\text{P}}),\hfill & \text{\hspace{1em}if }Q<Q_{\text{D}}*,\hfill \\ \left[-0.35\left(\frac{Q}{Q_{\text{D}}^{*}}-1\right)+\frac{1}{2}\right](M_{\text{T}}+M_{\text{P}}),\hfill & \text{\hspace{1em}if }Q>Q_{\text{D}}*,\hfill \end{array}$$(B.1)

where Q is the collisional energy per unit target mass and $Q_{\text{D}}^{*}$ is the catastrophic impact energy threshold per unit target mass required to fragment and disperse half of the target mass. Usually, $Q_{\text{D}}^{*}$ is a function of the target’s radius. However, it is important to remark here that in our model (following Morbidelli et al. 2009) $Q_{\text{D}}^{*}$ has to be calculated with an effective radius r_eff = 3(M_T + M_P)∕4πρ, being ρ the planetesimal density. Besides this, in our model the mass loss in the collision, which we define as (M_T + M_P − M_R), is distributed between the minimum mass bin considered and the mass bin corresponding to the biggest fragment M_F given by $$M_{\text{F}}=8\times {10}^{-3}\left[\frac{Q}{Q_{\text{D}}^{*}} e^{-{(Q/4Q_{\text{D}}^{*})}^2}\right] (M_{\text{T}}+M_{\text{P}}).$$(B.2)

As in Guilera et al. (2014), we note again that for some super-catastrophic collisions (which occur when M_R ≪ M_T + M_P), M_F > M_R. For such collisions we assume that M_F = 0.5M_R.

Unlike Morbidelli et al. (2009), the fragments are distributed following a power-law distribution given by (Kobayashi et al. 2011; Ormel & Kobayashi 2012) $$\frac{\text{d}n}{\text{d}m}\propto m^{-5/3},$$(B.3)

meaning that most of the mass is distributed in the larger fragments. In Guilera et al. (2014), we found that massive core formation is favored when the exponent of the power-law mass distribution is lower than 2, thus in this work we analyze the most favorable scenarios, taking the value of 5/3 for the exponent. To calculate the growth of the planet, we adopt a discrete planetesimal size distribution using 36 size (or mass) bins logarithmic equally spaced between 1 cm and 100 km, where initially all the solid mass is in the non-porous planetesimals of 100 km radius. However, when we calculate the planetesimal fragmentation process, we extrapolate the planetesimal size distribution two orders of magnitude below the minimum size $r_{\text{P}}^{\text{min}}$ of the main model to avoid the accumulation of spurious mass in the smaller fragments. Therefore, only the mass ejected from the collision distributed between the mass of the larger fragment and the minimum size considered ($r_{\text{P}}^{\text{min}}=1$ cm) is taken into account to calculate later the solid accretions rates; that is, we assume that the mass distributed below 1 cm is lost.

The feeding zone of the embryo extends to four Hill radii at either side of the embryo. Adopting 2500 radial bins along the protoplanetary disk guarantees that there are at least ten radial bins between R_P − 4R_H and R_P + 4R_H at the beginning of the simulation. We define the width of the fragmentation zone as twice the feeding zone, that is, eight times the embryo Hill radius at both sides of the embryo. The excitation of eccentricities and inclinations decay with the distance to the embryo, especially outside the feeding zone. Then, our definition of the fragmentation zone guarantees that collisions are well determined within this zone. The amount of bins increases as the core mass grows, for example, in the baseline case (see Sect. 3.1), when the planet reaches the cross-over mass, the fragmentation zone has ~ 600 radial bins.

B.1 Velocities and probabilities of collision regimes

Following Greenberg et al. (1991) we adopt three different regimes and their transitions, regime A: dominance by random motion; regime B: dominance by Keplerian shear motion; regime C: Keplerian shear dominance in a very thin disk. For the dispersion regime (dominance by random motion), where Keplerian behavior is unimportant, we adopt the impact rate given by Eq. (1).

The transition between regime A and regime B is given by $$\frac{(a_{\text{P}}+a_{\text{T}})}{2}\frac{(e_{\text{P}}+e_{\text{T}})}{2}=2.5R_{{\text{H}}_{\text{T}}},$$(B.4)

where $R_{{\text{H}}_{\text{T}}}$ is the Hill radius of the target, and a_P, e_P and a_T, e_T are the semi-major axis and eccentricities of the projectile and the target, respectively. Then, for larger values of e and a the system is dominated by random motion (regime A) and for smaller values of these orbital parameters the system is in the Keplerian shear regime (regime B). If a_pi_p < R_G, being R_G the target’s gravitational diameter and i_p the inclination of the projectile’s orbit, the particles are in regime C, wherein still dominates the Keplerian shear motion but the system is two dimensional.

For regimes B and C the relative velocity is given by $$v=0.58{(2{\mu }^{1/15}-1.27)}^{1/2}\Delta a,$$(B.5)

where μ = M_p∕M_⊙ and $\Delta a=2.5R_{{\text{H}}_{\text{T}}}$. The impact rate of regime B follows $$\text{Impact rate}{|}_{\text{B}}=\pi R_{\text{P}}^2{\left(1+b\frac{V_{\text{esc}}^2}{v^2}\right)}^{1/2}\frac{\sigma {\left(2.5R_{{\text{H}}_{\text{T}}}\right)}^{2}1.125\omega }{\frac{a_{\text{P}}+a_{\text{T}}}{2}4a_{\text{p}}{i}_{\text{p}}{\mu }^{2/5}{M}_{\text{P}}},$$(B.6)

where σ is the projectile’s surface density and ω is the Keplerian frequency using a = (a_P + a_T)∕2.

Finally, the impact rate of regime C is given by $$\text{Impact rate}{|}_{\text{C}}=R_{\text{P}}{\left(1+b\frac{V_{\text{esc}}^2}{v^2}\right)}^{1/2}\frac{\sigma {(2.5R_{{\text{H}}_{\text{T}}})}^{2}1.125\omega }{\frac{(a_{\text{P}}+a_{\text{T}})}{2}{\mu }^{2/5}{M}_{\text{P}}}.$$(B.7)

B.2 Catastrophic impact energy threshold

From SHP simulations, Benz & Asphaug (1999) found that $Q_D^{*}$ can be expressed by the functional form $$Q_{\text{D}}^{*}=Q_0{\left(\frac{R_{\text{T}}}{1 \text{cm}}\right)}^a+B\rho {\left(\frac{R_{\text{T}}}{1 \text{cm}}\right)}^b,$$(B.8)

being Q₀, B, a, and b parameters that depend on the properties of the material and on the impact velocity over the target, and ρ the density of the non-porous planetesimals (in this work we adopted ρ = 1.5 g cm⁻³). They performed simulations for non-porous basalts at 3 and 5 km s⁻¹, and for non-porous ices at 0.5 and 3 km s⁻¹. Later, Benz (2000), performed new SPH simulations to calculate $Q_{\text{D}}^{*}$ for non-porous basalts at low impact velocities (between 20 and 30 m s⁻¹) finding that targets impacted at such low impact velocities are weaker than targets impacted at greater velocities. Despite the fact that Benz (2000) did not provide a functional form for $Q_{\text{D}}^{*}$, we fit the results of their simulations adopting the same functional form proposed by Benz & Asphaug (1999). Table B.1 summarizes the values of the free parameters that determine the values of $Q_{\text{D}}^{*}$ used in this work for different types of materials and different impact velocities.

References

All Tables

All Figures

	Fig. 1 Catastrophic disruption thresholds for basalt targets at impact velocities of 20–30 m s−1 (Benz 2000), 3 and 5 km s−1 (Benz & Asphaug 1999). The blue points correspond to the discrete data extracted from Benz (2000).
In the text

	Fig. 2 Catastrophic disruption thresholds for icy targets at impact velocities of 0.5 and 3 km s−1 (Benz & Asphaug 1999).
In the text

	Fig. 3 Core masses (solid lines) and envelope masses (dashed lines) as function of time. Here, as in Figs. 5, 8, 9, 12, and 13, the red thick lines correspond to the baseline model, while the gray thin lines correspond to the formation of the planet without planetesimal fragmentation.
In the text

Fig. 4 Left panel: time evolution of the radial profiles of the eccentricity (solid lines) and inclination (dashed lines) for 100 km-sized planetesimals. Right panel: time evolution of the relative velocity between planetesimals with 100 km radius. In brackets we show the total mass of the planet at the given times. The black curve corresponds to the time (and mass) at which the planet reaches the cross-over mass.

In the text

	Fig. 5 Coremasses (solid lines) and envelope masses (dashed lines) as a function of the time. Gray lines: model without planetesimal fragmentation. Red lines: baseline model. Pink lines: model in which $Q_{\text{D}}^{*}$ is a functionof relative velocities.
In the text

	Fig. 6 Top panel: solid accretion rates for basaltic planetesimals and small particles of various radii, products of the collisional evolution of planetesimals, for the baseline model. Bottom panel: same as top panel but for the model that considers $Q_{\text{D}}^{*}$ as a function of the planetesimal relative velocities.
In the text

	Fig. 7 Time evolution of the relative velocity between basaltic planetesimals of the same size for the model that considers $Q_{\text{D}}^{*}$ as a function of the planetesimal relative velocities. Top panel: 1 m-sized planetesimals. Bottom panel: 1 km-sized planetesimals. In brackets we show the total mass of the planet at the given times.
In the text

Fig. 8 Core masses (solid lines) and envelope masses (dashed lines) as a function of time. Gray lines: model with no planetesimal fragmentation. Red lines: baseline model. Pink lines: model with $Q_{\text{D}}^{*}$ for ices at impact velocities of 3 km s−1. Blue lines: model where $Q_{\text{D}}^{*}$ is calculated by a linear combination between $Q_{\text{D}}^{*}$ for basalts at 3 km s−1 and $Q_{\text{D}}^{*}$ for ices at 3 km s−1.

In the text

	Fig. 9 Core masses (solid lines) and envelope masses (dashed lines) as a function of time. Gray lines: model with no planetesimal fragmentation. Red lines: baseline model. Pink lines: model with $Q_{\text{D}}^{*}$ for basalts at impact velocities of 3 km s−1 and where small particles with Stokes numbers less than or equal to unity are accreted as pebbles of second generation.
In the text

	Fig. 10 Ratio between the Hill radius and the scale height of small particles, with Stokes number less than or equal to unity, for the baseline model (solid lines) and for the simulation including pebble accretion (dashed lines) as a function of time.
In the text

	Fig. 11 Solid accretion rates for small particles, product of the collisional evolution of planetesimals, for the baseline model (solidlines) and for the simulation including pebble accretion (dashed lines) as a function of time. Small particles of these sizes have Stokes number less than or equal to unity.
In the text

	Fig. 12 Coremasses (solid lines) and envelope masses (dashed lines) as a function of time. Gray lines: model without planetesimal fragmentation. Red lines: baseline model. Pink lines: model with all the improvements in the calculation of $Q_{\text{D}}^{*}$ (the dependencies of $Q_{\text{D}}^{*}$ with the planetesimal impact velocities and compositions) and the pebble accretion rates.
In the text

Fig. 13 Core masses (solid lines) and envelope masses (dashed lines) as a function of time. Gray lines: model without planetesimal fragmentation. Black lines: model without planetesimal fragmentation but considering grain opacities reduction (G.O.R). Red lines: baseline model. Pink lines: model taking into account all the improvements in the calculation of $Q_{\text{D}}^{*}$ (the dependencies of $Q_{\text{D}}^{*}$ with the planetesimal impact velocities and compositions), pebble accretion rates, and considering reduced grain opacities.

In the text

	Fig. 14 Envelope mass as a function of the core mass for the two cases where planetesimal fragmentation is not considered. We plot the core mass until 10 M⊕ in order to show more clearly the differences at low-mass core values. The black line represents the case where grain opacities are not reduced, while the red line represents the case where grain opacities are reduced in the planet envelope.
In the text