© ESO, 2014

1. Introduction

According to the core-accretion model (Lissauer & Stevenson 2007), the formation of a giant planet occurs through a sequence of events: initially, the dust (particles of ~μm sizes) collapses onto the protoplanetary disk mid-plane. Then these particles agglomerate by different mechanisms, which are still debated, which leads to the formation of planetesimals (objects of between hundreds of meters and hundreds of kilometers). These planetesimals grow by mutual accretions until some bodies begin to separate themselves from the population of planetesimals (planetary embryos, objects with sizes of a few thousand kilometers). The embryo gravitational excitation, which is higher than that of the planetesimals, limits their growth, and embryos are the only bodies that grow by accretion of planetesimals. As the embryos grow, they bound a gaseous envelope. Initially, the planetesimal accretion rate is much higher than the gas accretion rate. However, when the mass of the envelope reaches a critical value on the order of the mass of the solid core, a gaseous runaway growth process starts. Finally – by mechanisms that are also still debated – the planet stops to accrete gas and evolves, contracting and cooling at constant mass.

Therefore, it is important to study the evolution of the population of planetesimals together with the process of accretion and planet formation. The evolution of the population of planetesimals is a complex phenomenon. Probably the most relevant phenomena are planetesimal accretion by the embryos, migration due to gas drag produced by the gaseous component of the protoplanetary disk, collisional evolution through gravitational excitations produced by embryos, and planetesimal dispersion and gap openings.

As embryos grow, they increase the relative velocities of the planetesimals causing planetesimal fragmentation. After successive disrupting collisions also called collisional cascade planetesimals become smaller. Inaba et al. (2003), Kobayashi et al. (2011), and Ormel & Kobayashi (2012) found that a significant amount of mass, which remains in small fragments as a product of the collisions between planetesimals, may be lost by migration due to gas drag. This means that the planetesimal fragmentation seems to play an important role in forming the cores of giant planets. Moreover, as embryos grow, they begin to bind the surrounding gas. Initially, embryo gaseous envelopes are less massive, but wide spread. These envelopes cause a loss of the planetesimal kinetic energy, significantly increasing the capture cross-section of the planets. The smaller planetesimals of the distribution suffer more from these effects. Accordingly, while smaller fragments have higher migration rates due to gas drag, they are more efficiently accreted by the planet. There is a strong competition between the time scales of migration and accretion of small fragments generated by planetesimal fragmentation. Therefore, it is important to study in detail whether the generation of small fragments, the products of planetesimal fragmentation, favors or inhibits giant planet formation.

In previous works (Guilera et al. 2010, 2011) we developed a model for giant planet formation that calculates their formation immersed in a protoplanetary disk that evolves in time. In these previous works the population of planetesimals evolved by the accretion of the embryos and by planetesimal migration. In this new work, we incorporated the fragmentation of planetesimals to study whether this phenomenon produces significant changes, which would make it a primary key to consider in the process of planetary formation.

This work is organized as follows: in Sect. 2 we introduce some improvements to our previous model; in Sect. 3, we explain in detail our planetesimal fragmentation model; Sect. 4 shows the results of the role of planetesimal fragmentation on the growth of an embryo located at 5 au; finally, in Sect. 5 we present the conclusions from our results.

2. Improvements to our previous model

Our model that describes the evolution of the protoplanetary disk is based on the works of Guilera et al. (2010, 2011) with some minor improves incorporated. We used an axisymmetric protoplanetary disk characterized by a gaseous and a solid component. The gaseous component is represented by a 1D grid for the radial coordinate, while the solid component is represented by a 2D grid, where one dimension is for the radial coordinate and the other one is for the different planetesimal sizes. Some quantities are only functions of the radial coordinate (R), such as the gas surface density Σ_g(R), while some others are also functions of the planetesimal sizes, such as the planetesimal surface density Σ_p(R,r_p).

2.1. Planetesimal size distribution

We changed the treatment of the planetesimal size distribution from our previous work. We now considered that the different radii are logarithmically equally spaced, so the j species is given by (1)where r_pN and r_p1 are the largest and smallest radii of the size distribution, and N is the number of size bins considered. If m_pj is the mass between m_{pj − 1 /2} and m_{pj + 1 /2}, and adopting that the planetesimal mass distribution is represented by a power law (d), m_pj is given by (2)where we use that m_pj = Δ^{3(j − 1)/(N − 1)}m_p1 with Δ = r_pN/r_p1. In the same way, the total mass is given by (3)The amount of mass (with respect to the total mass) corresponding to the j species is given by (4)Finally, the planetesimal surface density corresponding to the planetesimals of radius r_pj is obtained by multiplying p_j and the total surface density of solids. Then, we treated each planetesimal size independently. In this approach we can use only one planetesimal size (for this case p = 1) or a discrete number (N) of bins to approximate the continuous planetesimal size distribution.

2.2. Evolution of planetesimal eccentricities, inclinations, and velocity migrations

As in our previous works, we considered that the evolution of the eccentricities and inclinations of the planetesimals are governed by two main processes: the embryo gravitational excitations and the damping due to the gas drag.

The embryo stirring rates of the eccentricities and inclinations are given by (Ohtsuki et al. 2002) where M_P is the mass of the planetary embryo, M_⋆ is the mass of the central star, b is the full width of the feeding zone of the planetary embryo in terms of its Hill radius, and P_orb is the orbital period of the embryo. Finally, P_stirr and Q_stirr are functions of the planetesimal eccentricities and inclinations (for more details see Chambers 2006). However, this is a local approach. The gravitational excitation decreases with increasing distance between the planetary embryo and the planetesimals. Hasegawa & Nakasawa (1990) showed that when the distance from the planetary embryo is larger than about four times its Hill radius, the excitation over the planetesimals decays significantly. Therefore, we need to restrict this effect to the neighborhood of the planetary embryo. Using the EVORB code (Fernandez et al. 2002), we fit a a modulation function to reproduce the excitation over the quadratic mean value of the eccentricity of a planetesimal. We found that this excitation is well reproduced by (7)where Δ = R − R_P represents the distance from the planet (R_P is the planet radius orbit), R_H is the planetary embryo Hill radius, and f(Δ) guarantees that the eccentricity and inclination profiles of the planetesimals are smooth enough for a numerical treatment and that the planetary excitation on planetesimals is restricted to the embryo neighborhood.

On the other hand, the eccentricities and inclinations of the planetesimals are damped by the gaseous component of the protoplanetary disk. This damping depends on the planetesimal relative velocity with respect the gas, , and on the ratio between planetesimal radius and the molecular mean free path, λ. Adopting a gaseous disk mainly composed of molecular hydrogen (H₂), the last is given by (Adachi 1976) (8)where μ_H2 and d_H2 are the molecular weight and molecular diameter of the molecular hydrogen, and ρ_g is the volumetric gas density.

As in the recent work of Fortier et al. (2013), we considered three different regimes (Rafikov 2004; Chambers 2008)

• the Epstein regime: r_p<λ_H2;

• the Stokes regime: r_p>λ_H2 and Re<Re_trans; and

• the quadratic regime: r_p>λ_H2 and Re>Re_trans;

where is the Reynolds number and Re_trans = 20 is the transition between the Stokes and quadratic regimes (Rafikov 2004). The viscosity ν corresponds to the molecular viscosity, given by (9)where c_s is the local speed of sound.

Incorporating the different regimes is important because smaller fragments (products of planetesimal fragmentation) might be present in the Stokes or Epstein regimes. The three drag regimes can be characterized in terms of the stopping time given by (Chambers 2008) where ρ_p is the planetesimal density. The relative velocity between planetesimals and the gas is given by (13)where v_k is the Keplerian velocity and η = (v_k − v_g) /v_k is the ratio of the gas velocity to the Keplerian velocity.

The damping rates of the eccentricities and inclinations for each regime are given by (Rafikov 2004; Chambers 2008) where , with , and Finally, the evolution of the eccentricities and inclinations are given by solving the coupled equations by The gas drag also causes an inward planetary orbit migration. Then, the change rate of the major semi-axis is given by (22)

2.3. Oligarchic accretion regime

As in previous works, we considered that the embryos grow in the oligarchic regime. Assuming the particle-in-a-box approximation, the planetesimal accretion rate of the j species is given by (Inaba et al. 2001) (23)where M_C is the embryo core mass and P_coll is the collision probability between the planetesimal j species and the embryo (see Guilera et al. 2010, for the explicit expression of P_coll). The collision probability is a function of the embryo core radius (R_C), the embryo Hill radius, and the relative velocity between the planetesimal j species and the embryo, which is given by (24)When the embryo has a substantial envelope we have to incorporate the enhancement of capture cross-section of the embryo. As in previous works, we used the prescription given by Inaba & Ikoma (2003) to calculate the embryo enhanced radius , where they proposed replace for R_C in the expressions of the collision probability.

The feeding zone of the embryo is often defined as the ring around it where planetesimals can be accreted. We defined the width of the feeding zone as about four times the embryo Hill radius (at both sides of the embryo ). To do this, we integrate Eq. (23) over the radial grid (25)where ψ(R,R_P,R_H) is a normalization function that satisfies . In contrast with our previous work, we chose that (26)where Γ is the gamma function. With this new choice of ψ, , so the tail of the function has a negligible contribution in Eq. (25) and continues be smooth for a numerical treatment. We employed a Simpson rule to integrate Eq. (25) where at least ten radial bins between R_P − 4R_H and R_P + 4R_H are considered. Finally, the total planetesimal accretion rate is given by (27)The rest of the model is the same as described in detail in Guilera et al. (2010).

3. Planetesimal fragmentation

We incorporated a model of planetesimal fragmentation in our model of giant planet formation that is based on the BOULDER code (Morbidelli et al. 2009, and supplementary material). This code models the accretion and fragmentation of a population of planetesimals (because our model of giant planet formation starts with an embryo already formed, sourrunded by a swarm of planetesimals, i.e., in th oligarchic growth regime, we first took into account only the corresponding fragmentation prescriptions, see next sections).

According to this model, if is the specific impact energy per unit target mass (energy required to disperse 50% of the target mass) and Q is the collisional energy per unit target mass, the collision between a target of mass M_T and a projectile of mass M_P (with M_P ≤ M_T) gives a remnant of mass M_R that is given by (28)If M_R>M_T, the collision results in accretion. On the other hand, if M_R< 0, the target is fully pulverized and its mass is lost. In general, is function of the radius of the target. However, because the model considers that the mass of the remnant is a function of (M_T + M_P), for consistency, must be calculated using an effective radius which is given by (29)where ρ is the density of the planetesimals. The mass ejected from the collision, defined as (M_T + M_P − M_R), is distributed following a power-law mass distribution dn/ dm ∝ m^−p between the minimum bin mass considered and the bin mass corresponding to the larger fragment M_F given by (30)We found that for some hugely catastrophic collisions (when M_R ≪ M_T + M_P) the mass of the larger fragment is greater than the mass of the remnant. Therefore, we set M_F = 0.5M_R for these collisions.

The exponent p of the mass distribution is given by (31)where q is the exponent of the cumulative power-law distribution, and is given by (32)For the specific impact energy, we adopted the prescription given by Benz & Asphaug (1999). We used the prescription for basalts at impact velocities of 5 km s^-1 given by (33)using a planetesimal density of ρ_p = 1.5 g cm^-3. In this prescription we used the effective radius (given by Eq. (29)) instead of the planetesimal radius.

In our global model we consider the evolution by migration, accretion and fragmentation of a population of planetesimals of radii between 1 cm and (a free parameter). However, in a collisional regime the evolution of the population is ultimately governed by the size distribution of the smallest objects. This means that the truncation in r_p = 1 cm can generate the accumulation of spurious mass in the smaller fragments. To avoid this problem, we extrapolated the size distribution to a minimum fragment size of r_p = 0.01 cm, when we calculated the fragmentation process. In this way, the mass ejected from the collision is distributed between the mass bin corresponding to r_p = 0.01 cm and the bin mass corresponding to M_F, as mentioned above. This means that we are considering that the mass distributed below the mass bin corresponding to r_p = 1 cm is lost. Moreover, if the mass bin corresponding to M_R is below the mass bin corresponding to r_p = 1 cm, we considered that the target is fully pulverized.

The total number of collisions between targets j and projectiles i in a time Δt is given by (Morbidelli et al. 2009, and supplementary material) (34)where is the intrinsic collision probability, n_pj (r_pj) and n_pi (r_pi) are the numbers (radii) of targets and projectiles, respectively, and is the gravitational focusing factor. The time step Δ t is limited by a physical condition. Using that is the total number of collisions of the targets j given by (35)and defining τ_i as (36)we can write, (37)Then, cannot be greater than n_pj, so for our model we adopted that Δ t < 0.1 /τ_i. This condition implies that for our global model the time step cannot be greater than 10^-4 Myr for km, and 10^-5 Myr for km.

3.1. Implementation of the fragmentation model

The evolution of the surface densities of planetesimals obeys a continuity equation, (38)where R and r_p reference radial and planetesimal size dependencies, and ℱ are the sink terms. In our previous works, we only considered the planetesimal accretion by forming embryos as a sink term. With the incorporation of the planetesimal fragmentation, we introduce a new sink term in Eq. (38), which is solved with a full implicit method in finite differences.

To incorporate the fragmentation process in the global model, we defined a fragmentation zone for each embryo where we calculated the collisional process (Fig. 1). This zone extends to eight Hill radii on either side of the embryo (twice the feeding zone)¹. In this way we reduced the computational cost and safely guaranteed that collisions that produce fragmentation (craterization or catastrophic collisions) are within this zone, since they do not extend far away from the feeding zone. In fact, the excitations of the eccentricities and inclinations of the planetesimals (and hence the relative velocities) abruptly decay far away from the feeding zone (Eq. (21)).

Fig. 1 Schematic illustration of the fragmentation zone for an isolated planet in the radial grid. IEFZ (OEFZ) represents the inner (outer) edge of the fragmentation zone, while RP is the radial bin corresponding to the planet. The fragmentation zone extends to eight Hill radii on either side of the planet.

In each radial bin of the fragmentation zone, the eccentricities, inclinations, and surface densities for each planetesimal size are defined. As we mentioned above, to calculate the evolution of the system mass, we considered a size distribution between 1 cm and a , where initially the total mass of the system is distributed in planetesimals of radii . However, to model the fragmentation we extrapolated the planetesimal sizes (also the eccentricities, inclinations, and numbers of bodies) to planetesimals of radius r_p = 0.01 cm. In this way, we avoided spurious mass accumulation in the smaller planetesimals of the distribution.

The implementation methodology is as follows. The eccentricities, inclinations, and surface densities are defined for each size bin (between r_p = 1 cm and ) and for each radial bin in the fragmentation zone. From the surface densities, the number of bodies for each size bin and for each radial bin are calculated. With these data, we extrapolate the inclinations, eccentricities, and number of bodies for the corresponding size bins between r_p = 0.01 cm and r_p = 1 cm. Then, we take a target j belonging to the radial bin IEFZ (inner edge of the fragmentation zone, Fig. 1). We take a projectile i (with r_pi ≤ r_pj, i.e. from r_pi = 0.01 cm to r_pi = r_pj) of the radial bin IEFZ, and calculate if the orbits of the target and the projectile overlap (taking into account only the eccentricities of the target and projectile). If the orbits overlap, the number of collisions between targets j and projectiles i are calculated. With this information we can calculate how much mass the projectiles i disperse from targets j and how this mass (remnant plus fragments) is distributed in smaller planetesimals than target j². This is repeated for all projectiles belonging to the radial bin IEFZ. Then we move to the radial bin IEFZ+1 and repeat the process for all the projectiles that correspond to target j of the radial bin IEFZ. This process is repeated until the radial bin OEFZ (outer edge of the fragmentation zone) is achieved. This process is repeated for all targets j from radial bin IEFZ, that is, from r_pj = 1 cm to . Then we move to the radial bin IEFZ+1 and repeat the process for all the targets. The process is repeated until it reaches the radial bin OEFZ for the targets.

When the process is completed, we have the change in mass (loss and gain) for each planetesimal size and for each radial bin, the product of the planetesimal fragmentation. With this we can calculate the change in the surface densities for each planetesimal size inside the fragmentation zone.

4. Results

The protoplanetary disk is defined between 0.4 au and 20 au, using 2500 radial bins logarithmically equally spaced. We applied our model to study the role of planetesimal fragmentation on giant planet formation. We calculated the in situ formation of an embryo located at 5 au for different values of the minimum mass solar nebula (MMSN; Hayashi 1981), which is given by where Σ_p and Σ_g represent the planetesimal and gaseous surface densities, T is the temperature profile, and ρ_g is the volumetric density of gas at the mid-plane of the disk. The discontinuity at 2.7 au in the surface density of planetesimals is caused by the condensation of volatiles, often called snow line. For numerical reasons, and following Thommes et al. (2003), we spread the snow line with a smooth function, so the surface density of planetesimals is described by (43)We carried out two different sets of simulations. In the first one, we took into account that the planetesimal surface density evolved only by accretion of planetesimals by the embryo and for the orbital migration of planetesimals. In the other set of simulations the planetesimal surface density evolved by accretion, orbital migration, and planetesimal fragmentation.

For the gaseous component of the disk, Alexander et al. (2006) found that after a few Myr of viscous evolution the disk could be completely dissipated by photo-evaporation in a time-scale of 10⁵ yr. For simplicity, we considered that the gaseous component of the disk exponentially dissipated in 6 Myr, when photo-evaporation acts and completely dissipates it. So, we ran our models until the embryo achieved the critical mass (when the mass of the envelope equals the mass of the core) or for 6 Myr.

We started our simulations with an embryo of 0.005 M_⊕, which has an initial envelope of ~10^-13M_⊕, immersed in an initial homogeneous population of planetesimals of radius (we considered different values for : 0.1, 1, 10, 100 km). It is important to remark that our initial conditions correspond to the beginning of the oligarchic growth. Employing statistical simulations, Ormel et al. (2010) found that starting with an homogeneous population of planetesimals of radius r₀, the transition from the runaway growth to the oligarchic growth is characterized by a power-law mass distribution given by dn/ dm ∝ m^−p (with p ~ 2.5), between r₀ and a transition radius r_trans for the population of planetesimals, and isolated bodies (planetary embryos). This implies that most of the solid mass lies in small planetesimals. For simplicity, and considering the fact that most of the mass lies in the smaller planetesimals of the population, we used a single-size distribution instead of a planetesimal-size distribution to represent the initial planetesimal population. Consequently, our initial conditions are consistent with the oligarchic growth regime using as r₀.

In this work we aim to analyze how planetesimal fragmentation impacts on the process of planetary formation. Accretion collisions between planetesimals are important for studying the transition from planetesimal runaway growth to the oligarchic growth. In the oligarchic growth, embryos gravitationally dominate the dynamical evolution of the surrounding planetesimals. As embryos grow, they increase the planetesimal relative velocities and collisions between planetesimals result in fragmentation (erosive or disruptive collisions). For these reasons, we focused our analysis on fragmentation collisions. However, as we show in next sections, for some special cases the total planetesimal accretion rates are dominated by the accretion of very small fragments (r_p ~ 1 m). For these small fragments, collisions between them (and obviously with smaller fragments) result in accretion. Therefore, we also calculated coagulation between planetesimals for these special cases (see next sections for a detailed discussion of this topic).

As we mentioned in the previous section, when planetesimal fragmentation is considered (when it is not considered, we used a single-size distribution to represent the planetesimal population) we used a discrete size grid between 1 cm and to represent the continuous planetesimal size distribution where the fragments, products of planetesimal fragmentation, are distributed. The size step is given by Eq. (1). Because we used the same size step independently of the value of , this implies that we used 21 size bins for km, 26 size bins for km, 31 size bins for km, and 36 size bins for km. In all cases, the initial total mass of solids along the disk remains in the size bins corresponding to . The collisional evolution of planetesimals is the mechanism that regulates the exchange of mass between the different size bins.

In previous works (Guilera et al. 2010, 2011) we showed that the in situ simultaneous formation of solar system giant planets occurred on a time-scale compatible with observed estimates only if most of the solid mass accreted by the planets remains in small planetesimals (r_p < 1 km). Recently, Fortier et al. (2013) studied the role of planetesimal sizes in planetary population synthesis in detail. They found that including oligarchic growth makes the formation of giant planets using big planetesimals (r_p ~ 100 km) unlikely, and only if most of the mass of the system remains in small objects (r_p ~ 0.1 km) cores grow enough to form giant planets. However, small planetesimals have lower specific impact energies per unit target mass. This means that small planetesimals suffer catastrophic collisions at lower collisional energies.

The results of the two sets of simulations are summarized in Table 1. For smaller planetesimals (r_p = 0.1 km), the collisional evolution completely inhibited planetary formation. When planetesimal fragmentation was not taken into account, the embryo was able to achieve the critical mass for 2−10 values of the MMSN before the dissipation of the nebula³. Moreover, for disks more massive than four times the MMSN the times at which embryo achieves the critical mass are very short (≲0.5 Myr), while the critical masses are high (≳30 M_⊕). When we included planetesimal fragmentation, the picture drastically changed. For none of the analyzed cases the critical mass was reached. Moreover, even for the case of the more massive disk (10 MMSN), the planetary embryo was unable to achieve one Earth mass. In Fig. 2, we show the total accretion rate of planetesimals as a function of the core mass for km. When planetesimal fragmentation is incorporated, the total planetesimal accretion rate is completely inhibited despite the low values of the core masses. The small differences between the disks at which the values of the core mass planetesimal fragmentation start to decrease the accretion rates are caused by the amount of gas of the different disks. The more massive the disk, the higher the damping in the planetesimal eccentricities and inclinations, hence the planetesimal relative velocities are lower for the same value of the core mass.

Fig. 2 Total planetesimal accretion rates for three different disks as a function of the core mass. The case without planetesimal fragmentation is shown as dotted lines. For this case, the critical mass is achieved for the three different disks in less than 6 Myr. When planetesimal fragmentation is incorporated the total planetesimal accretion rates significantly drop. Despite the lower values of the core masses, the fragmentation process truncated the embryo growth.

These drops in the planetesimal accretion rates are due to a drastic diminution in the mean value of the surface density of planetesimals of radius 0.1 km in the feeding zone compared with the case without planetesimal fragmentation.

Fig. 3 Time evolution of the mean value of the surface density of planetesimals in the feeding zone for a disk 6 times more massive than the MMSN. For the case without planetesimal fragmentation (red solid curve), the total surface density corresponds to planetesimals of radius 0.1 km. This is not the case when planetesimal fragmentation is considered (solid black curve). However, the total surface density is dominated by the surface density of planetesimals of 0.1 km (dashed red curve). (Color version online.)

In Fig. 3, we plot the time evolution of the mean value of the surface density of planetesimals in the feeding zone for a disk six times more massive than the MMSN. When we included planetesimal fragmentation, the total planetesimal surface density significantly drops. It is important to remark that for the case without planetesimal fragmentation, the total planetesimal surface density corresponds to the surface density of planetesimals of radius 0.1 km. This is not the case when we considered planetesimal fragmentation. However, as we can see in Fig. 3, the surface density of planetesimals of radius 0.1 km is almost the total planetesimal surface density. This is because of two effects. First, most of the mass distributed in the fragments, products of the collisions, is lost below the size bin corresponding to 1 cm. As we mentioned above, we considered that the mass loss in collisions is distributed following a power-law mass distribution between fragments of r_p = 0.01 cm and the largest fragment of mass M_F. Accordingly, if the exponent p of the mass distribution (dn/ dm ∝ m^{− p}) is higher than 2, most of the mass is distributed in the smaller fragments, and if p < 2, most of the mass is distributed in bigger fragments. When p > 2, most of the mass is distributed in fragments smaller than 1 cm, so this mass is lost in our model. In Fig. 4, we show the values of q, the exponent of the cumulative power-law distribution (Eq. (32)), as a function of . We can see that except for values between , q is always lower than −3. This means that in general q < − 3, so p > 2 (see Eq. (31)) and most of the mass distributed in fragments is lost. It is important to remark that for these step distributions, the integration over the mass, between m = 0 and m = M_F, diverges. In the Boulder code, this problem is solved by using that these mass distributions are valid between m = m_t (m_t represents a transition mass) and m = M_F and using an ad hoc mass distribution with an exponent p = 11/6 between m = 0 and m = m_t. In this approach, most of the mass loss in collisions is distributed in the radial bin correspondig to r_t (the transition radius corresponding to m_t). However, we did not follow this approach. We truncated the mass distribution in a minimum mass, the one corresponding to r_p = 0.01 cm (we tested lower values than r_p = 0.01 cm and found analogous results). This is an alternative approach, adopting a fixed value of r_t = 0.01 cm and calculating the exponent p, such that the integration over the mass between m = 0 and m = M_F, converges⁴. In this way, we always guarantee for each collision that the mass distributed between r_p = 0.01 cm and the radial bin corresponding to planetesimals of mass M_F does not exceed the mass ejected from the collision (M_T + M_P − M_R).

Fig. 4 Exponent of the cumulative power-law distribution as a function of given by the fragmentation model. The value q = −3 implies that p = 2 (Eq. (31)). This means that the mass lost in a collision is distributed homogeneously among the fragments.

The second effect is that the remnants, products of the collisions between planetesimals of radius 0.1 km, are quickly pulverized. Because these planetesimals initially contain all the mass of the system, the pulverization of the remnants implies a high loss of mass. For example, in Fig. 5 we can see the radial profiles of the eccentricities and inclinations at different times at the embryo’s neighborhood. The gravitational perturbations of the planet increase the eccentricities and inclinations of the planetesimals near the planet’s location. From this profile we can analyze the relative velocities, and the ratio , when targets and projectiles belong to the same radial bin. Figure 6 shows the time evolution of the radial profiles for the relative velocities (top panel) and for the ratio (bottom panel) when targets and projectiles belong to the same radial bin. It follows from Eq. (28) that if the ratio corresponding to a collision is higher than ~2.5, the remnant of such a collision is pulverized. As we can see from Fig. 6, the collisions between planetesimals with radius of 0.1 km quickly become supercatastrophic and the mass of the remnants, their products, is lost.

Fig. 5 Time evolution of the radial profiles of the eccentricities (solid lines) and inclinations (dashed lines) in the embryo’s neighborhood for planetesimals with a radius of 0.1 km. With time, the planet’s gravitational excitation increases the eccentricities and inclinations of the planetesimals. It is also clear that planetesimals do not reach equilibrium values (β = i/e ≠ 0.5). (Color version online.)

Fig. 6 Time evolution of the radial profiles for the relative velocities (top) and for the ratio (bottom) for targets and projectiles of radius 0.1 km when both belong to the same radial bin. In contrast, the relative velocities are always lower than 1 km s-1, the ratio is much higher than unity at very early times. For this reason, collisions quickly become highly catastrophic at early times. (Color version online.)

Finally, in Fig. 7 we show the time evolution of the number of planetesimals and the planetesimal surface densities at the embryo’s radial bin. The number of planetesimals with a radius of 0.1 km is quickly reduced by the collisional evolution. Moreover, the generation of fragments does not compensate the loss of planetesimals with a radius of 0.1 km. In fact, the values of the planetesimal surface densities for planetesimals smaller than 0.1 km are always ≪1 gr cm^-2.

Fig. 7 Time evolution of the number of planetesimals (top) and planetesimal surface density (bottom) as a function of planetesimal radius at the embryo’s radial bin. The number of planetesimals evolves by accretion of the embryo, migration, and planetesimal fragmentation. The generation of fragments (produced by the collisional evolution) does not compensate for the large loss of planetesimals with a radius of 0.1 km. (Color version online.)

We found similar results for the others values of . Collisions between planetesimals of radius become supercatastrophic and significantly reduce the total planetesimal accretion rates. This effect, combined with the fact that most of fragment mass is deposited in size bins smaller than the one corresponding to r_p = 1 cm, prevented the formation of a core from reaching the critical mass.

Only for km and massive disks, we were able to form cores with masses greater than one Earth mass. This is because despite of higher relative velocities for bigger planetesimals, the ratio is lower compared to small planetesimals (Fig. 8). In Fig. 9, we plot the planetesimal accretion rates for km and a disk ten times more massive than the MMSN. With the red solid line, we plot the total planetesimal accretion rate. This accretion rate abruptly drops compared with the case without planetesimal fragmentation because of the lower accretion of planetesimals with a radius of 100 km. The accretion of fragments also does not compensate for the lower accretion of planetesimals with a radius of 100 km.

Finally, in Fig. 10 we show the time evolution of the radial surface density profiles of planetesimals with a radius of 100 km, for the case with and without planetesimal fragmentation. The profiles are the same at 0.5 Myr. At 1 Myr the planetesimal surface density around the planet’s location (5 au) is clearly lower for the case with planetesimal fragmentation. This is caused by planetesimal fragmentation, but not by planetesimal accretion by the embryo. We can see from Fig. 9 that at this time the accretion rate of planetesimals with a radius of 100 km is practically the same as that without planetesimal fragmentation. Eventually, the lower planetesimal surface density around the planet’s location increases even more, for the case with planetesimal fragmentation. Finally, at 4 Myr, the planetesimal surface density is almost zero around the planet’s location. We also can see that the mass-loss through planetesimal fragmentation is higher near the planet’s location and lower far away from the planet.

Fig. 8 Time evolution of the relative velocity profiles (top) and the ratio profiles (bottom) for targets and projectiles of radii 100 km. The profiles represent the case for km and a disk six times more massive than the MMSN. Despite the higher values of the relative velocities of planetesimals with rp = 100 km, the ratio is lower compared with the case for small planetesimals. This is because the specific impact energy () for planetesimals of rp = 100 km is approximately three orders of magnitude higher than that for planetesimals of rp = 0.1 km. (Color version online.)

Fig. 9 Planetesimal accretion rates as a function of time for an embryo located at 5 au in a disk ten times more massive than the MMSN, where initially all the solid mass is deposited in planetesimals of 100 km of radius. The black dashed line corresponds to the case without planetesimal fragmentation. The red solid line corresponds to the total planetesimal accretion rate with planetesimal fragmentation. The accretion rate of planetesimals with a radius of 100 km is plotted as a green solid line. We also plot the accretion rates for other planetesimal sizes. (Color version online.)

Fig. 10 Time evolution of the radial surface density profiles of planetesimals with a radius of 100 km for a disk ten times more massive than the MMSN. Solid lines represent the case with planetesimal fragmentation, dashed lines represent the case without it. For this last case, the planet achieves the critical mass at 4.07 Myr. When planetesimal fragmentation is considered, the planet does not achieve the critical mass. (Color version online.)

We found that our results are insensitive to the numbers of radial and size bins. We tested our results with twice radial and size bins and obtained analogous results.

4.1. Fragment distribution

As we mentioned in previous sections, most of the fragment mass produced by the collisions between planetesimals is deposited in size bins smaller than 1 cm, because the exponent p of the power-law of the fragment mass distribution is generally larger than 2. However, other works suggested that the exponent p should be in a range between 1 and 2. This implies that most of the fragment mass is deposited in larger fragments and the mass loss is much smaller than in our model. Kobayashi & Tanaka (2010) developed a similar fragmentation model (in a qualitative way about the outcome of a collision) where the fragment mass is also distributed following a power-law distribution (dn/ dm ∝ m^{− p}). In Kobayashi et al. (2010, 2011) and Ormel & Kobayashi (2012), this model was applied to study planetary formation adopting a value of p = 5/3. These authors found that in general planetesimal fragmentation inhibits the formation of massive cores. The formation of cores with masses greater than 10 M_⊕ is only possible for massive disks, and only if large planetesimals are considered (r_p ≥ 100 km).

We applied our model using a fixed exponent p = 5/3 for the power-law mass distribution of the fragments. We again ran the simulation for a disk ten times more massive than the MMSN using an initial population of planetesimals with a radius of 100 km (our best case). Using a fixed exponent p = 5/3, we found that planetesimal fragmentation favors the formation of a massive core. For this simulation, we found that the embryo achieved the critical mass at 3.61 Myr (~0.5 Myr early than without planetesimal fragmentation) with a core of 18.58 M_⊕ (~6.5 M_⊕ smaller than without planetesimal fragmentation).

In Fig. 11 we plot the time evolution of the total planetesimal accretion rates for the case without planetesimal fragmentation and p = 5/3, and the planetesimal accretion rates for different planetesimal sizes for p = 5 /3. Without planetesimal fragmentation, the planetesimal accretion rate of planetesimals with a radius of 100 km corresponds to the total planetesimal accretion rate. However, this is not the case when planetesimal fragmentation is considered. Between ~0.5 Myr and ~1 Myr the accretion rate of planetesimals with r_p = 100 km is slightly lower than that without fragmentation. However, the total planetesimal accretion rate is higher. This is because of the accretion of fragments, especially for the accretion of fragments between ~0.1 km and ~25 km. Then, the accretion rate of planetesimals of r_p = 100 km increases at ~1 Myr. This is because at this time the ratio becomes higher than unity, so the enhanced in the capture-cross section due to the embryo’s envelope for planetesimals with a radius of 100 km make the accretion of such planetesimals more efficient. The total planetesimal accretion rate remains higher than without fragmentation until ~2.5 Myr. This excess in the total planetesimal accretion rate is the reason that at ~2.5 Myr the planet has a core of ~12.5 M_⊕ (Fig. 13 top). At the same time, the planet without planetesimal fragmentation has a core of ~4.5 M_⊕. After 2.5 Myr, the total planetesimal accretion rate decreases because of planetesimal fragmentation until the planet reaches the gaseous runaway phase.

Fig. 11 Planetesimal accretion rates as a function of time for an embryo located at 5 au in a disk ten times more massive than the MMSN, where initially all the solid mass is deposited in planetesimals with a radius of 100 km. The dashed black curve represents the total planetesimal accretion rate for the case without planetesimal fragmentation. The red solid curve represents the total planetesimal accretion rate when p = 5/3 is a fixed value for all collisions. We also plot the planetesimal accretion rates for the different planetesimal sizes for p = 5/3. (Color version online.)

Fig. 12 Time evolution of the number of planetesimals adopting km and a disk ten times more massive than the MMSN. Filled circles represent the case where the mass loss in collisions is deposited in smaller fragments, while open circles represent the case where the mass loss in collisions is deposited in larger fragments. For this last case, the planet achieved the critical mass at ~3.61 Myr. (Color version online.)

In Fig. 12, we compare the time evolution of the number of planetesimals at the planet’s radial bin for the cases where the exponent of the power-law that represents the planetesimal mass distribution has a constant value of 5 / 3 and where this exponent is calculated via Eq. (31). There are many more fragments between ~0.1 km and ~25 km in the first case for 0.75 Myr, 1 Myr, and 2 Myr. This is because for this case, the mass loss in collisions is distributed in larger fragments. The accretion of these fragments favors the formation of a massive core. Then, for 3 Myr, there are fewer fragments between ~0.1 km and ~25 km than when the mass loss in collisions is distributed in smaller fragments. But this is because of the accretion of such fragments. Finally, for the first case the formation of the core occurs at ~3.61 Myr.

A similar behavior occurs for less massive disks. However, for a disk eight times more massive than the MMSN the planet did not reach the critical mass. After 6 Myr of evolution, the planet achieved a total mass of ~20 M_⊕ (~14 M_⊕ for the core and ~6 M_⊕ for the envelope, Fig. 13 bottom). Although this final core is smaller than the case without planetesimal fragmentation (Table 1), at 4 Myr the planet has a core slightly larger than 10 M_⊕ (in the case without planetesimal fragmentation the planet has a core of ~4 M_⊕ at the same time).

Fig. 13 Core masses (solid lines) and envelope masses (dashed lines) as a function of time for the cases with planetesimal fragmentation considering p = 5/3 (black lines) and without planetesimal fragmentation (red lines), for a disk ten times more massive than the MMSN (top) and a disk eight times more massive than the MMSN (bottom). For this last case, the planets did not achieve the critical mass after 6 Myr of evolution. (Color version online.)

Fig. 14 Top: planetesimal accretion rates for different planetesimal sizes, for a disk ten times more massive then the MMSN, in which p = 5/3 and km, as a function of time. The dashed black line represents the total planetesimal accretion rate for the case without planetesimal fragmentation (WF). The dashed red line represents the total planetesimal accretion rate for the case with planetesimal fragmentation in which p is calculated with Eq. (31), our basis model (BM). The solid red line represents the total planetesimal accretion rate for the case with planetesimal fragmentation with p = 5/3. Bottom: time evolution of the total planetesimal accretion rates (solid lines) and the accretion rates of fragments of 1 m (dashed lines) for the case where only planetesimal fragmentation is considered (red lines) and for the case in which planetesimal coagulation and fragmentation are considered (black lines). The plot corresponds to a disk ten times more massive then the MMSN and for the case of km. (Color version online.)

We also found that for small planetesimals, if most of the fragment mass is deposited in larger fragments the total planetesimal accretion rate becomes higher than in the case in which most of the fragment mass is deposited in smaller fragments. However, for these cases the total planetesimal accretion is always lower than for the case in without planetesimal fragmentation. We calculated again the simulations for a disk ten times more massive than the MMSN. However, for km, the situation is different. For these cases, the accretion rate corresponding to quickly drops and the total accretion rate is dominated by fragments of ~r_p = 1 m. However, collisions between these small fragments (and obviously with smaller fragments) are not disruptives but rather the outcome of a collision that results in an effective accretion. This means that in this scenario planetesimal coagulation is necessary. As an example, in Fig. 14 (top) we plot the planetesimal accretion rates for km. The accretion rate of planetesimals with a radius of radius 0.1 km drops significantly (green curve), and the total accretion rate is ultimately dominated by small planetesimals (r_p ~ 1 m, violet curve). The planetesimal accretion rates of fragments with ~r_p = 1 m become significant high values of the total planetesimal accretion rates, but this effect is fictitious because small planetesimals coagulate and form larger bodies. For this case, we stopped the simulation at 0.75 Myr because at this time the core had reached a mass of ~12 M_⊕. But again, these results might be fictitious. A planetesimal coagulation model is necessary in these cases.

This effect does not occur for large planetesimals, where the accretion of such small fragments does not significantly contribute to the total accretion rate. In fact, for large planetesimals, the total accretion rate is always dominated by the accretion of planetesimals of radius and the fragments that significantly contribute to the total accretion rate have radii larger than 0.1 km (collisions are disruptives for such fragments).

However, to be more rigorous with these intuitive analyse, we recalculated the simulations for km and km, but also incorporating coagulation between planetesimals in the model. Following the Boulder code, we considered that the outcome of a collision results in accretion if the mass of the remnant is larger than the mass of the target. We considered that when a coagulation between planetesimals occurs, the target and projectile are removed from their corresponding radial bins, and a new object of mass M_T + M_P is put in the radial bin corresponding to the target, that is, we considered perfect accretion. We pursued a very simple intention with this: we wish to analyze whether the coagulation between planetesimals modifies (or not) the results found in this section. Therefore, we considered for simplicity the size grid that represents the continuous planetesimal size distribution to be fixed. It is important to note that the computational costs are much higher.

For km, we found identical results for both cases, considering only planetesimal fragmentation or considering planetesimal coagulation and fragmentation. As we argued before, this is because small fragments contribute negligibly to the total planetesimal accretion rate (Fig. 11).

However, for km, we found that the incorporation of planetesimal coagulation significantly modified the results. In Fig. 14 (bottom), we plot the time evolution of the total planetesimal accretion rate and the accretion rate of fragments of 1 m for the case with only planetesimal fragmentation and for the case with planetesimal coagulation and fragmentation. For this last case, after 6 Myr of evolution the core achieved a mass of only ~3 M_⊕. In both cases, the accretion of fragments of ~1 m governed the total accretion rates, but the incorporation of planetesimal coagulation drastically lowers the accretion of such fragments. It is clear that when the accretion of very small fragments becomes important in the total planetesimal accretion rate, a full planetesimal collisional model (which includes coagulation and fragmentation) is needed. We will study this topic in detail by developing a full planetesimal collisional model in a future work, incorporating an adaptative size grid to study in more detail the collisional evolution of the planetesimal population and starting from the planetesimal runaway growth.

Finally, planetesimal coagulation did not quantitatively modify the results shown in Table 1 for low values of . For these cases, all the accretion rates for the different planetesimal sizes significantly decrease with time because the fact that most of the mass loss in collisions is distributed below the size bin corresponding to r_p = 1 cm.

Fig. 15 Factor as a function of the core mass for different planetesimal sizes. The plot corresponds to the case of a disk ten times more masive than the MMSN and km when planetesimal coagulation and fragmentation are considered. The black dashed line corresponds to the evolution of the term (RP/RH)1 / 2, where RP is the radius of the planet. (Color version online.)

4.2. Accretion of small fragments

We briefly discuss the accretion of small fragments (pebbles). Lambrechts & Johansen (2012) found that pebbles are accreted extremely efficiently by embryos. Pebbles with the appropriate Stoke number have a capture cross-section as large as the Hill radius of the embryo, even if the embryo does not have a gaseous envelope. When they pass within the Hill radius, they spiral down to the embryo’s physical radius because of gas drag. Lambrechts & Johansen (2012) found that the pebble accretion rates (in the Hill accretion regime) is given by (44)where v_H = Ω_kR_H, and Ω_k is the Keplerian frecuency. As the authors note, in the classical scenario of planetesimal accretion, planets do not accrete planetesimals from the Hill radius. Instead, planets accrete planetesimals from a fraction α^1/2 of the Hill radius, with α = R_C/R_H.

We used the planetesimal accretion rates of Inaba et al. (2001) given by (45)In the low velocity regime, the probability collision or small planetesimals is given by (46)because we also considered the enhanced radius due to the planet’s gaseous envelope . In terms of the pebble accretion rate, our planetesimal accretion rate (for small fragments) is given by (47)In Fig. 15 we plot the evolution of the term as a function of the core mass for different planetesimal sizes for the case of a disk ten times more masive than the MMSN and km when planetesimal coagulation and fragmentation are considered. For small fragments (r_p ≤ 1 m) the evolution is almost the same. This is because for these small fragments the enhanced capture radius becomes the planet’s radius⁵ for low values of the core mass. Due to the factor 5.65 in Eq. (47), our planetesimal accretion rates for r_p ≲ 1 m are higher than the pebble accretion rates of Lambrechts & Johansen (2012) when the planet’s core mass becomes larger than ~0.2 M_⊕. Despite these high accretion rates, these small fragments have a neglegible contribution in our models because of the low values for the corresponding surface densities (see for example Fig. 7). It is important to note that the probability collision for smaller fragments (r_p ≲ 1 m) always corresponds to the low-velocity regime.

Fig. 16 Top: time evolution of the radial profiles of the surface density of planetesimals of rp = 100 km when planetesimal fragmentation is considered with p = 5 / 3 and for a disk ten times more massive than the MMSN. Black dashed lines correspond to the isolated formation of a planet located at 5 au. Red solid lines correspond to the simultaneous formation of two planets located at 5 au and 6 au. The different lines correspond to different times: 0, 0.5, 1, 1.5, 2, and ~2.5 Myr. Eventually, the planetesimal surface density decreases in the neighborhood of both embryos. Bottom: same as top panel, but for fragments that most contribute into the total accretion rate, rp ~ 25 km. In this case, the profiles increase around 5 au for t = 0.5,1,1.5 Myr, but for t = 2, ~2.5 Myr profiles decrease because of accretion. The increments in the surface density around 6 au correspond to the profiles at times 1.5, 2, and ~2.5 Myr. (Color version online.)

Fig. 17 Core masses as a function of time for the isolated formation of a planet located at 5 au (black dashed line) and for the simultaneous formation of two planets located at 5 au (red solid line) and 6 au (blue solid line). Models correspond to km for a disk ten times more massive than the MMSN, and planetesimal fragmentation is considered with p = 5/3. The simultaneous formation stopped at ~2.5 Myr because the distance between planets became smaller than 3.5 mutual Hill radii. (Color version online.)

4.3. Simultaneous formation of two embryos

Finally, we analyzed the in situ simultaneous formation of two embryos. We aimed to study if the fragments generated by an outer embryo (which have an inward migration) favor the formation of an inner embryo. We only analyzed the case of a disk ten times more massive than the MMSN, where initially all the solid mass of the system is deposited in planetesimals of r_p = 100 km and for the case p = 5 / 3. We located the embryos at 5 au and 6 au. Both embryos initially have cores of 0.005 M_⊕ and envelopes of ~10^-13M_⊕. The simulation stopped at ~2.5 Myr. At this time, the embryo located at 5 au achieved a total mass of 13.20 M_⊕ (12.10 M_⊕ for the core and 1.10 M_⊕ for the envelope), while the embryo located at 6 au achieved a total mass of 2.61 M_⊕ (2.60 M_⊕ for the core and 0.01 M_⊕ for the envelope). The simulation stopped because the distance between the embryos became smaller than 3.5 mutual Hill radii. When two embryos are too close, their mutual gravitational perturbation may lead to encounters or collisions between them.

The time evolution of the radial profiles for the surface density of planetesimals of r_p = 100 km and r_p ~ 25 km (which are the two sizes that most contribute to the total planetesimal accretion rate, see Fig. 11, where the 25 km is the size of the fragments that contribute most) for the embryo located at 5 au are very similar (Fig. 16). Therefore, we do not expect a significant difference in the formation of this embryo. In Fig. 17 we plot the time evolution of the core mass for both embryos compared with the isolated embryo located at 5 au. We do not find a difference in the formation of the embryo located at 5 au between isolated and the simultaneous formation with an outer embryo, at least for the profile of the disk we analyzed here.

5. Conclusions

Employing the Boulder code, Morbidelli et al. (2009) found that the present size frequency distribution of bodies in the asteroid belt can be reproduced starting with an initial population of large planetesimals with sizes of between 100−1000 km. An initial population of large planetesimals is suggested by the works of Johansen et al. (2007), Cuzzi et al. (2008, 2010), and Youdin (2011). However, Fortier et al. (2009) demonstrated that the formation of massive cores able to achieve the critical mass to start the gaseous runaway phase, starting from large planetesimals, requires several million years, even for massive disks.

We studied the role of planetesimal fragmentation on giant planet formation. We developed a model for planetesimal fragmentation based on the Boulder code and incorporated it in our model of giant planet formation (Guilera et al. 2010, 2011), so that the planetesimal population evolved by planet accretion, migration, and fragmentation. We numerically studied how planetesimal fragmentation modified the formation of an embryo located at 5 au for a wide range of disk masses and planetesimal sizes.

We considered that initially all the solid mass of the system is deposited in planetesimals of radius , and that the mass loss in collisions is distributed by a power-law mass distribution between the largest fragment and the smallest size bin. The exponent of this power-law mass distribution was calculated by the model. We found that this exponent is generally larger than 2, so most of the mass is distributed in fragments smaller than 1 cm, which is the smallest size we considered for accretion. For this reason, most of the mass produced by the collisions was lost and fragments did not significantly contribute to the embryo growth.

We found that planetesimal fragmentation inhibits planet formation. Only for large planetesimals ( km), the mass of the embryo was larger than 1 M_⊕ for high-mass disks, but for all cases the core did not achieve the critical mass. However, if most of the mass loss in collisions was distributed in larger fragments, that is, when the exponent of the power-law mass distribution was lower than 2, we found that planetesimal fragmentation favored the relatively rapid formation of a massive core (larger than 10 M_⊕), as previously reported (Kobayashi et al. 2011; Ormel & Kobayashi 2012). The accretion of fragments between ~100 m and ~25 km increments the total planetesimal accretion rate, but this total accretion rate is always governed by the accretion of planetesimals of r_p = 100 km. For this case, we also analyzed whether the presence of an outer embryo modified the formation of an inner one. For this, we calculated the in situ simultaneous formation of two embryos located at 5 au and 6 au. In particular, we analyzed whether the fragment migration produced by the outer embryo affected the formation of the inner one. We did not find a difference between isolated and simultaneous formation for the embryo located at 5 au. But in this case, at ~2.5 Myr, the separation between the embryos became smaller than 3.5 mutual Hill radii, therefore the simulation was stopped. Possible mergers between massive embryos may lead to a rapid formation of massive cores.

On the other hand, Weidenschilling (2011) showed that the present size distribution observed in the asteroid belt can be also reproduced starting from planetesimals with a radius of ~0.1 km. Kenyon & Bromley (2012) concluded that the size distribution of the transneptunian objects can be reproduced starting from a massive disk composed of relative small planetesimals (r_p ≲ 10 km). However, for such small planetesimals, collisions quickly become highly catastrophic because of the low values for the specific impact energy. This means that targets are quickly pulverized and the total surface density of solids drops drastically. When the mass loss in collisions is distributed in smaller fragments, as the Boulder code predicts, planet formation is completely inhibited. When the mass loss in collisions is distributed in larger fragments, fragments of ~1 m ultimately govern the total planetesimal accretion rate. However, for these small fragments collisions result in accretion. Therefore, we repeated the simulations (for km) for the case

where the mass loss in collisions is deposited in larger fragments, but incorporated planetesimal coagulation. As we expected, planetesimal coagulation did not modify the process of planetary formation for km, because small fragments contribute only negligibly to the total planetesimal accretion rate. On the other hand, planetesimal coagulation significantly modified the results for km, inhibiting the formation of a massive core. We will analyze this case in detail, developing a full planetesimal collisional model in a future work.

Acknowledgments

We thank the anonymous referee for her or his comments, which helped us to improve the work. This work was supported by the PIP 112-200801-00712 grant from Consejo Nacional de Investigaciones Científicas y Técnicas, and by G114 grant from Universidad Nacional de La Plata.

References

All Tables

All Figures

	Fig. 1 Schematic illustration of the fragmentation zone for an isolated planet in the radial grid. IEFZ (OEFZ) represents the inner (outer) edge of the fragmentation zone, while RP is the radial bin corresponding to the planet. The fragmentation zone extends to eight Hill radii on either side of the planet.
In the text

Fig. 2 Total planetesimal accretion rates for three different disks as a function of the core mass. The case without planetesimal fragmentation is shown as dotted lines. For this case, the critical mass is achieved for the three different disks in less than 6 Myr. When planetesimal fragmentation is incorporated the total planetesimal accretion rates significantly drop. Despite the lower values of the core masses, the fragmentation process truncated the embryo growth.

In the text

Fig. 3 Time evolution of the mean value of the surface density of planetesimals in the feeding zone for a disk 6 times more massive than the MMSN. For the case without planetesimal fragmentation (red solid curve), the total surface density corresponds to planetesimals of radius 0.1 km. This is not the case when planetesimal fragmentation is considered (solid black curve). However, the total surface density is dominated by the surface density of planetesimals of 0.1 km (dashed red curve). (Color version online.)

In the text

	Fig. 4 Exponent of the cumulative power-law distribution as a function of given by the fragmentation model. The value q = −3 implies that p = 2 (Eq. (31)). This means that the mass lost in a collision is distributed homogeneously among the fragments.
In the text

Fig. 5 Time evolution of the radial profiles of the eccentricities (solid lines) and inclinations (dashed lines) in the embryo’s neighborhood for planetesimals with a radius of 0.1 km. With time, the planet’s gravitational excitation increases the eccentricities and inclinations of the planetesimals. It is also clear that planetesimals do not reach equilibrium values (β = i/e ≠ 0.5). (Color version online.)

In the text

Fig. 6 Time evolution of the radial profiles for the relative velocities (top) and for the ratio (bottom) for targets and projectiles of radius 0.1 km when both belong to the same radial bin. In contrast, the relative velocities are always lower than 1 km s-1, the ratio is much higher than unity at very early times. For this reason, collisions quickly become highly catastrophic at early times. (Color version online.)

In the text

Fig. 7 Time evolution of the number of planetesimals (top) and planetesimal surface density (bottom) as a function of planetesimal radius at the embryo’s radial bin. The number of planetesimals evolves by accretion of the embryo, migration, and planetesimal fragmentation. The generation of fragments (produced by the collisional evolution) does not compensate for the large loss of planetesimals with a radius of 0.1 km. (Color version online.)

In the text

Fig. 8 Time evolution of the relative velocity profiles (top) and the ratio profiles (bottom) for targets and projectiles of radii 100 km. The profiles represent the case for km and a disk six times more massive than the MMSN. Despite the higher values of the relative velocities of planetesimals with rp = 100 km, the ratio is lower compared with the case for small planetesimals. This is because the specific impact energy () for planetesimals of rp = 100 km is approximately three orders of magnitude higher than that for planetesimals of rp = 0.1 km. (Color version online.)

In the text

Fig. 9 Planetesimal accretion rates as a function of time for an embryo located at 5 au in a disk ten times more massive than the MMSN, where initially all the solid mass is deposited in planetesimals of 100 km of radius. The black dashed line corresponds to the case without planetesimal fragmentation. The red solid line corresponds to the total planetesimal accretion rate with planetesimal fragmentation. The accretion rate of planetesimals with a radius of 100 km is plotted as a green solid line. We also plot the accretion rates for other planetesimal sizes. (Color version online.)

In the text

Fig. 10 Time evolution of the radial surface density profiles of planetesimals with a radius of 100 km for a disk ten times more massive than the MMSN. Solid lines represent the case with planetesimal fragmentation, dashed lines represent the case without it. For this last case, the planet achieves the critical mass at 4.07 Myr. When planetesimal fragmentation is considered, the planet does not achieve the critical mass. (Color version online.)

In the text

Fig. 11 Planetesimal accretion rates as a function of time for an embryo located at 5 au in a disk ten times more massive than the MMSN, where initially all the solid mass is deposited in planetesimals with a radius of 100 km. The dashed black curve represents the total planetesimal accretion rate for the case without planetesimal fragmentation. The red solid curve represents the total planetesimal accretion rate when p = 5/3 is a fixed value for all collisions. We also plot the planetesimal accretion rates for the different planetesimal sizes for p = 5/3. (Color version online.)

In the text

Fig. 12 Time evolution of the number of planetesimals adopting km and a disk ten times more massive than the MMSN. Filled circles represent the case where the mass loss in collisions is deposited in smaller fragments, while open circles represent the case where the mass loss in collisions is deposited in larger fragments. For this last case, the planet achieved the critical mass at ~3.61 Myr. (Color version online.)

In the text

Fig. 13 Core masses (solid lines) and envelope masses (dashed lines) as a function of time for the cases with planetesimal fragmentation considering p = 5/3 (black lines) and without planetesimal fragmentation (red lines), for a disk ten times more massive than the MMSN (top) and a disk eight times more massive than the MMSN (bottom). For this last case, the planets did not achieve the critical mass after 6 Myr of evolution. (Color version online.)

In the text

Fig. 14 Top: planetesimal accretion rates for different planetesimal sizes, for a disk ten times more massive then the MMSN, in which p = 5/3 and km, as a function of time. The dashed black line represents the total planetesimal accretion rate for the case without planetesimal fragmentation (WF). The dashed red line represents the total planetesimal accretion rate for the case with planetesimal fragmentation in which p is calculated with Eq. (31), our basis model (BM). The solid red line represents the total planetesimal accretion rate for the case with planetesimal fragmentation with p = 5/3. Bottom: time evolution of the total planetesimal accretion rates (solid lines) and the accretion rates of fragments of 1 m (dashed lines) for the case where only planetesimal fragmentation is considered (red lines) and for the case in which planetesimal coagulation and fragmentation are considered (black lines). The plot corresponds to a disk ten times more massive then the MMSN and for the case of km. (Color version online.)

In the text

	Fig. 15 Factor as a function of the core mass for different planetesimal sizes. The plot corresponds to the case of a disk ten times more masive than the MMSN and km when planetesimal coagulation and fragmentation are considered. The black dashed line corresponds to the evolution of the term (RP/RH)1 / 2, where RP is the radius of the planet. (Color version online.)
In the text

Fig. 16 Top: time evolution of the radial profiles of the surface density of planetesimals of rp = 100 km when planetesimal fragmentation is considered with p = 5 / 3 and for a disk ten times more massive than the MMSN. Black dashed lines correspond to the isolated formation of a planet located at 5 au. Red solid lines correspond to the simultaneous formation of two planets located at 5 au and 6 au. The different lines correspond to different times: 0, 0.5, 1, 1.5, 2, and ~2.5 Myr. Eventually, the planetesimal surface density decreases in the neighborhood of both embryos. Bottom: same as top panel, but for fragments that most contribute into the total accretion rate, rp ~ 25 km. In this case, the profiles increase around 5 au for t = 0.5,1,1.5 Myr, but for t = 2, ~2.5 Myr profiles decrease because of accretion. The increments in the surface density around 6 au correspond to the profiles at times 1.5, 2, and ~2.5 Myr. (Color version online.)

In the text

Fig. 17 Core masses as a function of time for the isolated formation of a planet located at 5 au (black dashed line) and for the simultaneous formation of two planets located at 5 au (red solid line) and 6 au (blue solid line). Models correspond to km for a disk ten times more massive than the MMSN, and planetesimal fragmentation is considered with p = 5/3. The simultaneous formation stopped at ~2.5 Myr because the distance between planets became smaller than 3.5 mutual Hill radii. (Color version online.)

In the text