Free Access
Issue
A&A
Volume 632, December 2019
Article Number A14
Number of page(s) 16
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201936061
Published online 22 November 2019

© ESO 2019

1 Introduction

Understanding terrestrial planet formation is an ongoing challenge in planetary sciences. We know that planets are formed around stars as a natural by-product of the star formation process. We can broadly summarize it in different stages (Morbidelli et al. 2012). Once the star collapses, a disk of gas and dust forms. The dust settles into orbit about the host star via an accretion process in the disk. Then, as small grains or ices orbit, they collide at low relative velocity and stick, thereby forming aggregates (Weidenschilling 1977; Krijt et al. 2015).

While the entire mechanism is not yet fully understood, these aggregates eventually become bodies of a few kilometers across, called planetesimals. Once they reach this size, the planetesimals are decoupled from the gas and move on Keplerian orbits around the star. Mutual gravitational interactions and collisions then become important as their cross sections for accretion becomes larger than their physical size because their gravitational pull is sufficient to capture adjacent material. As long as the impacting velocities are lower than the escape velocity of the larger planetesimals, they can grow rapidly. When this happens, the planetesimals grow much faster than smaller bodies and lead to a new phase of accretion called runaway accretion (Kokubo & Ida 1996), resulting in planetary embryos or protoplanets of more than 50 km in size. Runaway growth breaks the uniformity of the velocity and spatial distributions of planetesimals (Kokubo & Ida 1998). The growth of larger embryos slows down when protoplanets grow large enough to increase the velocity dispersion of the nearby planetesimals, thus reducing the rate of accretion (Ida & Makino 1993). Under this scenario, several protoplanets are formed and dominate the planetesimals dynamical evolution. This stage is called oligarchic growth (Kokubo & Ida 1998). Finally, the last stage of planetary formation consists of an increase of close encounters that lead to giant impacts between protoplanets. The final configuration of the planetary system is defined at the end of this stage between 10 Myr and 150 Myr (Chambers & Wetherill 1998; Chambers 2001; Jacobson et al. 2014).

The N-body integrators are known in the study of the different stages of planetary formation. In particular, a great variety of this kind of integrators exists in the scientific community for different physics problems. The main goal of performing N-body simulations is to study the late stages of terrestrial planet formation. Several authors used this kind of integration to analyze the late stage of accretion for the Solar System (e.g., Chambers 2001; O’Brien et al. 2006; Raymond et al. 2009, to name a few) and for exoplanetary systems with different stellar types, and with or without gas giants (e.g., de Elía et al. 2013; Dugaro et al. 2016; Darriba et al. 2017; Zain et al. 2018; Sánchez et al. 2018, among others). In these studies, although the dynamical evolution was well represented, the collision treatment was simplified. Colliding planetary embryos perfectly merge into a new body conserving mass and water content. One of the reasons for this simple model is that the computation time in an N-body simulation scales with N2, where N is the number of particles in the system. We have to set a compromise between the number of bodies and performance. Without allowing fragmentation between planetary embryos, the number of bodies in the system is not increased. As collisions take place, the number of bodies drop and the computing time decreases.

Leinhardt & Stewart (2012, hereafter LS12) presented a complete analytic collision model for gravity-domained bodies. The range of outcomes that LS12 proposed includes partial accretion, partial erosion, catastrophic collisions, and hit-and-run encounters. These outcomes for each collision depend on the target size, projectile size, impact velocity, and impact angle. The LS12 authors define the catastrophic disruption criteria as the specific energy needed to disperse half of the involved mass in the collision. For this reason, the authors introduce two material parameters: as a coupling parameter and the dissipation parameter c*. With the analytic derivation and scaling laws presented in their work, LS12 described the transition between collision regimes and size-velocity distribution of the post-collision bodies. The analytic model developed by LS12 is a powerful tool that has two major advantages: it improves the physics of collisions in numerical simulations of planet formation and collisional evolution and it is easy to adapt to an N-body code.

On the other hand, Genda et al. (2012) studied the merging criteria for collisions of rocky embryos on terrestrial planet formation. The authors used a smoothed particle hydrodynamic method for the giant impacts and investigated the critical impact velocity, where transition between hit-and-run and merging takes place. The authors derived a simple formula for the normalized critical impact velocity vcr/, where is the two-body escape velocity, which depends on the mass ratio of the protoplanets and the impact angle. That said, a hit-and-run collision with an impact velocity lower than the critical value may result in a second collision leading to a perfect merge.

This improvement was adopted by several authors to study the Solar System and extrasolar systems. Chambers (2013) implemented a refined collisional algorithm into MERCURY based on the work developed by LS12. In that work he studied the final stage of terrestrial planet formation with N-body simulations, concluding that hit-and-run collisions are a common outcome of two colliding big bodies. The author came to this conclusion comparing the dynamical evolution of two identical sets of initial conditions, both considering and not considering fragmentation. He found two important differences between these two models. On the one hand, the planetary masses are lower when fragmentation is included. On the other hand, the excitation of the eccentricity over time is also lower in comparison with the model that considers only perfect merges.

Quintana et al. (2016) performed N-body simulations to study the terrestrial planet formation around a Sun-like star, using the modified version of MERCURY implemented on Chambers (2013). The authors compared the standard accretion model with the modified version of MERCURY (Chambers 2013), which allows fragmentation in the collisions between big bodies. The authors found that the collisional history of the planets that survive at the end of the integration differs significantly between the two models, although the overall masses and number of planets generated are comparable.

More recently, Mustill et al. (2018) carried out his own implementatation of the collisional algorithm developed by LS12 into the public version of MERCURY. In particular, these authors studied how this improvement in the treatment of the collisions affects the outcomes concerning the in situ planet formation on packed systems and instability scenarios. It is worth noting that, unlike the works carried out by Chambers (2013), Quintana et al. (2016), and Wallace et al. (2017), Mustill et al. (2018) introduced a factor of mass removal in the collisional algorithm to represent the material that is ground and then removed from the system by radiation forces.

With this in mind, it is evident that there is a need to have a numerical tool of our own that allows us to study, in a more detailed way, the collisional history of the bodies that make up a planetary system as well as its dynamical evolution. Thanks to the improvements in the collisional model proposed by LS12, we decided to move away from the perfect merging model and developed the D3 N-body code with a more realistic treatment of collisions between planetary embryos. This improvement allows us to carry out a detailed study about the final composition of the planets formed. In particular, we can study the water delivery in a more realistic way than in the classic models of accretion.

In this sense, Marcus et al. (2010) described two empirical models for the mantle stripping in differentiated planetary embryos after a collision. The authors used a simple differentiated structure for the terrestrial embryos, assuming an iron core and an ice-silicate mantle, and performed a series of SPH (Smoothed particle Hydrodynamics) simulations with the GADGET (GAlaxies with Dark matter and Gas intEracT) code (Springel 2005). In particular, Marcus et al. (2010) varied the projectile-to-target ratio, impact angle, and number of particles of each body to investigate how the giant impacts affect the final abundance of water in such bodies. From this, they concluded that the more violent the collisions, the more mantle is lost, and that water is more easily removed than the silicates.

More recently, Dvorak et al. (2015) obtained interesting results concerning to water content retained in significant fragments after a collision. From SPH simulations, these authors studied the outcome of an impact between two bodies with different values of impact velocity and impact angle. One relevant conclusion of their research suggests a significant water loss for faster and/or less inclined collisions.

The investigations developed by Marcus et al. (2010) and Dvorak et al. (2015) suggested that incorporating a realistic model of transport and removal of volatiles in an N-body code may lead to reduced water contents of the resulting terrestrial-like planets in comparison with those derived from classical models that assume perfect mergers.

The main goal of the present research is to study the physical and dynamical properties of terrestrial-like planets and water delivery in the habitable zone (HZ) using N-body simulations that incorporate fragmentation and hit-and-run collisions. According to this, this paper is therefore structured as follows. In Sect. 2, we present the main properties of the D3 N-body code integrator and the collisional model introduced by LS12. In Sect. 3, we describe the application and initial conditions used to carry out our simulations. In Sect. 4, we show our results and carry out a detailed analysis of all simulations and the planets remaining in the HZ. Finally, we discuss such results within the framework of current knowledge of planetary systems and the limitations of the model in Sect. 5.

2 N-body code: classic integrators

In this section, we describe some basic ideas of the symplectic theory and the general properties of the numerical integrator. We also present the different regimes implemented in the numerical integrator.

2.1 Overview

The solution to the N-body problem is obtained through the resolution of the equations of motion that derive from the Hamiltonian of the system. In effect, the Hamiltonian of a system of N + 1 particles of masses mi (i = 0, ..., N), subject solely to the action of their mutual gravitational attraction, is given by (1)

where qi = (xi, yi, zi) and are the position and momentum vectors of the ith body (of mass mi) with respect to an inertial reference frame, and G the universal gravitational constant.

However, for the particular problem of a planetary system, where most of the mass is concentrated in one body (e.g. the Sun in the Solar System), it is more convenient to divide the Hamiltonian as follows: (2)

where HKep is the part of the Hamiltonian that describes the Keplerian motion of the particles around the massive body and Hint is the part that describes the interaction between each pair of particles, except with the main body (i.e. the star).

Depending on the coordinates chosen, the Hamiltonian can be divided in different ways. Wisdom & Holman (1991), for instance, used Jacobian coordinates. This system constitutes a reference frame in which the position and momentum of each body are considered with respect to the center of mass of all the bodies with indices lower than a certain previously defined value. Symplectic integrators based on this division of the Hamiltonian are known as symplectic integrators of mixed variables, i.e., mixed variable symplectic methods (Saha & Tremaine 1992).

An alternative procedure proposed by Duncan et al. (1998) is to define the position with respect to the main body and the momentum with respect to the barycenter of the system. This set of variables is defined as heliocentric democratic variables. By means of a canonical transformation, it is possible to transform the original set of coordinates and momenta (q, p) into a new set of heliocentric coordinates (Q) and barycentric momenta (P).

Implementing this transformation, and adopting this new set of variables, the Hamiltonian from Eq. (1) could be rewritten in the following way (3)

Since in Eq. (3) the coordinates Qi are heliocentric, which corresponding to the central body vanishes, i.e. Q0 = 0. Consequently, the Hamiltonian is independent of the coordinate Q0, therefore the momentum P0 is a constant of motion, and its contribution to the Hamiltonian is not considered. It is worth noting that N indicates the number of particles excluding the star and m0 indicates its mass. Therefore, we can group the terms from Eq. (3) as (4)

where

each one of these terms describing the following:

  • HKep: the Keplerian motion of the particles around the massive body,

  • Hint: the interaction between each pair of particles (except with the main body),

  • H: the main body barycentric momentum, which appears due to the chosen set of variables.

Following Chambers (1999), by solving Hamilton’s equations of motion, the general solution to the rate of change of a given function, u = u(q, p), over a time τ, starting from its initial value u(0), is given by (8)

where H is the Hamiltonian operator. As shown in the aforementioned paper, if an operator H is written as the sum of two operators (i.e., H = HA + HB), then that in Eq. (8) can be split as the consecutive product of the two operators, (9)

By splitting Hint and H in half, we obtain a second-order integrator, and Eq. (8) becomes (10)

From the operator defined in Eq. (10), we can construct a five-step integration scheme as follows:

(i) the coordinates remain fixed and each body receives an acceleration from the other bodies (but not from the main body) that modifies its momentum over a time interval τ∕2; (ii) the momenta remain fixed, and each body undergoes a displacement in its position in the amount ; (iii) each body evolves around a Keplerian orbit (with the same central location and the same central mass) over the whole time interval τ; (iv) repeat step II; and (v) repeat step (i).

However, if two bodies are too close to each other, the corresponding term in Hint becomes large enough so HKep ceases to be the dominant term and the error increases substantially. For this reason, the previous scheme cannot resolve situations of close encounters.

One of the solutions proposed to address this problem (Duncan et al. 1998) consists in dividing the perturbative terms of Hint and assigning each of these a different integration step size, so that the strongest perturbations have the smallest step sizes. The resulting integrator is purely symplectic, although complicated to carry out, and also does not retain the high speed of the basic symplectic method. An alternative solution is proposed by Chambers (1999), who constructed a hybrid algorithm that mixes both symplectic and non-symplectic components to retain properties of both. As we mentioned before, we need to make the term Hint small again (so that HKep remains dominant). One way to do this is to transfer the close encounter term in Hint to HKep for the duration of the close encounter. For example, if the bodies α and β have a close encounter, the terms from Eqs. (5) and (6) result in

With this scheme, HKep can no longer be solved analytically. However, we can solve it numerically using a conventional integrator (briefly described later), and the remaining bodies that do not have close encounters can be analytically integrated. We needed to keep in mind that doing this brings amodification to the original Hamiltonian and the integrator no longer remains purely symplectic. In order to make the hybrid integrator symplectic, we had to make sure that no term is transferred between the different parts of the Hamiltonian. Following Chambers (1999), we were able to do this by dividing each term of interaction between HKep and Hint so that the corresponding part in Hint is always kept small, while the part in HKep is evaluated only during close encounter as follows

The form of the function K must be such that, when the separation |QiQj| between the bodies i and j is large, K must tend to one, while it tends to zero when such a separation is small. On the other hand, the transition has to be smooth enough. Following Chambers (1999), we adopted the function,

where (15)

and rcrit it is a free parameter that indicates the critical switching distance. Generally is a multiple of the mutual Hill radius of the involved bodies.

This ensures that that |Hint|≤|HKep|, even during a close encounter (for a detailed explanation, refer to Chambers 1999). Nevertheless, the choice of K is an open debate on how the integrator keeps its symplecticity (Hernandez 2019; Rein et al. 2019) about the switchover function and its expression is discussed in Sect. 5. Thus, the second-order hybrid integrator proceeds as follows:

(i) The coordinates remain fixed and each body receives an acceleration from the other ones (but not from the main body), weighted in the situations of close encounter by a factor K, which modifies its momentum over a time interval τ∕2; (ii) the momenta remain fixed, and each body undergoes a displacement in its position in the amount over all the N bodies (except the central star); (iii) the bodies that are not in close encounter move in a Keplerian orbit around the main body over a time interval τ. For those that are in close encounter, the Keplerian terms and the close encounter terms are weighted by (1 – K) and integrated numerically over a time interval τ; (iv) repeat step (ii); and (v) repeat step (i).

It is worth noting that the accretion problems that we focus in this work take place in step (iii), where the collision conditions are given as a result of the numerical integration.

Following Chambers (1999), we adopted the same numerical integrator for solving the close encounters. Because of its robustness, speed and precision, the Bulirsch–Stöer method is an excellent choice for this kind of numerical problems.

2.2 Collisions in classic integrators

Perfect merging is the simplest collision model used in simulations of planetary accretion studies so far, either to study the Solar System (Chambers 2001; O’Brien et al. 2006; Raymond et al. 2009) as well as extrasolar systems (de Elía et al. 2013; Dugaro et al. 2016; Darriba et al. 2017; Sánchez et al. 2018; Zain et al. 2018). This collisional model does not have physical and geometric parameters involved in the outcome of the collision: two bodies that collide result in a body with a mass equal to the sum of the masses of the target and the projectile.

Recently, LS12 performed high-resolution simulations of collisions between planetesimals. There, they found a wide range of possible outcomes: merging, cratering, super-catastrophic disruption, and hit-and-run events. The authors derived useful scaling laws that determined the transition between the different regimes. The different collision outcomes were determined by physical properties such as impact velocity, mass ratio of the bodies, and a geometric parameter like the impact angle.

It is important to remark that in the D3 N-body code embryos are treated as big bodies, meaning that they interact with any other body in the system. Planetesimals and fragments, on the contrary, are treated as small bodies, that is, they do not interact with each other, only feeling the gravitational pull from the embryos and the central star.

2.3 Collisional model: different type of collisions

In this section we briefly describe the regimes derived by LS12 and were used in this work. After this description, we present the algorithm to identify and calculate the different giant impacts.

Collision outcomes can be described in terms of the impact energy per unit mass Q and a critical value , where is defined as the specific energy per unit mass required to disperse half of the total colliding mass. Depending on the collision, it may result in a large body and several fragments. In this case, we reference as Mlr to the mass of the largest remnant of the mentioned collision, and the remaining mass is distributed in equal-mass fragments. According to Chambers (2013), the snumber of fragments generated is directly related to a parameter that indicates the minimum permitted fragment mass, Mmin. On the one hand, this value has to be high enough so the number of fragments generated are not too high and, thus, slow down the simulations. On the other hand, if the value is too high, the fragments generated are unrealistic. Therefore, the value Mmin must be set considering a compromise between number of fragments generated, performance and realism of the model.

We considered a collision between two planetary bodies as shown in Fig. 1. A target of mass Mt with a radius Rt and a projectile of mass mp with a radius rp, where Mt > mp, collide with an impact velocity vi and an impact angle θ, where θ is the angle formed between the line connecting the centers of the two bodies and the projectile velocity vector relative to the target. The bulk density for both the target and projectile are ρt and ρp, respectively. We describe the different types of collisions as follows:

  • Perfect merging: This type of collision is the most commonly used by the classic integrators so far. No fragments are generated and the remaining body contains the total colliding mass.

  • Partial accretion: When Mlr is calculated and is higher than the mass of the target but lower than the total colliding mass (i.e., MtMlrMt + Mp), the outcome is considered a partial accretion. The remaining mass, if it higher than the minimum fragment mass Mmin, is distributed in fragments. Otherwise, the collision is considered as a perfect merging.

  • Erosive collision: In this case, Mlr is lower than the mass of the target. The projectile is completely destroyed and, in addition to the remaining mass of the target, is distributed in fragments.

  • Super-catastrophic collision: In giant impacts, where the impact energy is so high that the Mlr is less than 10% of the total colliding mass, the outcome is considered a super-catastrophic collision, and the remaining mass is distributed in fragments.

  • Hit-and-run: This outcome occurs when both the target and the projectile remain intact. No fragments are generated.

  • Graze and merge: Depending on the velocity of a hit-and-run collision, this may result in a second collision leading to a perfect merge; there is no generation of fragments.

  • Erosive hit-and-run: Named as such because during a hit-and-run collision, the projectile may be partially or totally disrupted and the target comes out unaffected. The mass stripped out of the projectile is distributed in fragments.

With the different types of collisions described, in the next section we address the modification that has to be made to an N-body numerical integrator, to implement the collision classification.

thumbnail Fig. 1

Schematic representation of a collision between two planetary bodies. The target, with a mass Mt and a radius Rt collides with a projectile of mass mp and radius rp at an impactvelocity vi. The angle θ represents the angle formed between vi and the line connecting the two centers.

Open with DEXTER

2.4 Modification in the numerical integrator

We must have detailed information concerning each collision to calculate the results obtained by LS12 and apply the classification of the different regimes. More precisely, we need the position of the bodies, impact velocity, and impact angle at the time of the collision within a certain tolerance.

Following Chambers (1999), the collisions are detected with a pre-checker, based on a Hermite interpolation scheme. When this happens, following Mustill et al. (2018), we redid the integration in an iterative fashion, by successively dividing the integration step by half, until the bodies are in collision or the integration step is lower than 10−5 days. In thisway, we acquired information accurate enough to be able to resolve the collision and classify it with the regime scheme of LS12.

2.5 Collisional algorithm

The collisional algorithm implemented in this work is based on LS12 (for a detailed explanation, see the appendix in the aforementioned paper). This algorithm flow can be summarized in a series of steps, as detailed below

Step 1: when a collision is detected, because one of the bodies, namely the target, is necessarily an embryo (since interaction between planetesimals are neglected), the code first enquires about the body type of the second body involved, i.e., the projectile. If it is a planetesimal (or a fragment), the code assumes a perfect merge, preserving the total mass and momentum. Then, it proceeds to analyze the next collision, starting again on step 1. If the projectile is an embryo, it proceeds to the next step.

Step 2: in effect, the collision occurs between two embryos, the mutual escape velocity is calculated by means of the following expression , where We changed the definition of the mutual escape velocity proposed by LS12 and instead we used that proposed by Mustill et al. (2018). A comparison between the impact velocity vi and the escape velocity vesc is made. If vivesc, then the collision is also assumed to be a perfect merge, conserving the total mass and momentum, as in step 1. In the case vi > vesc, we proceed to the next step.

Step 3: a comparison between the impact velocity vi and the escape velocity vesc is made. If vivesc, then the collision is also assumed to be a perfect merge, conserving the total mass and momentum, as in step 1. In the case vi > vesc, we proceed to the next step.

Step 4: we calculated the impact parameter b = sinθ and the critical value bcrit = Rt∕(Rt + rp). Depending on the relation between b and bcrit the collision may result in either a grazing or a non-grazing impact. If b > bcrit the impact was grazing and we skipped to step 8. Otherwise, if b < bcrit we proceed to calculate the impact energy per unit mass, Q, and the catastrophic disruption criterion, . The value of Q is given by (16)

where μ = Mt mp∕(Mt + mp). Then, is calculated as follows (17)

where (18)

where α is the fraction of the projectile that intersects the target in an oblique impact, a measure of how energy and momentum from the projectile are coupled to the target, c* an energy dissipation factor within the target, ρ1 = 1 g cm−3, γ = mpMt, and Rc1 is the spherical radius of a body with mass Mtot = Mt + mp with density ρ1 (Stewart & Leinhardt 2009). According to such a work, in this case, and c* adopt values of one-third and 1.8, respectively.

Step 6: with the values of Q and we check, by simple arithmetic comparison, if the collision is in the super-catastrophic regime (Q > 1.8 QRD). This result affects on how to compute the mass of the largest remnant. Following LS12, if the collision is in such a regime, we calculate Mlr as (19)

and the mass to distribute in fragments is mfrag = MtotMlr. Step 7: In the case , the collision is not in the super-catastrophic regime, and the mass of the largest remnant is computed as (20)

By comparing the largest remnant’s mass Mlr and the target’s mass Mt, the collision may result in:

  • 1.

    if MlrMt, it’s a partial erosion,

  • 2.

    if Mlr > Mt, it can result in either a partial accretion, if the number of fragments is larger than 0, or a perfect accretion, if mfrag < Mmin, for which no fragments are generated.

Step 8: Grazing impacts: for grazing impacts, we checked whether the collision is a possible hit-and-run. If it is not classified as a hit-and-run, in order to classify the collision as super-catastrophic, partial erosion, or partial accretion of perfect accretion, we followed the same criteria used for non-grazing impacts, described on steps 5 through 7.

Step 9: in case the collision is classified as a possible hit-and-run, Genda et al. (2012) investigated the critical impact velocity vcr, that establishes the boundary between the pure hit-and-run and merging impacts. This velocity is given by (21)

where Γ = (1 − γ)∕(1 + γ) and Θ = 1 − sinθ. The fitting parameters are c1 = 2.43, c2 = −0.0408, c3 = 1.86, c4 = 1.08, and c5 = 2.5. These authors found that the normalized does not depend on Mt but rather on γ and the impact angle θ, where is the mutual escape velocity.

We applied the formula derived by Genda et al. (2012) and classified the hit-and-run at low velocity impacts as merging impacts. Chambers (2013) coined the term for this type of collision as graze-and-merge, meaning that a grazing impact at low velocity derives in a secondary collision where the two bodies merge into one, like a two-stage perfect merging.

Then, if the impact velocity satisfies that vi < vcr, we classifiedthe collision as a graze-and-merge.

Step 10: if the collision effectively derives in a hit-and-run (i.e., vi > vcr), the mass of the target is modified in a negligible amount, so considered Mlr = Mt (Asphaug et al. 2006; Genda et al. 2012) and we had to calculate the critical disruption energy for the reverse impact on the projectile (see, Sect. 4.2 of LS12). Thus, basically, the roles of the target and projectile are inverted. According to LS12, we set the reverse variables as follows: (22)

where η is the fraction of the target that interacts with the projectile (for a detailed explanation on how to get this fraction see Sect. 4.2 of LS12), and (23) (24) (25)

where Rc1 is given by (26)

The impact energy is (27)

Following Chambers (2013), the mass of the largest remnant for the projectile is calculated as follows (28)

Then, we have to calculate , and distribute it in equal-mass fragments. In the case mfrag < Mmin, we considered the collision as a pure hit-and-run, therefore the masses of both the target and projectile remain unaffected.

This algorithm described in the previous steps is used every time that a collision is detected. However, the computing time required to rerun the collision and classify it is negligible compared to the overall integration time.

3 Applications: setup

In this section, we describe the scenarios of work that are used to test our code.

For the present work, we focus on the study of a dynamical scenario aimed at exploring the sensitivity of the results to the presence of massive perturbers. The methodology to define the physical properties of the planetary embryos to be used in our model is explained as follows.

3.1 Protoplanetary disk: model

We describe the surface density profile that represents a relevant parameter of the properties of the planetary disk. Following Lynden-Bell & Pringle (1974) and Hartmann (1998), we adopted our model of protoplanetary disk based on the evolution of a thin Keplerian disk only ruled by the gravity of a point-mass central star. The gas-surface density profile Σg (R) is given by (29)

where R is the radial coordinate in the disk mid-plane, γ the exponent that represents the surface density gradient, Rc a characteristic radius, and a normalization constant.

If we integrate Eq. (29) over the total disk area and assuming axial symmetry, can be written as a function of Rc, γ, and the mass of the disk Md by (30)

where (31)

In the same way, we define a solid-surface density profile Σs (R) given by (32)

where is a normalization constant, and ηice is a parameter that represents an increase in the amount of solid material due to the condensation of water beyond the snow line.

In the present study, we assumed a central star with a mass M = 1 M and a solar metallicity ([Fe/H] = 0). From this, we considered that the relation between the gas and solid surface densities is given by , where z0 is the primordial abundance of heavy elements in the Sun and has a value of z0 = 0.0153 (Lodders et al. 2009). Moreover, we adopted a characteristic radius Rc of 25 au and an exponent γ = 0.9, which are in agreement with the median values derived from observations of different disks studied by Andrews et al. (2010) in the 1 Myr Ophiuchus star-forming region. Moreover, we adopted a disk mass of Md = 0.01 M, which is consistent with observational studies developed by several authors, such as Andrews et al. (2010) and Testi et al. (2016). Finally, we remark that the snow line is assumed to be located at 2.7 au in the present model, according toIda & Lin (2004). Following Lodders (2003) and Lodders et al. (2009), we assumed that the parameter ηice adopts values of 1 and 2 inside and beyond the snow line, respectively.

The increase in the amount of solid material from the condensation of water beyond the snow line produces a radial compositional gradient in the protoplanetary disk. Thus, we propose that the initial fraction of water by mass of the material that composes such a disk is a function of the radial coordinate in the mid-plane R and is given by

It is worth noting that this distribution considers that the inner region of the system was populated with water-rich material from the outer regions during the gaseous phase associated with the evolution of the disk, which is consistent with that proposed by Raymond & Izidoro (2017).

Finally, we define the HZ of the system as the region around a star in which a planet could retain liquid water on its surface. In this sense, Kopparapu et al. (2013) established inner and outer limits for the HZ around stars of different spectral types. In particular, the investigation developed later in the present work uses the optimistic estimates for a solar-type star derived by those authors, where the inner (outer) edge is located at 0.75 (1.7) au.

In the following section, we describe the necessary parameters for performing the N-body simulations of our research.

3.2 N-body simulations: physical and orbital parameters

For the development of our simulations, we have to specify the physical and orbital parameters of the bodies involved. For that reason, we start by specifying the mass distribution for the planetary embryos.

From Eq. (32), we can determine the mass distribution once the gas in the disk is fully dissipated. The extension of the region of study is between 0.5 au < R < 5.0 au.

Following Kokubo & Ida (2000), we can determine the mass of an embryo growing in the oligarchic growth regime located at a distance R from the central star as (33)

where the factor p represents the ratio of the mass of the embryos to the total mass of the system, δRH is the orbital separation of two consecutive planetary embryos of mass M in terms of their mutual Hill radii given by (34)

and δ is a parameter randomly chosen following a uniform distribution. The value for δ is generated each time the mass and spacing of an embryo are calculated, ranging between 5 and 10.

Combining Eqs. (32) and (34) in Eq. (33), we can derive a function that relates the mass of each embryo with the distance R as follows (35)

Our scenarios of work only consider planetary embryos for which, we set the value of the factor p = 1.0, assuming that there are no planetesimals left at the end of the gaseous phase.

The first embryo is assigned to have an initial semimajor axis a1 = 0.5 au and a mass of M1 = 0.036 M (Eq. (35)). We calculated the initial semimajor axes and masses for the remaining planetary embryos as follows (36) (37)

We repeated this procedure until the outer limit of our region of study (5.0 au) is reached, obtaining ~ 50 embryos and a total mass distributed of 12.8 M. Figure 2 illustrates the distribution of mass of each embryo as a function of the distance (R) from the central star once the gas is fully dissipated. As for the initial orbital parameters, we set the values of eccentricity and inclination randomly with a maximum value of 0.02 and 0.5°, respectively. As for the argument of pericenter ω, longitude of ascending node Ω, and the mean anomaly M, the values were taken randomly between 0° and 360°.

To have a dispersive scenario that favors collisions between planetary embryos, we added two giant planets with masses and physical densities analogous to those of Jupiter and Saturn. Moreover, such planets are assumed to be located on their current orbits, which have semimajor axes of 5.2 and 9.5 au, and eccentricities of 0.0489 and 0.0565 for Jupiter and Saturn, respectively.

The essential part of this investigation is to study how the treatment of collisions changes the dynamical behavior and the evolution of a planetary system for a given set of initial parameters. We perform a total of 46 simulations using the D3 N-body code, 23 of which take into account a realistic collision prescription and 23 simulations that only consider perfectly inelastic collisions between planetary embryos.

As we mentioned in Sect. 2.3, the value of Mmin is a compromise between computational time and realism. From this, we adopted the value of 0.018 M. A discussion about Mmin and how thevariation of this value could affect the overall integration is carried out in Sect. 5. The time step used in this work is two days, which is shorter than 1/40th of the orbital period of the innermost planetary embryo in our simulations. We considered that a embryo collided with the central star if the distance is lower than 0.1 au. We adopted this nonrealistic stellar radius to avoid numerical errors in bodies with very small perihelion. Moreover, a body is considered ejected of the system if the distance is greater than 1000 au. Finally, the integration time for all simulations was of 200 Myr.

The main results derived from our investigation are presented in the next section.

thumbnail Fig. 2

Example of a mass distribution of planetary embryos from one simulation as a function of the initial semimajor axis at the end of the gaseous phase. The size of the points are scaled with the mass of each embryo: the different colors illustrate the initial fraction of water by mass for the embryos. In fact, orange (dark blue) indicates a fraction of 10−4 (0.5) water by mass. The sky blue shaded region represents the HZ for a solar-type star, while the dashed blue line illustrates the snow line assumed in our model.

Open with DEXTER

4 Simulation results

In this section, we present a detailed study about the general results obtained from numerical simulations that model the evolution of the dynamical scenario presented in Sect. 3. In particular, we carry out a comparative analysis between N-body simulations that incorporate fragmentation and hit-and-run collisions and those that assume that all impacts lead to perfect mergers. This investigation focuses on the physical and dynamical properties of the terrestrial-like planets produced in the HZ of the systems under study, analyzing the role of the fragments in their evolutionary histories.

thumbnail Fig. 3

Histogram of percentage of collisions for different regimes implemented in the D3 N-body code for the scenario under study.

Open with DEXTER

4.1 General analysis

The proposed scenario was constructed with the goal of studying the dynamical evolution of a planetary system under formation, subject to the perturbations of two giant planets. As we mentioned in Sect. 3, this was carried out by means of a more realistic treatment of the collisions between planetary embryos. This implementation might lead to a difference in the collision history of a planet during its formation. In order to better comprehend this, we analyzed the sensitivity in the collision type distribution.

In Fig. 3 we show the percentage of each collision regime for the mentioned scenario. For this figure, only the collisions between massive bodies were taken into account, ignoring those between planetary embryos and fragments, which are always considered as a perfect accretion.

We can observe that the majority of collisions are clumped in two big groups, i.e. perfect merge collisions (~ 25%) and hit-and-run encounters (~44%). This outcome is consistent with the work of Kokubo & Genda (2010), who found that 49% of the collisions between planetary embryos are hit-and-run impacts. Also, our result is consistent with Chambers (2013), who concluded that, when a refined treatment of collisions is included, about 42% of the total number of collisions embryo-embryo results in hit-and-run encounters. In addition to this, we obtained that a considerable percentage of collisions are not either perfect mergers nor hit-and-run impacts. More specifically, this percentage is about ~ 30%, including partial accretion, erosive, and super-catastrophic collisions. This wide range of different collisions leads us to consider how these results affect the amount of bodies generated and the overall evolution of the system.

We can observe in Fig. 4 the total number of bodies as a function of the integration time for a simulation considering only perfectly inelastic collisions and another simulation that assumes a realistic collision treatment. The chosen simulations are representative of each set of runs (fragmentation enabled and fragmentation disabled). The cyan curve shows the decrease of the number of planetary embryos in that simulation represented by the blue curve, but without the fragments generated ineach collision.

A natural result is that simulations that include fragmentation take more time for the bodies to decrease in comparison with simulations with perfect accretion. When fragmentation is activated, the time to get half of the initial number of planetary embryos is approximately ~ 11 Myr, while that time decreases down to ~5 Myr in runs without fragmentation. This result is a direct consequence of the large number of hit-and-run collisions, which leads to the embryos taking more time to collide in comparison with the classical model where this type of encounter isnot possible. With this in mind, it is evident that the accretion timescales for systems including a realistic collision prescription are longer than those that only consider perfect mergers, such as presented in Chambers (2013).

It is worth noting that the total amount of fragments generated in the simulations range between 50–130. Moreover, we remark that 2–22 fragments were generated from super-catastrophic collisions, while the partial accretions and erosive impacts produced between 1–21 fragments. These numbers could be misleading at first. In fact, a basic assumption is that a super-catastrophic collision should generate more fragments than other impact types, since a minimum of 90% of the colliding mass is available to be distributed. However, the number of generated fragments is sensitive not only to the collision regime, but also to the masses of the bodies involved in it. For instance, in the same run, a partial accretion between two bodies with masses of 0.92 M and 0.88 M generates 21 fragments, while a super-catastrophic collision between two smaller embryos with masses of 0.16 M and 0.12 M only produces 11 fragments.

In light of these results, we are interested in studying how this realistic collision treatment affects the final architecture of the system and the mass evolution of the resulting planets throughout the entire simulation. Figure 5 illustrates the final architecture for each N-body simulation using the fragmentation improvement model (right panel) and the classical model of perfect accretion (left panel) after 200 Myr of evolution. The filled black circles represent the final planets and the point size is scaled with their masses, as shown on top of each panel for 0.1 M, 0.5 M, 1 M, and 2.5 M. The sky blue shaded area indicates the limits of the HZ. We can observe that, in general terms, the number of planets for both sets of runs remains similar, which is consistent with Chambers (2013), who showed that the final architecture seems insensitive to the collisional model. However, it is important to remark that, although the number of final planets does not differ significantly between the two models, an increase in the number of systems with more than two planets is more evident in simulations that incorporate fragmentation and hit-and-run impacts. In fact, 11 runs with fragmentation produce systems with three or more planets, while only five systems have three planets in runs without fragmentation. This result should be interpreted carefully, since a longer integration time associated with simulations with fragmentation might lead to systems with a lower number of final planets.

A distinctive feature observed in almost all simulations of the right panel of Fig. 5 indicates the significant existence of low-mass planets around the inner edge of our region of study in comparison with other zones of the system. This feature is not observed in simulations without fragmentation, which are illustrated in the left panel of Fig. 5. This result is also observed in Fig. 6, which represents the final masses of the planets formed in runs with fragmentation and without it as a function of the semimajor axis.

In general terms, Fig. 6 shows that the maximum values associated with the overall mass distribution of the final planets formed in simulations without fragmentation tend to be greater than those produced in runs with a realistic collision treatment, reaching 2.95 M. This result should not be surprising since planets formed in simulations with fragmentation could lose mass in each collision while in simulations that assumed that all impacts lead to perfect mergers, each collision that occurs makes the planetary embryo grow. In the same way, the minimum values observed in the overall mass distribution illustrated in Fig. 6 are determined by the blue circles, which are associated with the planets formed in runs with fragmentation and reach 0.08 M. It is important to remark that the difference between the minimum and maximum values of the planetary masses associated with runs with and without fragmentation is enhanced between 0.5 au and 1.0 au. In such a region, we can observe a large concentration of bodies with masses greater than1 M for planets formed in simulations without fragmentation, while the planetary masses derived in runs with fragmentation are distributed below those values. On the other hand, in the region comprised between 1 au and 2 au, the masses of the largest planets produced in both sets of simulations are comparable. These results are consistent with those observed by Chambers (2013), who found similar results concerning the final masses of the surviving planets concluding that, when fragmentation is enabled, more low-mass objects are formed, and that the effects of fragmentation over the masses of largest planets are less notorious.

Finally, Fig. 6 shows that there is a significant amount of planets produced in runs with and without fragmentation in the HZ of the system, which is located between 0.75 and 1.7 au. The study concerning the evolution, survival, and physicalproperties of such a planets is one of the most important scopes of this work. In the following section, we describe this topic inmore detail.

thumbnail Fig. 4

Number of remaining bodies vs. integration time. The red curve corresponds to a simulation that lead to perfect mergers. The blue and cyan curves represent a simulation with the realistic collision treatment. In particular, the cyan curve only shows the decrease of planetary embryos. Both curves (red and blue) are representative for each set of simulations.

Open with DEXTER

4.2 Habitable zone planets

The study of the physical and dynamical properties of the planets formed in the HZ of the system is very important owing to the astrobiological interest of that kind of planets, we can observe in Fig. 7 the mass as a function of the semimajor axis of all planets that survive in the HZ of the systems of study. In particular, the red circles illustrate those planets that result from the N-body simulations that assume that all collisions lead to perfect mergers, while the blue circles show the planets produced from those N-body simulations that incorporate collisional fragmentation and hit-and-run collisions. In general terms, the N-body experimentswith fragmentation and hit-and-run collisions form HZ planets with somewhat lower final masses. In fact, the simulations without (with) fragmentation produce HZ planets whose final masses range from 0.42 M (0.34 M) to 2.68 M (2.31 M).

It is worth to mentioning that the planets that are formed and that do not belong to the HZ present different physical properties. In particular, the range of masses of the planets in the inner region of the HZ (a < 0.75 au) is between 0.08 M and 1.1 M in runs with fragmentation, while for the simulations in which we only consider collisions as perfect mergers the range of mass between 0.39 M and 1.5 M. As for the outer region of the HZ (1.7 au < a < 3.0 au), the planets present a masses in the range of 0.1 M–2.8 M in runs with fragmentation and 0.18 M–2.95 M for simulations without it. We can observe that the upper limit of the masses increases as we move outside from the central star.

It is important to remark that all our simulations produce two different kinds of planets in the HZ, depending on the initial location of their accretion seeds. In simulations without fragmentation the accretion seed is defined as the largest body in each collision (Raymond et al. 2009). For simulations with fragmentation the accretion seed is defined as the largest body in each collision where only one embryo survives. In the case in which the two embryos survive (pure and erosive hit and run), both of these keep their original seed. According to this, we refer to those planets whose accretion seed starts the simulation inside (beyond) the snow line of the system as class A (class B) HZ planets. From the distribution assumed in Sect. 3.1, the class A HZ planets have very low primordial water contents by mass, while the class B HZ planets have highly significant primordial water contents. Our study produces 6 (11) and 17 (14) class A and class B HZ planets, respectively, in N-body simulations with (without) fragmentation. This information can be found summarized in Table 1.

The blue and red circles illustrated in the left panel of Fig. 7 show the final mass of the class A HZ planets, which result from N-body experiments with and without fragmentation, respectively. Such planets show a broad range of final masses, which reach a minimum and a maximum value of 0.34M (0.42 M) and 1.87 M (2.30 M) in N-body simulations with (without) fragmentation, respectively. According to this, the minimum and maximum values of the final mass of the class A HZ planets resulting from N-body runs with fragmentation are about 20% smaller than those produced without fragmentation.

In order to understand such differences in the final masses, we carried out an analysis of the collisions involved in the evolutionary history of the class A HZ planets formed in N-body experiments without fragmentation and in those that incorporate collisional fragmentation and hit-and-run collisions.

In runs without fragmentation, the 11 class A HZ planets underwent a total of 70 perfect mergers. While one of these did not experience any collision during 200 Myr, and the two other such planets underwent between 1 and 2 perfect mergers, the other 8 class A HZ planets experienced more than 7 perfectly inelastic collisions, which led to significant values in the relative growth of their masseswith respect to the initial values.

It is important to mention that the collisions that are taken into account are only those with the surviving embryos. We do not count the possible collisions that the projectiles may have had.

This result can be observed in the top and right panel of Fig. 8, which illustrates the fraction contributed by the initial mass (violet) and perfect mergers with embryos (green) to the final mass of the class A HZ planets. According to this, the 8 most massive class A HZ planets with M ≳ 1 M produced in N-body simulations without fragmentation reach more than 90% of their final masses from perfectly inelastic collisions with planetary embryos over 200 Myr of evolution.

In runs with fragmentation, the 6 class A HZ planets produced in such simulations experience a total of 57 collisions, of which 9 are perfect mergers with planetary embryos, 27 perfect mergers with fragments, 6 partial accretions, 2 erosive impacts, 9 hit-and-runs, and 4 erosive hit-and-runs. According to this, the number of partial accretions and perfect mergers with planetary embryos is relatively low in the evolutionary history of the class A HZ planets in comparison with the number of perfect mergers with generated fragments. In fact, partial accretions and perfect mergers with planetary embryos (fragments) represent about 10 and 16% (47%) of the total number of collisions experienced by those planets, respectively. However, it is important to remark that the perfect mergers with generated fragments play a secondary role in the growth of class A HZ planets. In fact, as the top left panels of Fig. 8 shows through sky blue boxes, the fragments contribute with less than 30% to the final mass of those planets. An important difference observed with respect to the simulations without fragmentation is associated with the relative growth of the class A HZ planets more massive than 1 M. In fact, as shown in the top and left panel of Fig. 8 in violet boxes, the initial masses of such planets represent more than 70% of their final masses, which indicates that their relative growths are significantly less than those associated with the class A HZ planets with M ≳ 1 M produced in runs without fragmentation. Our results seem to suggest that the most massive class A HZ planets in simulations with fragmentation require relatively high initial masses.

Similarly, the blue and red circles represented in the right panel of Fig. 7 illustrate the final mass of the class B HZ planets, which are formed from N-body experiments with and without fragmentation, respectively. In this case, the resulting planets show a more narrow range of final masses. In fact, the least and most massive class B HZ planets produced in runs with (without) fragmentation have 0.97 M (1.17 M) and 2.31 M (2.68 M), respectively. In this case, the minimum and maximum values of the final mass of the class B HZ planets produced from N-body runs with fragmentation are about 15% smaller than those obtained without fragmentation. Again, it is necessary tocarry out an analysis of the collisions undergone by the class B HZ planets produced in our N-body simulations. In runs without fragmentation, the 14 class B HZ planets experience a total of 27 perfect mergers with embryos. While 2 of these do not undergo any collision over 200 Myr of evolution, the other 12 class B HZ planets experience between 1 and 4 perfectly inelastic collisions. According to that observed in the bottom and right panel of Fig. 8 through violet boxes, the initial masses of the class B HZ planets produced in runs without fragmentation represent more than 50% of their final masses. From this, the relative growth of such planets is less significant than that experienced by the class A HZ planets with M ≳ 1 M formed in simulations without fragmentation.

In runs with fragmentation, the 17 class B HZ planets formed in such simulations undergo a total of 93 collisions, of which 33 are perfect mergers with planetary embryos, 34 perfect mergers with fragments, 9 partial accretions, 1 erosive impact, 6 hit-and-runs, and 10 erosive hit-and-runs. This information summarized in Table 2. In this case, the partial accretions and the perfect mergers with planetary embryos (fragments) represent about 10 and 35% (37%) of the total number of collisions undergone by the class B HZ planets, respectively. As commented on the class A HZ planets, the contribution of the generated fragments to the final mass of the class B HZ planets is also secondary in this case, since they provide less than 10%, which is shown in the bottom and left panel of Fig. 8 through sky blue boxes. It is worth mentioning that, 3 of the class B HZ planets formed in runs with fragmentation do not undergo any perfect merger with an embryo over 200 Myr of evolution: however the other 14 experience between 1 and 5 of such perfectly inelastic collisions, which is comparable to that described for the class B HZ planets produced in simulations without fragmentation. From this and thelow number of collisions that lead to mass losses, the fraction contributed by the initial mass and perfect mergers with embryos to the final mass of the class B HZ planets is similar in simulations with and without fragmentation. This comparative analysis can be observed from the results illustrated in the bottom panels of Fig. 8, where the contribution provided by the initial mass and the perfect mergers with embryos is represented in violet and green boxes, respectively. It is worth mentioning that for the class A HZ planets in simulations with and without fragmentation their feeding zones are fully contained in the region inner to the snow line. While for class B HZ planets, their feeding zones are mainly contained in the outer region of the snow line, although there’s a minor contribution from the inner region. This result, affects to the final water content of the final bodies.

The final water content of the class A and class B HZ planets plays a key role in understanding their potential astrobiological interest. To quantify the final fraction of water by mass of such planets, we tracked their collisional histories throughout the entire evolution. For N-body simulations without fragmentation, the treatment is very simple since all collisions are assumed to result in perfect mergers, which conserve the mass and water content of the interacting bodies. For N-body experiments that incorporate collisional fragmentation and hit-and-run collisions, the treatment is more complex since we must specify the amount of water acquired by the largest remnant and the fragments generated in each event. To do this, we define the following two models, which are based on prescriptions proposed by Marcus et al. (2010):

  • Model 1: The fraction of water of the largest remnant of each collision is determined by assuming that the mass that escapes is water from the projectile (first) and the target (second), and then rocky material from such bodies in the same order.

  • Model 2: The fraction of water of the largest remnant of each collision is determined by assuming that the mass that escapes is water (first) and rocky material (second) from the projectile, and then water and rocky material from the target in the same order.

We note that Model 1 leads to the largest remnant with a lower fraction of water by mass than that derived from Model 2. However, it is worth mentioning that Model 1 produces a distribution of fragments with a greater amount of water than that generated from Model 2. Thus, the final water content of the HZ planets strongly depends on the class of impacts that they undergo throughout the entire evolution.

The left panel of Fig. 9 illustrates the final fraction of water by mass as a function of the final mass of the class A HZ planets formed in N-body simulations with and without fragmentation, which are represented by blue and red circles, respectively. In particular, the open and filled blue circles show the results derived using Models 1 and 2 to determine the water transfer in runs with fragmentation, respectively.

As shown in the left panel of Fig. 9, all class A HZ planets produced in N-body runs without fragmentation conserve theinitial value of 10−4 associated with the fraction of water by mass over 200 Myr of evolution, since all collisions are assumed to result in perfect mergers and their feeding zones are only associated with the region of the system inside the snow line. Moreover, we note that the primordial water contents of 8 of the 11 class A HZ planets formed in runs without fragmentation represent less than 10% of their final water contents. This can be inferred from the top right panel of Fig. 8. In fact, since the feeding zones of such planets are restricted to the region of the system inside the snow line, the contributions to the mass fraction normalized to the final mass are equivalent to the contributions to the water fraction normalized to the final water content. Thus, in general terms, in our N-body experiments without fragmentation, the main source of water of the class A HZ planets are through collisions with embryos, instead of from its primordial water content.

The analysis is more complex for the class A HZ planets that result from simulations with fragmentation. In such a case, the final fractions of water of the class A HZ planets depend on the model adopted to track the evolution of water. In fact, it is important to remark that the perfect mergers with generated fragments represent about 47% of the total number of collisions experienced by the class A HZ planets. Thus, while the fragments play a secondary role in the final masses of such planets, they can be important in determining the final fraction of water of such planets.

As we have already mentioned above, the left panel of Fig. 9 particularly illustrates the final fraction of water by mass as a function of the final mass of the class A HZ planets formed in simulations with fragmentation. According to this, the final fractions of water by mass derived from Model 1 show significant differences with respect to those obtained using Model 2. In fact, Model 1 offers a wide range of results. Indeed, such a model can lead to slight increases of ~ 2.4 × 10−4 and significant decreases of ~6.5 × 10−6 in the final fraction of water by mass of the class A HZ planets with respect to their initial fraction of water of 10−4. Moreover, Model 1 can also keep the final fraction of water by mass of some of such planets close to the initial value of 10−4. On the contrary, Model 2 shows more conservative results concerning the final fraction of water. From such a model, four of six class A HZ planets show a final fraction of water by mass very close to the initial value of 10−4. The other two planets undergo a slight decrease in the fraction of water by mass over 200 Myr with respect to the initial fraction, reaching a minimum final value of ~5 × 10−5, which is associated with the less massive class A HZ planet.

According to this, in general terms, the water distribution between the largest remnant and the generated fragments in a given collision proposed by Model 2 offers results are comparable to those obtained in simulations without fragmentation concerning the final fraction of water by mass of the class A HZ planets. In this sense, Model 1 shows more diverse results regarding the final fraction of water by mass of such planets, leading to slight increases, significant decreases, and even similar fractions with respect to those derived in runs without fragmentation.

It is very interesting to determine the main source of water of the class A HZ planets formed in simulations with fragmentation for each of the two models proposed to track the evolution of water. In order to analyze this, the left and right panels of Fig. 10 illustrate the contribution of the primordial water content, and that provided by collisions with planetary embryos and generated fragments to the final fraction of water by mass of the class A HZ planets, which is derived making use of Models 1 and 2, respectively.

From this, it is possible to observe that, in general terms, the role of the generated fragments is very important in the final fraction of water by mass of the class A HZ planets when Model 1 is used. In fact, the fragments represent the main source of water in three of six of such planets, providing the totality of their final contents. In two of six class A HZ planets, about 40% of their final water contents is acquired by perfect mergers with fragments. These examples are very important since they extend the significance of the fragments in cases in which the planets receive a relevant contribution of water from collisions with planetary embryos. Finally, it is worth noting that the fragments only play a secondary role in the final water content of the most massive class A HZ planet, in which the primordial water content represents about 70% of its final content. Beyond this, in general terms, it is necessary to note that the primordial water content does not play a relevant role in the final water content of the class A HZ planets formed in runs with fragmentation when Model 1 is adopted.

The conclusion concerning the contributions of the primordial water, planetary embryos, and generated fragments to the final water content of the class A HZ planets produced in simulations with fragmentation is very different when Model 2 is used. In this particular case, the general result shows that the generated fragments play a minor role in the final fraction of water by mass of such planets in comparison with that described using Model 1. In fact, this is clearly observed if we focus the analysis on those three class A HZ planets that acquire the totality of their final water contents by perfect mergers with fragments using Model 1. If Model 2 is adopted, the fragments only provide about 10% (30%) of the final water content of two (one) of those three class A HZ planets. In general terms, it is important to note that the primordial water content has a relevant contribution to the final water content in four of six class A HZ planets, which represents another important difference with that derived from Model 1.

Our results show that the main source of water and the final fraction of water by mass of the class A HZ planets produced in simulations with fragmentation are very sensitive to the model adopted to determine the distribution of water between the largest remnant and the generated fragments in a collision. According to this, it is very necessary to specify a realistic model of volatile transport to determine in detail the water abundances of the formed terrestrial-like planets from N-body experiments.

The situation is very different when the final water content of the class B HZ planets is analyzed. In fact, as illustrated in the right panel of Fig. 9, the final fractions of water by mass of the class B HZ planets produced in runs with fragmentation are represented by open (filled) blue circles and range from 15 (25) to 50% (50%) when Model 1 (2) is used. From this, we derive two important conclusions. On the one hand, in general terms, the final fraction of water by mass of the class B HZ planets does not strongly depend on the model adopted to determine the distribution of water in a given collision. Indeed, 12 of 17 of such planets reach almost the same final fraction of water by mass with both models. On the other hand, our results show that the collisional fragmentation is not a barrier to the formation and survival of water worlds in the HZ of the system. Another interesting result represented in the right panel of Fig. 9 indicates that the final fractions of water by mass of the class B HZ planets produced in runs without fragmentation show similar values to those obtained in simulations with fragmentation. This point is very important since the results of simpler N-body simulations without fragmentation concerning the water content of the so-called water worlds should provide a good first approximation to the real data.

thumbnail Fig. 5

Final outcomes at 200 Myr of 46 N-body simulations carried out with the D3 code. Left panel: 23 planetary systems that resulted from runs without fragmentation, while right panel: 23 systems produced from runs that incorporate fragmentation and hit-and-run collisions. The circle size scales with the planet mass according to that indicated in the top region of the figure. The sky blue area observed in both panels illustrates the HZ associated with a solar-type star (Kopparapu et al. 2013).

Open with DEXTER
thumbnail Fig. 6

Final mass as a function of the semimajor axis of the planets formed in simulations when fragmentation is activated (blue circles) and without fragmentation (red circles) fragmentation over 200 Myr of evolution.

Open with DEXTER
thumbnail Fig. 7

Final mass as a function of the semimajor axis of the class A (left panel) and class B (right panel) HZ planets produced in simulations with (blue circles) and without (red circles) fragmentation over 200 Myr of evolution.

Open with DEXTER
Table 1

Summary of amount of bodies formed for class A and class B HZ planets and their range of masses and final water fraction, in runs with and without fragmentation.

thumbnail Fig. 8

Distribution of the final mass for the surviving planets in the HZ. The violet bars represent the mass fraction corresponding to the initial mass of the body. The green and sky blue bars represent the mass fraction contribution from partial (or total) accretion of embryos and fragments, respectively. Top panels: the results for the class A HZ planets for Frag (left) and No Frag (right) simulations. Bottom panels: top panels but for the class B HZ planets. The height of each bar is normalized to the final mass of each planet.

Open with DEXTER
Table 2

Total amount and type of collisions ocurred for class A and class B HZ planets.

thumbnail Fig. 9

Final fraction of water by mass as a function of the final mass of the class A HZ planets (left panel) and class B HZ planets (right panel). The red filled circles illustrate the values associated with those planets formed in simulations without fragmentation, while the open and filled blue circles represent the values obtained from Models 1 and 2 in runs with fragmentation, respectively. In both panels, the black solid line represents the initial water fraction for class A and class B HZ planets.

Open with DEXTER
thumbnail Fig. 10

Same as Fig. 8 for the water contribution. The color code is the same as in Fig. 8. The panels represent the water contribution to the class A HZ planets using Model 1 (left panel) and 2 (right panel), respectively.

Open with DEXTER

5 Discussion and conclusions

In the present research, we analyze the formation and evolution of terrestrial-like planets around solar-type stars, in systems that host two giant planets with physical properties and orbital parameters similar to those associated with Jupiter and Saturn. In particular, our investigation is based on a comparative analysis of 46 N-body simulations, which are carried out using an own numerical code called D3 N-body code. This numerical tool incorporates fragmentation and hit-and-run collisions according to the prescriptions described in the works developed by Leinhardt & Stewart (2012), Genda et al. (2012), Chambers (2013), and Mustill et al. (2018). To develop our study, 23 out of 46 N-body simulations are carried out assuming a realistic treatment for collisions, while the other 23 N-body experiments are developed considering that all collisions lead to perfect mergers. In particular, our study focuses on the physical and dynamical properties of the terrestrial-like planets and water delivery in the HZ of the system.

It is worth noting that, unlike previous work such as Chambers (2013), Quintana et al. (2016), and Wallace et al. (2017) in our N-body simulations we did not include planetesimals for studying the evolution of a planetary system. Our protoplanetary disks are only composed only by planetary embryos. We understand that, to have real initial conditions and to deermine the mass ratio of planetesimals to embryos, it is necessary to have a detailed evolution of the gaseous stage. We propose a simple model. In this model we can compare the physical properties of terrestrial-like planets in N-body simulations with and without fragmentation.

Our comparative analysis between simulations with and without fragmentation derives results consistent with those obtained by Chambers (2013), concerning the general percentage of the different collisional regimes, temporal evolution of the number of objects of the system, final mass of the terrestrial-like planets, and final planetary architecture.

The physical properties of the terrestrial-like planets formed in the HZ of systems that result from simulations with and withoutfragmentation show significant differences in several aspects. In particular, our conclusions vary for class A and class B HZ planets, which show accretion seeds initially located inside and beyond the snow line, respectively.

In general terms, planetary embryos represent the main source of mass and water of the class A HZ planets produced in runs without fragmentation. In particular, the class A HZ planets more massive than 1 M experience very significant increases of mass and water with respect to their initial values. Beyond this, the final fraction of water of all class A HZ planets formed in simulations without fragmentation is equal to their initial value of 10−4, since all collisions are treated as perfect mergers and their feeding zones do not include water-rich embryos, initially located beyond the snow line.

As for the simulations that incorporate fragmentation and hit-and-run collisions, the generated fragments play a secondary role in the mass of the class A HZ planets, providing less than about 30% of their final value. In particular, the class A HZ planets more massive than 1 M undergo small increases of mass with respect to their initial values, since the primordial mass contributes with more than 70% to the final mass.

It is important to note that the role of the fragments in the final water content of the class A HZ planets, as well as their final fractions of water, strongly depend on the model adopted to determine the distribution of water in a given particular collision. On the one hand, if water from the projectile and the target escape first, the generated fragments play a very important role in the final water content of the class A HZ planets, while the primordial water content does not represent a significant fraction of the final content. Moreover, the final fraction of water of such planets covers a wide range of values, which show slight increases, significant decreases, and even very low deviations with respect to the initial fraction of 10−4. On the other hand, if water and rocky material from the projectile escape first, the role of the fragments in the final water content of the class A HZ planets is less significant, while the contribution of the primordial content of water plays a primary role. Moreover, the final fractions of water of such planets show values close to the initial fraction of 10−4, which is comparable with those fractions obtained for the class A HZ planets formed in runs without fragmentation.

As for the class B HZ planets, the conclusions are very different. First, the physical properties of such planets formed in runs with and without fragmentation are very similar. In fact, the contribution of the generated fragments to the final final water content of the class B HZ planets is negligible regardless of the model adopted to distribute the water after a given collision. This can be attribute to the fact that there are no erosive impacts in water rich embryos, leaving water-rich fragments that can contribute significantly to the final water content of class B HZ planets. On the other hand, in regards to the final mass of the planets the generated fragments play a secondary role contributing in some cases up to 10% of the final mass of the planet. The most important result derived for this kind of planets indicates that the incorporation of the fragmentation in the collisional algorithm does not prevent the formation and survival of water worlds in the HZ of our systems.

The model that we adopted to carry out our investigation has some limitations that should be considered. First, the mass of each fragment has the same constant value in all our simulations. The selection of a lower value for such a mass would lead to more generated fragments per collision, and this could affect our results. In this sense, Wallace et al. (2017) analyzed the role of the fragmentation in the rocky planet formation at small distances from the central star, and investigated the dependence of their results on the choice of the mass of fragments in the collisional algorithm. While these authors did not find significant discrepancies in N-body simulations developed using different values for the fragment mass, we consider that it would be interesting to investigate the role of such a mass value in the final physical properties of the planets formed in the HZ of the system.

The models adopted from Marcus et al. (2010) to determine the water content of the final planets have to be looked at carefully. We use two models for mantle stripping in a differentiated body for modeling the water loss in a collision between two planetary embryos. Although this is true for water-rich bodies, we have to be careful with bodies with a low percentage of water. One caveat of this model is that, in reality, the water may not actually appear as an icy mantle, but rather inside the silicate mantle and atmosphere. In consequence, the model would no longer be valid and the results may change abruptly. Mantle stripping models that could be used for low or high water content have to be specified to achieve a good description of the final water contents. Creating more complex models for water loss by mantle stripping is in development to include this in future works.

On the other hand, we consider that the treatment of the smaller collisional fragments is a point of study to be checked in future works. In fact, such as described in Chambers (2013) and Wallace et al. (2017), our numerical algorithm conserves the total mass of the interacting bodies in a given collision, which is distributed between the largest remnant and the generated fragments. However, the model adopted by Mustill et al. (2018) incorporates a factor of mass removal assuming that most fragments are ground to smaller sizes in a collisional cascade, and then removed by radiation forces before they can be accreted on planetary embryos.

As for the switchover function K, Rein et al. (2019) investigated an inconsistency found in Chambers (1999) between the published function and that used in the available code. Rein et al. (2019) investigated various switching functions, in particular, two extreme cases: one infinitely differentiable function and the Heaviside function. Their results show that an infinitely differentiable switching function does not perform much better over one close encounter than the polynomial function used by MERCURY. These author concluded that, besides the discrepancy between the published switchover function, the MERCURY code is indeed a switching integrator.

Finally, we note that the final water contents of the HZ planets produced in our runs are derived by post-processing N-body simulations. Future works should analyze the evolution of the water content of such planets from N-body experiments that track the removal and transfer of volatiles after each collision. Prescriptions derived from smoothed particle hydrodynamics experiments, such as those shown by Dvorak et al. (2015) and Burger et al. (2018), should be included in future planetary formation N-body simulations to propose a more realistic model associated with the water evolution during the accretion stages.

The incorporation of such considerations in future N-body experiments will allow us to obtain a better understanding concerning the physical properties of the terrestrial-like planets that compose thewide diversity of planetary systems of the Universe.

Acknowledgements

We thank the anonymous referee for helpful comments that improved this paper. We want to thank Dr. Alexander J. Mustill for useful comments about the numerical recipes implemented in our code. This work was partially financed by Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) through PIP 0436/13, and Agencia de Promoción Científica, through PICT 2014-1292, PICT 201-0505 and PICT 2016-2635. Moreover, we acknowledge the financial support by Facultad de Ciencias Astronómicas y Geofísicas de la Universidad Nacional de La Plata (FCAGLP-UNLP) and Instituto de Astrofísica de La Plata (IALP) for extensive use of their computing facilities.

References

  1. Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241 [NASA ADS] [CrossRef] [Google Scholar]
  2. Asphaug, E., Agnor, C. B., & Williams, Q. 2006, Nature, 439, 155 [NASA ADS] [CrossRef] [Google Scholar]
  3. Burger, C., Maindl, T. I., & Schäfer, C. M. 2018, Celest. Mech. Dyn. Astron., 130, 2 [NASA ADS] [CrossRef] [Google Scholar]
  4. Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  5. Chambers, J. E. 2001, Icarus, 152, 205 [NASA ADS] [CrossRef] [Google Scholar]
  6. Chambers, J. E. 2013, Icarus, 224, 43 [NASA ADS] [CrossRef] [Google Scholar]
  7. Chambers, J. E., & Wetherill, G. W. 1998, Icarus, 136, 304 [NASA ADS] [CrossRef] [Google Scholar]
  8. Darriba, L. A., de Elía, G. C., Guilera, O. M., & Brunini, A. 2017, A&A, 607, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. de Elía, G. C., Guilera, O. M., & Brunini, A. 2013, A&A, 557, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Dugaro, A., de Elía, G. C., Brunini, A., & Guilera, O. M. 2016, A&A, 596, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Duncan, M. J., Levison, H. F., & Lee, M. H. 1998, AJ, 116, 2067 [NASA ADS] [CrossRef] [Google Scholar]
  12. Dvorak, R., Maindl, T. I., Burger, C., Schäfer, C., & Speith, R. 2015, Nonlinear Phenomena in Complex Systems (Switzerland: Springer Nature), 18, 310 [Google Scholar]
  13. Genda, H., Kokubo, E., & Ida, S. 2012, ApJ, 744, 137 [NASA ADS] [CrossRef] [Google Scholar]
  14. Hartmann, L. 1998, Accretion Processes in Star Formation (Cambridge: Cambridge University Press.) [Google Scholar]
  15. Hernandez, D. M. 2019, MNRAS, 486, 5231 [NASA ADS] [CrossRef] [Google Scholar]
  16. Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388 [NASA ADS] [CrossRef] [Google Scholar]
  17. Ida, S., & Makino, J. 1993, Icarus, 106, 210 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  18. Jacobson, S. A., Morbidelli, A., Raymond, S. N., et al. 2014, Nature, 508, 84 [NASA ADS] [CrossRef] [Google Scholar]
  19. Kokubo, E., & Genda, H. 2010, ApJ, 714, L21 [NASA ADS] [CrossRef] [Google Scholar]
  20. Kokubo, E., & Ida, S. 1996, Icarus, 123, 180 [NASA ADS] [CrossRef] [Google Scholar]
  21. Kokubo, E., & Ida, S. 1998, Icarus, 131, 171 [NASA ADS] [CrossRef] [Google Scholar]
  22. Kokubo, E., & Ida, S. 2000, Icarus, 143, 15 [NASA ADS] [CrossRef] [Google Scholar]
  23. Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131 [NASA ADS] [CrossRef] [Google Scholar]
  24. Krijt, S., Ormel, C. W., Dominik, C., & Tielens, A. G. G. M. 2015, A&A, 574, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Leinhardt, Z. M., & Stewart, S. T. 2012, ApJ, 745, 79 [NASA ADS] [CrossRef] [Google Scholar]
  26. Lodders, K. 2003, ApJ, 591, 1220 [NASA ADS] [CrossRef] [Google Scholar]
  27. Lodders, K., Palme, H., & Gail, H.-P. 2009, Landolt Börnstein (Berlin: Springer), 712 [Google Scholar]
  28. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
  29. Marcus, R. A., Sasselov, D., Stewart, S. T., & Hernquist, L. 2010, ApJ, 719, L45 [NASA ADS] [CrossRef] [Google Scholar]
  30. Morbidelli, A., Lunine, J., O’Brien, D., Raymond, S., & Walsh, K. 2012, Ann. Rev. Earth Planet. Sci., 40, 251 [NASA ADS] [CrossRef] [Google Scholar]
  31. Mustill, A. J., Davies, M. B., & Johansen, A. 2018, MNRAS, 478, 2896 [NASA ADS] [CrossRef] [Google Scholar]
  32. O’Brien, D. P., Morbidelli, A., & Levison, H. F. 2006, Icarus, 184, 39 [NASA ADS] [CrossRef] [Google Scholar]
  33. Quintana, E. V., Barclay, T., Borucki, W. J., Rowe, J. F., & Chambers, J. E. 2016, ApJ, 821, 126 [NASA ADS] [CrossRef] [Google Scholar]
  34. Raymond, S. N., & Izidoro, A. 2017, Icarus, 297, 134 [NASA ADS] [CrossRef] [Google Scholar]
  35. Raymond, S. N., O’Brien, D. P., Morbidelli, A., & Kaib, N. A. 2009, Icarus, 203, 644 [NASA ADS] [CrossRef] [Google Scholar]
  36. Rein, H., Hernandez, D. M., Tamayo, D., et al. 2019, MNRAS, 485, 5490 [NASA ADS] [CrossRef] [Google Scholar]
  37. Saha, P., & Tremaine, S. 1992, J. R. Astron. Soc. Can., 86, 291 [NASA ADS] [Google Scholar]
  38. Sánchez, M. B., de Elía, G. C., & Darriba, L. A. 2018, MNRAS, 481, 1281 [NASA ADS] [CrossRef] [Google Scholar]
  39. Springel, V. 2005, MNRAS, 364, 1105 [NASA ADS] [CrossRef] [Google Scholar]
  40. Stewart, S. T., & Leinhardt, Z. M. 2009, ApJ, 691, L133 [NASA ADS] [CrossRef] [Google Scholar]
  41. Testi, L., Natta, A., Scholz, A., et al. 2016, A&A, 593, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Wallace, J., Tremaine, S., & Chambers, J. 2017, AJ, 154, 175 [NASA ADS] [CrossRef] [Google Scholar]
  43. Weidenschilling, S. J. 1977, Ap&SS, 51, 153 [NASA ADS] [CrossRef] [Google Scholar]
  44. Wisdom, J., & Holman, M. 1991, AJ, 102, 1528 [NASA ADS] [CrossRef] [Google Scholar]
  45. Zain, P. S., de Elía, G. C., Ronco, M. P., & Guilera, O. M. 2018, A&A, 609, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1

Summary of amount of bodies formed for class A and class B HZ planets and their range of masses and final water fraction, in runs with and without fragmentation.

Table 2

Total amount and type of collisions ocurred for class A and class B HZ planets.

All Figures

thumbnail Fig. 1

Schematic representation of a collision between two planetary bodies. The target, with a mass Mt and a radius Rt collides with a projectile of mass mp and radius rp at an impactvelocity vi. The angle θ represents the angle formed between vi and the line connecting the two centers.

Open with DEXTER
In the text
thumbnail Fig. 2

Example of a mass distribution of planetary embryos from one simulation as a function of the initial semimajor axis at the end of the gaseous phase. The size of the points are scaled with the mass of each embryo: the different colors illustrate the initial fraction of water by mass for the embryos. In fact, orange (dark blue) indicates a fraction of 10−4 (0.5) water by mass. The sky blue shaded region represents the HZ for a solar-type star, while the dashed blue line illustrates the snow line assumed in our model.

Open with DEXTER
In the text
thumbnail Fig. 3

Histogram of percentage of collisions for different regimes implemented in the D3 N-body code for the scenario under study.

Open with DEXTER
In the text
thumbnail Fig. 4

Number of remaining bodies vs. integration time. The red curve corresponds to a simulation that lead to perfect mergers. The blue and cyan curves represent a simulation with the realistic collision treatment. In particular, the cyan curve only shows the decrease of planetary embryos. Both curves (red and blue) are representative for each set of simulations.

Open with DEXTER
In the text
thumbnail Fig. 5

Final outcomes at 200 Myr of 46 N-body simulations carried out with the D3 code. Left panel: 23 planetary systems that resulted from runs without fragmentation, while right panel: 23 systems produced from runs that incorporate fragmentation and hit-and-run collisions. The circle size scales with the planet mass according to that indicated in the top region of the figure. The sky blue area observed in both panels illustrates the HZ associated with a solar-type star (Kopparapu et al. 2013).

Open with DEXTER
In the text
thumbnail Fig. 6

Final mass as a function of the semimajor axis of the planets formed in simulations when fragmentation is activated (blue circles) and without fragmentation (red circles) fragmentation over 200 Myr of evolution.

Open with DEXTER
In the text
thumbnail Fig. 7

Final mass as a function of the semimajor axis of the class A (left panel) and class B (right panel) HZ planets produced in simulations with (blue circles) and without (red circles) fragmentation over 200 Myr of evolution.

Open with DEXTER
In the text
thumbnail Fig. 8

Distribution of the final mass for the surviving planets in the HZ. The violet bars represent the mass fraction corresponding to the initial mass of the body. The green and sky blue bars represent the mass fraction contribution from partial (or total) accretion of embryos and fragments, respectively. Top panels: the results for the class A HZ planets for Frag (left) and No Frag (right) simulations. Bottom panels: top panels but for the class B HZ planets. The height of each bar is normalized to the final mass of each planet.

Open with DEXTER
In the text
thumbnail Fig. 9

Final fraction of water by mass as a function of the final mass of the class A HZ planets (left panel) and class B HZ planets (right panel). The red filled circles illustrate the values associated with those planets formed in simulations without fragmentation, while the open and filled blue circles represent the values obtained from Models 1 and 2 in runs with fragmentation, respectively. In both panels, the black solid line represents the initial water fraction for class A and class B HZ planets.

Open with DEXTER
In the text
thumbnail Fig. 10

Same as Fig. 8 for the water contribution. The color code is the same as in Fig. 8. The panels represent the water contribution to the class A HZ planets using Model 1 (left panel) and 2 (right panel), respectively.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.