© ESO 2020

1 Introduction and motivation

The origin of water on Earth, but also on terrestrial planets in general, is a long-standing problem in planetary science, especially for planets with (current) orbital distances where protoplanetary disk material is believed to have been practically dry. All four terrestrial planets in the Solar System are examples thereof. Against those odds, Earth managed to acquire (and retain) a substantial amount of water during its formation, at least that stored in some form on the surface and in the crust, referred to as an Earth-ocean¹, and there is evidence for considerable water inventories on early Venus and Mars as well. In the case of Earth, even the current total water content is highly uncertain, due to suspected large reservoirs in the mantle, estimated to contain up to ten Earth-oceans (e.g., Marty 2012). Alternative theories, where growing planets inside the ice line (snow line) acquired their water content from local disk material via direct incorporation into growing planetesimals, suffer various difficulties (see, e.g., Izidoro et al. 2013; O’Brien et al. 2018; D’Angelo et al. 2019). The same is true for the hypothesis of delivery by comets, due to nonmatching isotopic fingerprints (D/H ratios), but more severely because the amount of water that can be dynamically delivered is significantly too low (Morbidelli et al. 2000). Therefore, dynamical delivery from an external source at the current location of the (outer) asteroid belt (or at least from the same source material) appears to be the most likely explanation at present. From the point of view of habitability, liquid water is believed to be one of the few absolutely necessary conditions for life. Excessive water content however can also have adverse effects for habitability (e.g., Noack et al. 2016), can impede sea-floor volcanism, and the formation of high-pressure ices at the bottom of global oceans may cut off connection to the underlying silicates. Therefore, an intermediate water content (like Earth) appears to be the most favorable for habitability. The increasing number of known low-mass exoplanets has unraveled a broad compositional diversity, ranging from probably iron-dominated, Mercury-like bodies to sub-Neptunes with extended envelopes. Proper characterization of their physical structure requires knowledge of their bulk composition, which is shaped by collisions on all scales, up to (and particularly) giant impacts between planet-sized bodies (e.g., Marcus et al. 2009, 2010; Inamdar & Schlichting 2016; Bonomo et al. 2019). For all these reasons it is important to understand the origin and evolution of water and other volatiles, which is intimately linked to their fate in collision events.

Several theories for the growth from dust to planetary embryos have been put forward, from the classical scenario of planetesimal- followed by oligarchic growth (e.g., Kokubo & Ida 2002; McNeil et al. 2005), to more recent alternatives like direct accretion from pebbles (see Johansen & Lambrechts 2017, for an overview). However, all these theories cumulate in a common final phase of accretion, where embryos and remaining planetesimals gradually grow into planets. This stage begins roughly when the gas disk disappears and the remaining smaller bodies cannot provide enough dynamical friction to keep embryos on separated orbits. It is characterized by chaotic interaction and giant collisions between embryos and/or protoplanets, which shape the final characteristics of planets. Indirect evidence for such events has been found in young extrasolar systems (e.g., Wyatt & Jackson 2016), mainly as observable infrared excess of warm dust attributed to debris of recent collisions (e.g., Meng et al. 2014; Morlok et al. 2014; Leinhardt et al. 2015; Bonomo et al. 2019). A direct consequence of late-stage accretion is radial mixing, which allows the inwards transport of water-rich material from beyond the ice line. For the Solar System, dynamical as well as isotopic evidence indicates that most of Earth’s water originates from the asteroid belt region, or at least from the same precursor material. This basic means of water delivery to growing terrestrial planets is also the underlying working hypothesis in this study. A potentially crucial influence comes from giant planets, which form quickly (within a few Myr), and can excite the orbits of planetary building blocks even if they remain on static orbits themselves and do not run intodynamical instabilities or migrate larger distances (as proposed in the framework of Grand Tack models for example; Walsh et al. 2011; O’Brien et al. 2014). It stands to reason that similar mechanisms of water transport are at work during the formation of many extrasolar systems despite large variations in individual histories.

Late-stage accretion is commonly studied by means of N-body simulations, where mostly the sole gravitational interaction of up to a few thousand bodies is followed for typically a few hundred Myr (e.g., Raymond et al. 2006; Kokubo et al. 2006; Morishima et al. 2010; Chambers 2013; Fischer & Ciesla 2014; O’Brien et al. 2006, 2014; Quintana & Lissauer 2014; Quintana et al. 2016, as some of the more recent work). Planets grow in principle by collisions, and at least for the fate of water, simple perfect merging has been and is still assumed almost exclusively (e.g., Haghighipour & Raymond 2007; Raymond et al. 2004, 2007; Izidoro et al. 2013; Fischer & Ciesla 2014; O’Brien et al. 2014; Quintana & Lissauer 2014), even though volatiles can be expected to be particularly affected by collisional processes. More sophisticated collision outcome models based on semi-analytical scaling laws have been developed, first distinguishing only perfect accretion and hit-and-run (Kokubo & Genda 2010; Genda et al. 2012), and later including a more extended range of outcomes (Leinhardt & Stewart 2012; Stewart & Leinhardt 2012; Leinhardt et al. 2015), and also an approximate treatment of compositional changes following Marcus et al. (2009, 2010). These models have been applied in various studies on planetary accretion (e.g., Chambers 2013; Bonsor et al. 2015; Carter et al. 2015, 2018; Quintana et al. 2016). However, they seem to be currently not suited to also account for minor constituents like surface volatile layers, or the more complicated behavior of hit-and-run with individual (water) losses and transfer between the colliding bodies, as tested by Burger et al. (2018). A very recent approach is the use of machine learning to predict collision outcomes, but this has not yet progressed beyond basic parameters like the mass accretion efficiency (Cambioni et al. 2019). A combination of dynamical and impact simulations to study the flux of water-rich material to habitable zone planets has been applied by Bancelin et al. (2017), albeit only at the end point of the dynamical evolution of inward-scattered objects. Direct coupling of N-body computations and individual collision simulations – the approach we use in this work – has rarely been applied so far (Genda et al. 2011, 2017), and only to follow nonvolatile material (iron and silicates).

While assuming complete accretion (perfect merging) of the colliding bodies is typically a good approximation for large-enough mass ratios, it is often not justified for more similar-sized bodies, with a diverse range of possible outcomes between (partial) accretion, hit-and-run, and (partial) erosion (see, e.g., Leinhardt & Stewart 2012). There have only been a few studies on water transfer and loss in collisions of large, gravity-dominated bodies. In most cases, only suites of individual, separate encounters were investigated, from planetesimal- to embryo-sized (Maindl et al. 2014, 2017; Dvorak et al. 2015), up to planet-sized (Canup & Pierazzo 2006; Burger & Schäfer 2017; Burger et al. 2018, 2019) and even in the super-Earth regime (Marcus et al. 2010). These studies showed that water losses in a single event can reach tens of percent, but often considered only combined losses and did not distinguish the fate of individual survivors. The relevant collision parameter space is large, consisting at least of impact velocity, impact angle, the mass ratio, and the combined mass, but also further ones, such as the composition of the colliding bodies. For this reason, such studies are naturally forced to fix one or more parameters. Unlike the other collision regimes, in hit-and-run there are two large post-collision bodies, where not only water losses but also transfer (such as from a water-rich impactor to a water-poor target) are important (Burger et al. 2018, 2019). Hit-and-run collisions can also be very transformative and strip especially the smaller body of its outer layers, such as atmosphere, ocean, or parts of the mantle (Asphaug 2010), or result in large-scale degassing of the interior due to pressure unloading (Asphaug et al. 2006). Since the bodies continue on separated trajectories after hit-and-run, it is imperative to track them both independently (Emsenhuber & Asphaug 2019), which rules out approaches where collisions are evaluated only as a post-processing step (e.g., after perfect merging runs).

Even though we focus on the collisional aspects of water transport, the physical state of mantle, surface, and atmosphere can all influencewater inventories and their evolution. Large collisions, in combination with heating by short-lived radioisotopes and differentiation, can lead to local or even global magma oceans. Existing water is either vaporized into a steam atmosphere, or gradually outgassed during magma ocean solidification, after which it is susceptible to atmosphericloss processes and impact erosion (Schlichting et al. 2015; Schlichting & Mukhopadhyay 2018, and references therein). Current models for coupled magma ocean and atmosphere evolution on approximately Earth-sized planets (Hamano et al. 2013; Lebrun et al. 2013; Salvador et al. 2017) suggest that beyond a critical distance (~0.7 au in the Solar System) rapid magma ocean solidification within a few Myr and subsequent water ocean formation are possible. Smaller bodies can experience considerable volatile losses also farther from the star before re-condensation (Odert et al. 2018). However, if most water is gradually transported inwards from larger orbital distances, the likelihood of retention also by smaller protoplanets is increased. With typical time-spans between giant collisions of several Myr (Sect. 6.5), repeated surface solidification and water ocean formation appear plausible, at least for already sufficiently large bodies in the habitable zone and beyond.

In the following we first describe the applied physical model including the initial conditions in Sect. 2, while aspects of its numerical realization are discussed separately in Sect. 3. An overview of the computed scenarios is provided in Sect. 4. The presentation of our results is split into two logical parts, where collisional water transport from a system-wide point of view is covered in Sect. 5, followed by a detailed analysis focused on potentially habitable planets in Sect. 6. A combined discussion in Sect. 7 and a summary in Sect. 8 conclude the paper.

2 The physical model

Our simulations focus on the gravitational and collisional dynamics of the late-stage accretion (giant impact) phase of terrestrial planet formation, and are set in an already gas-free environment. The initial conditions represent the transition from oligarchic growth to late-stage accretion, with a bimodal system of a relatively small number of approximately Mars-sized oligarchs (planetary embryos) embedded in a large number of smaller bodies (see, e.g., Kokubo & Ida 2002). Embryos are initially still on dynamically cold orbits due to dynamical friction with small bodies (in the following referred to as “planetesimals”) and possibly previous dampening by the gas disk, but become quickly excited by the influence of giant planets and by self-stirring. In general, all bodies interact gravitationally, except for planetesimals with each other. Planetesimal–planetesimal interactions can be either switched off entirely (also collisions), or they interact only in close encounters (gravitationally, and they can collide), which is further described in Sect. 3.1.

We include a 1 M_⊙ star, and already formed nonmigrating Jupiter- and Saturn-like gas giants at 5.2 and 9.6 au. These start on either circular orbits (similar to the Nice model, Tsiganis et al. 2005, but with their current semi-major axes), or eccentric orbits with today’s eccentricities of those planets (~0.05), and are referred to as “circular Jupiter and Saturn” (cJS) and “eccentric Jupiter and Saturn” (eJS), respectively. Similar configurations have been frequently used in the literature to study different giant planet architectures (e.g., O’Brien et al. 2006; Fischer & Ciesla 2014; Quintana & Lissauer 2014). In these works, eccentric gas giants have generally proved detrimental for water delivery (by dynamically removing water-rich material), at least in the absence of remaining gas (Morishima et al. 2010). We note that it is not our principal aim to study (or reproduce) the formation of the Solar System in detail, but rather to investigate the implications of realistic treatment of collisional water transport during late-stage accretion in general.

2.1 Initial conditions

Our main scenarios include an initial debris disk of approximately Mars-sized planetary embryos and approximately Moon-sized “planetesimals”, extending from 0.5 to 4 au. Both of these populations follow the same disk profile, with the surface density in solids $${\Sigma }_{\text{solid}}(r)={\Sigma }_0{\left(\frac{r}{1\hspace{0.17em}\text{au}}\right)}^{-\alpha }.$$(1)

We use Σ₀ = 10 g cm⁻² (1.125 × 10⁻⁶ M_⊙ au⁻²) and α = 1.5, consistent with the Minimum Mass Solar Nebula, in all scenarios, and do not consider discontinuities across the ice line. These parameters result in a total disk mass in solids of ~6 M_⊕.

The generation of planetary embryos follows the concept of the isolation mass M_iso = 2πabΣ_solid(a), which is simply the total mass available in some disk annulus with width b at semi-major axis a. The initial mass fraction in embryos is independent of the location in the disk and they are placed at separations of several times their mutual Hill radii $$R_{\text{H},\text{m}}=\frac{a_1+a_2}{2}{\left(\frac{m_1+m_2}{3\hspace{0.17em}{M}_0}\right)}^{1/3},$$(2)

with semi-major axes a₁ and a₂, masses m₁ and m₂, and the mass of the star M₀. The precise parameter setups for all our scenarios are described in Sect. 4. The above introduced disk structure (Eq. (1)) corrected for the actual mass fraction in embryos, along with the isolation mass and their mutual spacings, constrains the distribution of embryos to a single unique solution consisting typically of a few tens of bodies. For the chosen disk architecture the masses and separations of embryos both increase with semi-major axis. The remaining fraction of disk-solids is assumed to be still in equal-mass planetesimals that have not yet been accreted. The supposition that planetesimals follow the same surface density profile as embryos (Eq. (1)) results also in increasing separations with semi-major axis. Initial eccentricities and inclinations are drawn from Rayleigh distributions, with scale factors σ = 0.005 for embryos, and σ = 0.05 for planetesimals, meaning almost co-planar, low-eccentricity orbits. The three remaining orbital elements (argument of perihelion, longitude of ascending node, mean anomaly) are distributed uniformly over [0, 2π].

The embryos and planetesimals themselves generally have a differentiated three-layer structure with an iron core, a silicate mantle, and a water shell (cf. Fig. 1). Initially, all bodies are assigned an equal core mass fraction of 0.25, but water mass fractions (WMF) vary with distance to the star, following a step function: $$\text{WMF}=\{\begin{array}{cc}{10}^{-5}& \hspace{1em}\hspace{1em}a<2\hspace{0.17em}\text{au}\\ {10}^{-3}& \hspace{1em}\hspace{1em}2\hspace{0.17em}\text{au}<a<2.5\hspace{0.17em}\text{au}\\ 0.05& \hspace{1em}\hspace{1em}a>2.5\hspace{0.17em}\text{au}.\end{array}$$(3)

This particular distribution is based on meteorite data and has been regularly used in previous studies (e.g., Raymond et al. 2004, 2006, 2007; Haghighipour & Raymond 2007; Fischer & Ciesla 2014; Quintana & Lissauer 2014).

2.2 Collision model

Our collision model is based on dedicated hydrodynamical simulations among differentiated embryos and planetesimals. The bodies do not possess atmospheric components, and therefore the water inventory is represented by a surface liquid or solid layer, where all material is modeled as a fluid without material strength. We have chosen such properties not because of their simplicity (we could have included solid-body behavior quite easily and this usually is also not much more computationally expensive), but rather due to the expected broad diversity of material properties in the simulated environment. It is not only the current physical state (mass, composition, thermal state, and distance to the star), but especially also its history (thermal and orbital evolution, previous collisions) that would have to be modeled in order to reliably assess the current state of a given body. We believe that a fluid description is the best compromise, as long as those factors are not accounted for additionally.

In general, all embryos and planetesimals can collide once they come closer than the sum of their (uncompressed) radii, hence we consider only real (physical) collisions. In addition there is also the possibility to switch off interactions between planetesimals, meaning that they cannot collide or interact otherwise (see Sect. 3.1). Each collision results in either 0, 1, or 2 post-collision bodies, which allows us to capture all main collision outcome regimes (total disruption, accretion, erosion, hit-and-run), while still assuring that the number of bodies does not increase. The outcome is determined by computing the largest (most massive) and second-largest gravitationally bound aggregates of collision fragments (spatially connected material). These each remain in the simulation if their mass is above a lower threshold (0.5−1.5 × 10²¹ kg≅1−3 × 10⁻⁴ M_⊕, depending on the scenario), and are discarded when below. Gravitational aggregates above 1.5 ×10²³ kg≅2.5 × 10⁻² M_⊕ are treated as fully interacting embryos or protoplanets, and those below are included as planetesimals². The simulation then continues with the aggregates’ positions, velocities, masses, and compositions. Collisionally induced rotation is currently not tracked. It occasionally happens that two bodies evolve into a (mutually) bound configuration (if the apocenter distance of their relative orbit is smaller than their mutual Hill radius) but do not physically collide. In those situations the bodies are still assumed to accrete (computed via perfect merging) in order to avoid computational slow down by a long-term bound configuration. In collisions of protoplanets, embryos, or planetesimals with the star or one of the giant planets, the smaller body is completely accreted by the larger one. This happens in approaches closer than ~3 × R_G for the giant planets, and closer than ~5–10 × R_⊙ (depending on the scenario) for the star, with radii of the giant planet R_G and the Sun R_⊙. Due to numerical reasons, bodies are also merged into the star if their semi-major axis falls below a threshold of ~0.23 au.

In addition to this model we also performed simulations with simple perfect merging in all collisions for comparison. The physical model and the initial conditions are the same, except for the treatment (or the lack thereof) of collisions.

Fig. 1 Collision geometry of two differentiated bodies in a frame centered on the target. The impact angle α ranges from 0° for head-on to 90° for just grazing. The impact velocity is |v0|. The common convention is to label the larger body as the target, but in some situations the smaller object is the one followed (e.g., on its growth path towards a planet) and is thus considered (and referred to as) the target.

2.3 Ejections

In the dynamical environment of our scenarios a considerable number of bodies are ejected from the system, or at least evolve onto orbits with large perihelia where they very likely would no longer participate in terrestrial planet formation. We remove bodies from the simulation if they have distances to the star above 50 au and are on unbound (hyperbolic) orbits, and also those that are still bound but have perihelia above ~11 au.

3 The numerical hybrid model

Our simulations are based on three numerical pillars: long-term N-body integrations, smooth particle hydrodynamics (SPH) collision simulations, and the interface between them.

3.1 N-body numerics

For the N-body integrations we use the REBOUND package (Rein & Liu 2012), which offers several different integrators and a highly modularized design, which helped to incorporate it into our hybrid framework. The bulk of integrations are carried out with the symplectic Wisdom-Holman integrator WHFast (Rein & Tamayo 2015), using modified democratic heliocentric (WHDS) coordinates (Hernandez & Dehnen 2017), and a timestep of always less than 1 day. The embryo population is fully self-interacting, while two different models for the behavior of planetesimals are used:

• “pp0” – planetesimals never interact with each other (neither gravitationally nor collisionally), only with the other types of bodies;

• “pp1” – planetesimals interact gravitationally with each other in close encounters, but not elsewhere, and they are allowed to collide with each other.

Model pp1 is intended to allow small bodies some self-interaction (e.g., mutual dynamical friction), while it simultaneously provides gravitational focusing for more realistic planetesimal–planetesimal collisions. In order to accurately resolve close encounters, the integration switches to the IAS15 integrator (Rein & Spiegel 2015), a fifteenth order Gauß-Radau scheme, whenever two bodies come close. The changeover distance is chosen as ~2–3 × R_H,m for all encounters except those between planetesimals, where a somewhat lower value of ~2 × R_H,m is used, with the mutual Hill radius R_H,m (Eq. (2)). Whether protoplanets, embryos, and planetesimals actually (physically) collide is determined only once they approach each other to within approximately 10–20 times the sum of their radii, which has proven suitable such that their (theoretical)two-body orbit is no longer significantly disturbed. The pericenter distance of the (analytically computed) relative two-body orbit is then compared to the sum of the bodies’ uncompressed radii, which are also used to include potential “tidal collisions” (see, e.g., Asphaug 2010), even though they may still also collide physically due to tidal deformation. If on a colliding trajectory then the simulation is continued until they have approached to the desired starting distance of the SPH simulation, set to 2.5 times the sum of their radii, which still allows the approximate buildup of tidal deformation before contact. At this point the SPH simulation is started (see Sect. 3.2), while the whole N-body system is additionally integrated more accurately with IAS15 for its duration (with the colliding bodies replaced by a merged placeholder), before SPH results are reinserted and the N-body integration continues. Using this procedure allows us to precisely control the timing of the N-body integration and to have a well-defined starting distance for SPH.

We took great care to constantly monitor and check all conservation quantities and the consistency of results in all simulations. This is particularly important for the long-term N-body integration, including encounter handling, as well as the handover to and from the SPH simulations. As a good indicator for overall accuracy, the relative error in total energy is on the order of 10⁻⁵ for the whole evolution over several hundred Myr (naturally excluding collisions and ejections or removals, where for example mass is removed from the system).

3.2 Numerics of SPH collisions

Our hybrid approach requires a large number of dedicated hydrodynamical simulations performed with our SPH code miluphcuda. It allows for fast and efficient computation by utilizing highly parallel graphics processing unit (GPU) hardware and has already been applied in various studies of collision and impact physics, including nonviscous fluids (e.g., Burger et al. 2018; Malamud et al. 2018), solid bodies (e.g., Maindl et al. 2014; Burger & Schäfer 2017), granular media (e.g., Schäfer et al. 2017), and porous continua (e.g., Haghighipour et al. 2018). A comprehensive description of the code is given in Schäfer et al. (2016), and we provide only a brief summary of the most relevant features here.

The colliding bodies themselves are fully self-gravitating and in addition the star and all other bodies in the underlying N-body simulation are included as point masses to ensure a consistent environment. However, the kinematic state of the other bodies at the end of the SPH run is taken from the (usually more accurate) N-body integration which computes the whole system in parallel (see Sect. 3.1). SPH runs start once two bodies are found to be on collision trajectories and have approached to the desired initial distance (2.5 times the sum of their radii). The total simulation time is at least eight times the time until impact (roughly the collision timescale), which is intended to allow all main post-collision fragments to form and clearly separate. The thermodynamic response of simulated materials is modeled by the nonlinear Tillotson (1962) equation of state (EoS), which is applicable over a wide range of conditions, but is still relatively simple and efficient (see Melosh 1989, for a summary). The reason for not using a more sophisticated EoS is similar to the rationale for a general fluid description for all materials (see Sect. 2.2) – the rather large uncertainties in the actual physical state of bodies before the collision. In addition, it has been shown that especially for lower-density regions very high particle numbers are required for accurate results, where outcomes often depend on resolution rather than precise modeling by the EoS (Reinhardt & Stadel 2017). We use Tillotson parameters for iron from Melosh (1989) and for basalt and water(-ice) from Benz & Asphaug (1999).

The differentiated bodies (iron core, silicate mantle, and water shell) are set up in an already relaxed state by using self-consistent hydrostatic radial profiles. Those can be computed very rapidly and eliminate the need for lengthy numerical relaxation. We refer to Burger et al. (2018) for details on this procedure. A numerical issue that requires special care is the spatial extension (thickness) of water layers. Water mass fractions in the initial conditions range from 10⁻⁵ to 0.05 (Sect. 2.1) and can and do naturally decrease further in collisions. We apply a lower WMF limit of 10⁻⁹, below which objects are considered practically dry and their WMF is set to zero. However, actual water layers are regularly thin, and would be made up of only a few SPH particles. In order to maintain a reasonable extent in terms of SPH numerics and resolution, we artificially increase simulated water layers to a thickness of at least approximately three SPH particles where necessary (by replacing mantle mass by water). This increase is corrected for after the simulation based on the assumption that relative water loss and transfer (e.g., in hit-and-run) is approximately independent of the actual thickness of the water layer. Therefore, putting α times as much water mass on a body means that about α times as much (as would actually have been) is lost or transferred. This approximate scale-invariance has been tested and confirmed by Burger et al. (2018). The underlying physical reasoning is that if some part of a body’s water shell is affected (lost, transferred, etc.) then the fate of this part is relatively independent of its thickness, where the error naturally decreases towards thinner surface layers. The respective correction factor (1/α) can be directly applied to the WMF of post-collision bodies, because the WMF is directly proportional to the (absolute) water content (for a fixed total mass). Since the colliding bodies can have vastly different WMFs (either prior to the artificial increase, after it, or both), it is necessary to additionally distinguish between projectile and target water and to keep track of the origin of SPH particles. However, for these individual reservoirs the above considerations hold again (Burger et al. 2018). This procedure leads to typical WMFs in SPH runs of around 0.05–0.2.

The post-collision analysis of identifying gravitationally bound objects and their dynamics and composition comprises two main steps. First, a friends-of-friends algorithm identifies all directly (spatially) connected clumps of SPH particles (we refer to them as “fragments”), before gravitationally bound “aggregates” of those fragments are identified in the second step. Starting from a seed aggregate consisting of the most massive fragment and all those mutually gravitationally bound to it, an iterative algorithm adds/removes fragments to/from this seed aggregate (its barycenter) until convergence. After identification of the largest aggregate, this algorithm is run once again to find the second-largest one, starting from the most massive fragment not yet incorporated into the largest aggregate. We refer to Burger et al. (2018) for more details on this procedure. For hit-and-run, the largest aggregate usually originates in the target (the more massive pre-collision body) and the second-largest one in the projectile. Each of the two aggregates is then either included as embryo (protoplanet), as planetesimal, or is otherwise discarded (see Sect. 2.2). If kept, the N-body run continues with the aggregate’s barycenter position and velocity.

All SPH simulations use between 25 000 and 75 000 SPH particles, depending on the mass ratio of the colliding bodies, with 25 k particles for mass ratios down to 1:25, 50 k between 1:25 and 1:100, and 75 k below 1:100. This ensures a reasonable resolution even for scenarios involving large targets and comparatively small projectiles, while saving computing time in more similar-sized collisions. The runtime depends mainly on particle number and impact velocity (which determines the overall simulated time), where a single SPH simulation, including all pre- and post-processing, takes on the order of an hour to complete.

4 The scenarios

Our scenarios, all summarized in Table 1, consist in general of a disk of embryos and planetesimals in the terrestrial planet region (0.5–4 au), embedded in a Solar System-like architecture, with a 1 M_⊙ star and two giant planets resembling Jupiter and Saturn. The latter are either on the current orbits of those planets with today’s eccentricities (eJS) or initialized on circular orbits (cJS). The surface density profile, total mass, and extent of the embryo + planetesimal disk are in principle the same in all simulations (see Sect. 2.1).

For our regular scenarios, the mass fraction of the disk initially in embryos is set to 0.7, leaving 0.3 in planetesimals. With an embryo spacing of 10 mutual Hill radii, this gives a system of 35 embryos (between ~0.06 and 0.25 M_⊕), and 150 equal-mass (~0.012 M_⊕ ~ M_Moon) planetesimals with approximately 1/10 of the average embryo mass. These parameters provide a significant fraction of the total mass initially still in planetesimals combined with planetesimals that are considerably smaller than embryos in order to have two distinct populations, while also maintaining acceptable runtimes. In addition to the two giant planet architectures (cJS and eJS), the two different models for planetesimal–planetesimal interactions (pp1 and pp0, Sect. 3.1) were applied. Due to the extremely stochastic nature of late-stage accretion, we performed five simulations with the SPH collision model for each of these parameter combinations, and three additional ones with perfect merging (PM), for comparison. This gives a total of 32 simulations for the regular scenarios (upper part of Table 1).

In addition, we also consider scenarios without planetesimals, that is, starting with a pure embryo disk, to study the principal influence of the planetesimal population on the dynamical as well as the collisional evolution. With the same disk parameters this results in 29 embryos, with masses between ~0.1 and 0.4 M_⊕. While for PM runs, planetesimals cannot form at all now, this is still possible with our collision model, but since usually a few planetesimals are around at most we do not expect a significant difference between pp0 and pp1. Therefore, we allow such interactions (pp1) in all embryo-only scenarios. As for the regular runs we performed five simulations with the SPH collision model and three with PM, for both the cJS and eJS architectures, which gives a total of 16 scenarios here (lower part of Table 1). We note that the embryo-only initial conditions are similar to those in the few existing and comparable direct hybrid simulations (Genda et al. 2011, 2017), where our disk spans a wider range in semi-major axis and thus includes more embryos than in these studies.

The total disk mass in the regular runs amounts to 6.1 M_⊕, including a total of 0.1 M_⊕ of water, equal to 345 Earth-oceans³ (O_⊕). Except for ~0.1 O_⊕ inside 2 au, all remaining water is located in the outer disk initially. These figures are in general the same for the embryo-only scenarios but can differ somewhat due to the initial configurations’ setup algorithm (see Sect. 2.1), with a slightly smaller disk mass (5.9 M_⊕), and a somewhat higher initial water inventory (375 O_⊕). The latter difference is found almost only in the outer disk beyond 2 au. Each scenario was computed for at least 200 and at most 700 Myr – until a system of final planets on approximately stable orbits had formed. The runtime per simulation with regular (embryos + planetesimals) initial conditions was on the order of 2–4 months, where those with the SPH collision model took at least twice as long as with PM. The main reason for that is not directly the extra computing time of the SPH runs (on the order of a week of combined GPU time), but the considerably slower decrease in the number of bodies in the N-body simulation (see Fig. 9). In each of the main scenarios about 150–200 collisions were treated with SPH, in addition to approximately 50 PM collisions, mainly with the star and the inner giant planet. For the embryo-only scenarios these figures are about 40–50 SPH-treated collisions, in addition to about 10 with PM. Therefore, altogether we have simulated around 4000 collision scenarios with dedicated SPH simulations.

5 Results I: system-wide evolution and water transport

In this first part we primarily present results from a system-wide point of view, while the accretion histories and water delivery to potentially habitable planets are presented in more detail in Sect. 6. An overview of all scenarios and their proper names used throughout this paper can be found in Table 1. Since this is the first time a hybrid code with dedicated collision simulations is applied to debris disks with embryos and an additional planetesimal population, we ran scenarios including planetesimal–planetesimal interactions (pp1) and also without them (pp0) in order to compare this novel aspect. However, for the most part, the outcomes of pp0 and pp1 are fairly similar, particularly compared to variations among other parameters (further discussed in Sect. 7.4). Therefore, we often combine pp0 and pp1 results in the following, unless stating otherwise. A similar policy is applied to the embryo-only simulations, which are only referred to when relevant to highlight effects of omitting the initial planetesimal population, while we focus on the (probably more realistic) regular scenarios otherwise (but see Sect. 7.3 for discussion). Results for embryo-only scenarios are includedin Tables 2–4, and in Fig. 3, but not in the other figures.

For quantifying water contents we mostly use absolute values (in easy-to-interpret units of O_⊕) instead of WMFs, with the exceptions of Figs. 2, 3, and 6 because using WMFs is common in this kind of plot in the literature. Absolute amounts, for example of water delivered to potentially habitable planets, are useful in relation to the known initially available supply, but also because the thickness of (surface) water layers is usually better indicated by the absolute inventory⁴. In addition, the limitations of our collision model suggest that some neglected refractory debris is likely re-accreted at some point, which adds mass and could therefore change the WMF but not directly the absolute water mass.

To enhance clarity we apply a uniform color coding and symbols for the different parameter combinations in all plots throughout the paper. Results from the SPH collision model are always shown as crosses or solid lines, and those from PM runs are indicated by circles or dashed lines. cJS scenarios are color-coded in black (pp1) and dark-green (pp0), and eJS simulations are shown in red (pp1) and orange (pp0).

Fig. 2 Final systems for all regular scenarios (cf. Table 1). Each column represents a particular giant planet architecture (cJS or eJS) and mode of planetesimal–planetesimal interaction (pp1 or pp0). Upper five rows: results based on the SPH collision model, lower three: PM comparison runs. The sizes of planets and their iron cores (black) are proportional to mass1∕3. The region corresponding to our definition of potentially habitable planets (see Sect. 6.1) is indicated as a gray bar, and the Solar System planets are plotted as open circles as a reference.

Fig. 3 Final systems for all embryo-only scenarios (cf. Table 1). The columns represent cJS and eJS architectures, while planetesimal–planetesimal interactions are always included (pp1). Upper five rows: results based on the SPH collision model, lower three: from PM comparison runs. Colors, symbols, and their sizes have the same meanings as in Fig. 2.

5.1 Final terrestrial planet systems

The final planetary systems are illustrated in the a–e plane in Figs. 2 and 3 for the regular and embryo-onlyscenarios, respectively. The upper five rows show results from our SPH collision model (mostly referred to simply as “our model” in the following), with the perfect merging (PM) comparison runs below. The main properties of formed planets as a function of semi-major axis are plotted in Fig. 4, and basic data and statistics on all scenarios are listed in Table 2. The final planets’ orbital parameters are always averaged values over the last few ten Myr to smooth out naturally occurring fluctuations (oscillations) caused by mutual perturbations. One of the first things that comes to mind when looking at Figs. 2 and 3 is the well-known highly stochastic nature of late-stage accretion, irrespective of whether our model or simple PM is used. The sample size of 5 + 3 scenarios (our model + PM) per parameter set (imposed by computational limits) is naturally insufficient to derive robust statistics, but general trends, and an assessment of the degree of scatter can nevertheless be obtained. Even though this is not our main intention, it is noteworthy that the formed systems rarely resemble the inner Solar System as a whole, but individual planets do, a well-known characteristic of the stochastic formation process.

Between one and four planets are formed in all scenarios, where between zero and two of them (one in the majority of runs) follow our definition of potentially habitable, that is they end up between 0.75 and 1.5 au (see Sect. 6). We count leftover embryos as planets, but not leftover planetesimals. Most systems end up with either two or three planets – and only rarely four – which is somewhat lower than what was often found in previous studies where roughly three to four planets were formed. A principal factor could be our relatively massive debris disk, with Σ₀ = 10 g cm⁻² and extending from 0.5–4 au (otherwise often Σ₀ = 8 g cm⁻² and/or a cutoff at 2 au), where it has been shown that lighter disks typically produce more planets (Kokubo et al. 2006; Izidoro et al. 2013). Besides that, a better-resolved population of small bodies (planetesimals) can provide more efficient eccentricity damping by dynamical friction (O’Brien et al. 2006; Raymond et al. 2006; Chambers 2013), generally leading to a higher number of planets (Chambers 2001; Quintana et al. 2016). With our collision model a relatively large amount of small-scale collisional debris is often neglected, which would otherwise also enhance dynamical friction and lower eccentricities. In contrast to that are results from Chambers (2013) and Quintana et al. (2016), who use the semi-analytical collision outcome model by Leinhardt & Stewart (2012); Leinhardt et al. (2015), which includes collisional debris as additional bodies (once sufficiently massive) that can contribute to dynamical friction. With this model these latter authors obtain a considerably higher number of planets (up to 5) than with PM, unlike in our results. A potential further influence in simulations with our model might be the substantially lengthened accretion times (see Sect. 5.4), which could lead to more ejections and removals of bodies due to longer chaotic interaction. From the current data it is not clear which of these mechanisms is mainly responsible for the discrepancy in the number of formed planets and which results are more realistic.

The masses of the planets formed with our model are significantly lower than with PM (Table 2), which is also clearly visible in the first panel of Fig. 4. The latter plot also illustrates the strong decline of mass beyond ~1.5 au (presumably due to the disturbing influence of the giant planets), consistent with earlier studies (e.g., Chambers 2001, 2013; Fischer & Ciesla 2014; Genda et al. 2017), and the high number of planets formed around 0.5 au, the inner edge of the initial disk. These planets end up much less massive with our model than with simple PM, which is also clearly visible in Fig. 2 and the first panel of Fig. 7, and often have highly elevated core mass fractions (CMF), as evident from the fourth panel of Fig. 4. This has to be viewed in light of our numerical model however, where small-scale debris and its potential re-accretion is not further tracked. This interesting outcome, possibly a robust way to produce Mercury-like planets, is further discussed in Sect. 7.1. Most other planets show only a moderate increase in CMF to around 0.3, consistent with earlier studies including fragmentation but complete mass conservation (e.g., Chambers 2013).

The mass and composition of the final systems generally depend on the initial inventory, which is reduced by material lost to the star or to the giant planets, removed or ejected from the system, or neglected as debris after SPH simulations. With our model, from the initial ~6 M_⊕, approximately 2 M_⊕ are accreted into planets with cJS versus ~1.5 M_⊕ with eJS. In the PM comparison runs these figures are ~3.5 M_⊕ for cJS versus ~3 M_⊕ for eJS (Table 2). The decisive influence of the giant planet architecture has been studied comprehensively, mostly with the classical static orbits (e.g., Raymond et al. 2004, 2006, 2007; O’Brien et al. 2006; Morishima et al. 2010; Izidoro et al. 2013; Fischer & Ciesla 2014; Quintana & Lissauer 2014), but also for migrating gas giants as in Grand-Tack scenarios (e.g., O’Brien et al. 2014). Our results confirm the detrimental influence of eccentric giant planets on the total surviving mass, also with our model. On the other hand, the large variations in final water content found in earlier studies can be reproduced with PM, but are considerably relaxed with our model. In the PM runs, the amount of surviving water (system-wide) is at most a few tens of O_⊕ for eJS, but ~100–150 O_⊕ with cJS, while with our model this figure is surprisingly similar for eJS and cJS with average values of ~30–40 O_⊕ (cf. Table 2 and Fig. 5), which is further discussed in Sect. 5.3.

Fig. 4 Various properties of final planets in the regular scenarios, as a function of semi-major axis. Mass and core mass fraction (first and last panel) are plotted linearly, while water mass fraction and total water mass (middle panels) are log-plots. The region of potentially habitable planets (see Sect. 6.1) is indicated as a gray bar.

5.2 Radial mixing and fate of available building material

The time evolution of the available building material (~6.1 M_⊕) and total water inventory (~345 O_⊕) is plotted in Fig. 5. The total mass as well as water inventory decrease much more rapidly for eJS (red + orange) than for cJS (black + dark-green) for several Myr, which is also clearly visible in Fig. 6. However, with our model especially the water contents converge to similar values eventually, further elaborated in Sect. 5.3. To also track radial mixing, the time evolution of the mass in four semi-major axis bins is plotted in Fig. 7, and correspondingly for water in Fig. 8. The second panels in these plots outline the region where potentially habitable planets might form. Differences between cJS and eJS at the very beginning are clearly visible in the fourth panels, where strong dynamical forcing of the giant planets with eJS (red + orange) immediately (~5–10 Myr) removes the majority of the material from that region (cf. Fig. 6). The disturbing influence of the giant planets in our scenarios is mainly due to the 3:1 and 5:2 mean-motion resonances and the ν₆ with the Saturn-like planet with eJS, which force bodies in the outer disk onto high-eccentricity orbits and often drive them into the star, a giant planet, or lead to ejection. In the cJS simulations (black + dark-green) dynamical forcing acts as well, but is less pronounced and sets in more slowly, therefore considerably more water typically survives beyond 2.5 au and slowly diffuses inwards. With eJS, a general decline of material also between 1.5 and 2.5 au can be observed within the first 10–20 Myr, while with cJS it remains rather constant or even increases initially. Radial mixing of water happens in all scenarios, and takes between a few and ~10 Myr to reach the 1.5–2.5 au region, and on the order of 10–30 Myr to arrive between 0.75 and 1.5 au. In earlier studies, similar timescales of around 10Myr for the onset of water accretion onto inner planets have been found (Raymond et al. 2004, 2006; Fischer & Ciesla 2014). The later evolution of both mass and water content is generally marked by a slow and steady decline caused by further losses of bodies from the system, or collisional losses in SPH simulations for scenarios with our model.

The amount of material lost outwards, either by ejection or removal, or collision with a giant planet, is remarkably constant throughout all scenarios between 2.2 and 2.6 M_⊕ (Table 2), which is more than a third of all available material. Comparable results were found by O’Brien et al. (2006) with strong dynamical friction, and Izidoro et al. (2013) reported losses of more than 50% in a similar dynamical environment (but perhaps including collisions with the star). For water it is even the majority of the initial inventory that is typically lost to ejections, removals, or the giant planets, about 250–300 O_⊕ in the regular scenarios. Those losses are only considerably lower in the PM-cJS runs with ~200 O_⊕, presumably due to the combination of less dynamical scattering and significantly faster accretion times (see Sect. 5.4), which prohibit further gradual losses (compared to SPH-cJS). Collisions with the star account for about 0.5–1 M_⊕ with eJS (PM and our model alike), but only about a third of that with cJS, consistent with earlier studies (e.g., Raymond et al. 2004; O’Brien et al. 2006). The difference is even larger for water, with an average of several tens of O_⊕ with eJS, but typically only around 1/10 of that with cJS, presumably also caused by strong eccentricity forcing of outer disk material with eJS.

A third source of material loss is due to our collision model (Sect. 2.2), where material that is not included in the further tracked bodies (gravitationally bound aggregates) is neglected. With cumulative amounts of between ~1 and 1.5 M_⊕ (see Table 2) these “losses” are surprisingly high, and while a considerable fraction of the refractory component may be re-accreted at some later point if it were tracked, the majority of the water content is likely lost. Similar values have been found in comparable hybrid simulations by Genda et al. (2017). The background and rationale for ignoring this debris is that we are foremost interested in the fate of water, where small-scale fragments are probably often devolatilized to a large degree anyway after going through a high-energy collision. The iron mass fraction of the neglected debris is approximately constant and around ~0.1, clearly confirming that collisional erosion strips preferentially outer layers, while core material is lost only in the most destructive events. The total amount of water lost in SPH collisions is on average tens of O_⊕ with cJS, but several times less with eJS, probably because there is simply less water included in eJS-scenario collisions due to the high early losses of water-rich material (cf. Fig. 5). These higher collisional losses with cJS contribute significantly to eventually reduce the system-wide water content to similar levels as in eJS runs (see Sect. 5.3).

Fig. 5 Time evolution of the total mass (upper panel) and total water mass (lower panel) for all regular scenarios.

5.3 Water inventories of final planetary systems

The water contents of all formed planets are in principle illustrated in Figs. 2 and 3 via their WMF. The absolute amount of water (in O_⊕) is plotted in the third panel of Fig. 4 and is given in Table 2. Here we discuss global water inventories of the formed systems, while potentially habitable planets in particular are addressed in Sect. 6.

All scenarios begin with a total of 345–375 O_⊕ (~0.1 M_⊕), located almost exclusively beyond 2 au. In the regular scenarios and with our collision model an average of ~30–40 O_⊕ survivesthe accretion phase (system-wide), similar for cJS and eJS alike (cf. Table 2 and second panel in Fig. 5). Considering only the region of potentially habitable planets (0.75–1.5 au), average water contents with cJS are still only about a factor of two above eJS. The lack of a larger difference is surprising and unexpected at first sight. Previous studies with similar dynamical environments and initial conditions, but without any dedicated collision model (hence assuming PM), have consistently found that with giant planets on circular orbits (akin to cJS), planets emerge much more water-rich than in eJS-like scenarios (e.g., Raymond et al. 2004; O’Brien et al. 2006; Quintana & Lissauer 2014; Fischer & Ciesla 2014). Our PM runs reflect these results perfectly (e.g., Fig. 4), while the similar water contents resulting from a realistic collisional evolution appear to be an interplay of different loss mechanisms and presumably also lengthened accretion, as elaborated below. However, the scatter in system-wide water content is much larger for eJS, ranging between ~0.5 and 98 O_⊕, but only between ~7 and 70 O_⊕ with cJS. Strong resonant forcing by the giant planets, mainly of water-rich bodies in the outer disk, acts particularly with eJS (see Sect. 5.2). Nevertheless, the bulk of this effect is not to deliver water-rich material to the inner disk, but to efficiently (within ~5–10 Myr) remove it from the system. The consequences are that only a few water-rich bodies take part in the further evolution with eJS, which is clearly visible in Fig. 6 (left column), where especially the individual fate of water-rich embryos with their large inventories⁵ can be decisive. This leads to high stochasticity and thus the large scatter in final water contents with eJS, and even more so for the subset of potentially habitable planets (Fig. 8 and Table 2). With cJS, dynamical removal of outer disk material is much weaker and leaves more water-rich bodies to diffuse inwards (Fig. 6, middle and right columns), making water delivery more robust. One manifestation of these differences between cJS and eJS are water-rich embryos that remain stranded in the outer disk (see scenarios SPH01, SPH17, SPH18, PM07 and PM12 in Fig. 2). Whether an individual embryo is transported inwards and delivers water to growing planets, survives as a stranded embryo, or is removed from the system often makes a tremendous difference with eJS, but usually has much less of an effect in cJS runs. Another example are the highly differing (average) water contents of systems formed in the PM-eJS scenarios (PM07 to PM12, see Table 2). Closer inspection revealed that the large apparent difference between the respective pp1 and pp0 scenarios is merely due to the high stochasticity of water delivery with eJS, without an actual systematic difference. Therefore, an average over all six scenarios would be more appropriate (and less misleading).

Initial dynamical removal of large quantities of water-rich material happens in all scenarios, and is particularly strong with eJS. In PM simulations this point seems to be the decisive split between cJS and eJS, where for cJS accretion into typically large, water-rich planets proceeds quickly and prohibits much further losses, while for eJS only so much water-rich material is left to be distributed among the forming planets (cf. Figs. 2 and 3). While the initial removal phase is similar to PM, realistic treatment of collisions changes the subsequent evolution drastically, where several tens of O_⊕ are removed in the form of continuous collisional losses with cJS, but only a fraction of it with eJS (Table 2), presumably simply because less water is available (after initial removal). The considerably prolonged growth phase compared to PM (Sect. 5.4) can be assumed to also allow for more losses of water-rich material. This can readily explain the increased amount of outward losses (collisions with the giant planets or ejections and removals) with SPH-cJS, which is much closer to the high values of SPH-eJS runs than to PM-cJS scenarios (Table 2).

A general observation is that WMFs never increase above the initial maximum value of 0.05 (second panel in Fig. 4), confirming earlier results (Marcus et al. 2010; Burger et al. 2018). A decrease below the initial minimum of 10⁻⁵ is of course possible (with our model), but actually happens only in relatively few cases, presumably also because at least some water-rich material is frequently accreted at some point and not necessarily due to inefficient collisional water losses. Another interesting outcome is that absolute water masses are rather unaffected by semi-major axis, except for the innermost part of the disk around 0.5 au (third panel in Fig. 4). This indicates efficient mixing of water also into the region of potentially habitable planets (see Sect. 5.2). For the WMFs of the final planets on the other hand, the situation is more complicated (second panel in Fig. 4), where the decreasing masses of planets formed in the outer disk lead to higher WMFs, also for similar absolute inventories (see also, e.g., Raymond et al. 2007). This is where one has to be careful with notions of “water-rich”, “water-poor”, and so on, as discussed at the beginning of Sect. 5.

Fig. 6 Simulation snapshots in the a–e plane for three different runs (columns). Typical evolution with eJS in the left column (SPH15-eJS-pp1) vs. cJS in the middle (SPH06-cJS-pp0) vs. cJS with PM on the right (PM06-cJS-pp0). Symbol sizes are proportional to mass1∕3. Each row represents an equal time in all scenarios: up to 1 Myr the evolution appears similar, except some more dynamical excitation with eJS (left). Within only ~5 Myr however, the great majority of water-rich material is removed with eJS, but remains in the system with cJS and starts to diffuse inwards. At 50 Myr, accretion is still actively ongoing with the SPH collision model (even though less so with eJS in this particular scenario), while with PM (right) it is almost over, except for some remaining planetesimals. We note that the cJS scenarios (middle and right columns) are pp0 runs, where some more leftover planetesimals are common. After ~100 Myr, body numbers are declining in the SPH scenarios (left and middle columns), but accretion and associated collisional evolution is still not over. We also note the damping effect on the eccentricity of the giant planet with eJS, consistent with similar studies (e.g., O’Brien et al. 2006). These scenarios represent typical behavior and were not specifically selected, except that the left and middle runs are those where the potentially habitable planets in Figs. 18 and 17 form.

Fig. 7 Time evolution of the mass in four semi-major axis bins for all regular scenarios. The curves have been time-averaged (smoothed) over ~3 Myr.

Fig. 8 Time evolution of the water mass (in O⊕) in four semi-major axis bins for all regular scenarios. We note the log y-axis in the upper two panels. The curves have been time-averaged (smoothed) over ~3 Myr.

Fig. 9 Total number of bodies (embryos or protoplanets + planetesimals) over time for all regular scenarios.

5.4 Length of the accretion phase

Our model allows for 0, 1, or 2 further-tracked post-collision bodies, therefore their number is still monotonically declining, but considerablyslower than with PM. This suggests longer accretion times in general, and indeed, when the time of the last collision (excluding star and giant planets) is used as a measure, accretion times roughly double compared to PM (Table 2). Interestingly, the same holds if the last collision only among embryos and protoplanets is considered, which is probably a more robust measure since it excludes the (in principle indefinite) accretionary tail and rather signifies the last giant collision. The mean time of the last collision only among embryos and protoplanets is between 80 and 100 Myr in PM runs, but between 150 and well over 200 Myr with our model. These differences are even more pronounced in the embryo-only scenarios, presumably due to the absence of damping by dynamical friction, which results in more energetic collisions with a lower likelihood of accretion. Looking at the growth curves of individual (potentially habitable) planets (Fig. 15, lower panel) confirms that planets formed with PM finish accretion considerably sooner than those assembled with our model. Another aspect is the frequency of collisions, which is naturally declining over time and will be further discussed in Sect. 6.5.

Figure 9 shows the time evolution of the total number of bodies, where accretion times with our model seem to be clearly lengthened and body numbers are roughly twice those in PM runs over long periods of time. Similar results have also been found in previous studies where collisional fragmentation was included (Chambers 2013; Quintana et al. 2016), based on the semi-analytical model by Leinhardt & Stewart (2012); Leinhardt et al. (2015). However, in these studies the number of embryos (protoplanets) declined approximately equally, and the authors concluded that hit-and-run and fragmentation do not actually lengthen the bulk growth of planets, but that higher body numbers are mostly due to low-mass planetesimals (collision fragments). Our results indicate a different behavior. Figure 10 shows that the number of embryos + protoplanets (i.e., excluding planetesimals) is also significantly higher compared to PM for most of the time, and a similar plot (not included here) for the embryo-onlyscenarios shows the same behavior. The growth times for potentially habitable planets to reach 50% (t₅₀) and 90% of their final mass (t₉₀), summarized in Table 3, also confirm that our collision model indeed results in roughly doubled accretion times compared to simple PM. Chambers (2013) on the other hand found similar values for t₅₀ when comparing their fragmentation model to PM, possibly because of differences of the applied collision and fragmentation models. While we include two post-collision bodies at most and neglect the remaining material not gravitationally bound to them, Chambers (2013) and Quintana et al. (2016) conserve mass by adding additional fragments on (slightly) escaping trajectories if necessary. The additional dynamical friction of those fragments might lead to more gentle subsequent collisions and therefore explain their faster decline in embryo + protoplanet numbers (see also O’Brien et al. 2006). However, the lowest allowed planetesimal/fragment mass in their model is ~ 5× 10⁻³ M_⊕ (for lower masses they simply assume PM), while in our model this is ~ 2× 10⁻⁴ M_⊕, which allows fragment formation in a wide mass range where theirs already assumes PM. The reason for their faster decline in embryo + protoplanet numbers may also be found somewhere else in the differing methodology, where a priori our actual hydrodynamic simulations for each individual collision scenario are probably more reliable than semi-analytical scaling laws, for example in predicting hit-and-run outcomes. Eventually it remains unclear which results are more accurate, and further work is required to gain a better understanding of this important issue.

Fig. 10 Number of embryos + protoplanets (i.e., excluding planetesimals) over time for all regular scenarios, split into cJS (upper panel) and eJS (lower panel) for better visibility. The respective “PM” curves have a slight y-offset of −0.2 for the samereason.

5.5 Global collision characteristics

The regular scenarios each include about 150–200 collisions among embryos, protoplanets, and planetesimals with our model (30–45 among embryos and protoplanets alone), and roughly 100 collisions in the PM runs (about 40–50 in the embryo-only scenarios, and only ~15–20 in their PM versions). These numbers are relatively robust and show only little variation on average (even when the different planetesimal–planetesimal interaction models pp1 and pp0 are compared). A detailed list of collision number statistics specifically for the growth histories of potentially habitable planets is included in Table 4 and discussed in Sect. 6.4.

In the regular scenarios about 100 collisions (out of the 150–200) result in hit-and-run, more than half of all events. In the embryo-only runs even more than 75% are hit-and-run, presumably because this outcome is more likely among more similar-sized bodies, and due to the lack of dynamical friction (to reduce impact velocities). The latter is also supported by the fact that in the regular runs, where dynamical friction is at work, the hit-and-run rate also in collisions among embryos and protoplanets alone is not substantially above 50%. Studies in similar dynamical environments found results in broad agreement: Kokubo & Genda (2010) reported 49% hit-and-run (embryo-only simulations), Chambers (2013) found that 42% of collisions among embryos and protoplanets are hit-and-run (they also had planetesimals), Quintana et al. (2016) found that 31% of impacts that lead to Earth-like planets are hit-and-run (they also included planetesimals), and Genda et al. (2017) reported 52% in their hybrid code (embryos only). We note that the definition of hit-and-run and its demarcation to other outcome regimes is not unambiguous, and here and in the following we use the term for all collisions that resulted again in two bodies according to our model (Sect. 2.2), independently of their mass ratio. This can result in (post-collision) mass ratios that would no longer be deemed hit-and-run by most definitions (Asphaug 2010; Genda et al. 2012), and can therefore somewhat overstate its occurrence rate.

Figure 11 shows the collision phase space of impact velocity versus impact angle. Impact velocity is normalized to the mutual escape velocity $v_{\text{esc}}=\sqrt{2GM/r}$, with the combined mass M and the sum of the bodies’ radiir. The impact angle α is defined as illustrated in Fig. 1, with α = 0° for head-on. As expected, impact velocities are always around or greater than v_esc and are rarely above 10 × v_esc. The highest values (and therefore the most destructive/erosive events) occur preferentially in the pp1 simulations, and are likely encounters between planetesimals, which are more likely on dynamically hot orbits. This is also the reason why pure embryo/protoplanet collisions (larger symbols) occur more frequently at lower velocities. Besides that, there seem to be no significant differences between the parameter setups (colors). Due to the limited number and fraction of planetesimals, dynamical friction is likely stronger in reality (see, e.g., O’Brien et al. 2006; Raymond et al. 2006; Chambers 2013), which would further damp eccentricities of embryos and protoplanets, presumably leading to somewhat gentler collisions (lower velocities).

Figure 12 illustrates the impact velocity as a function of time. The initially dynamically cold systems remain calm within the first ~1 Myr, and collisions preferentially result in accretion. This could be considered an artifact of the initial configuration, but also as the natural result of residual damping by the former gas disk and small body population. After the initial phase, orbits quickly become excited within the first few Myr, and appear to reach maximum excitation by ~10 Myr. The apparent decline at later times is probably rather due to the lower number of bodies, before indeed most highly excited planetesimals are eventually accreted or lost at the latest times. The last or at least one of the latest protoplanet impacts during Earth’s accretion is often linked to the formation of the Moon, where besides the classical (canonical) scenario (e.g., Canup et al. 2013), various different collision configurations are discussed, including formation ofmoons by hit-and-run (Reufer et al. 2012). Even though we cannot make any particular predictions, our data suggest that encounters between protoplanet-sized objects are very common in the later phases of the accretion process (Fig. 12), covering a wide range of impact velocities and angles in parameter space.

Again for the same sample of collisions, Fig. 13 illustrates (combined) water losses over time. This quantity is frequently found in the literature, but does not account for more complex processes like water transfer in hit-and-run, for example from a water-rich to a dry body. Not surprisingly, losses increase with time, with a naturally strong correlation to the increase of impact velocities in Fig. 12. There appears to be a noticeable clustering of water losses around 1, corresponding to encounters where the colliding bodies are either completely gravitationally dispersed, or at least highly erosive events which desiccate them entirely. The majority of collisions however, end up with combined water losses below ~20–30%, with a sizeable fraction even below 10%. Therefore, (average) losses are actually relatively low, but it is only the cumulative effect – perhaps including low-probability chance events – on the way to final planets that is of true significance. Following these evolutionary tracks for potentially habitable planets is the subject of Sect. 6.

Fig. 11 Impact velocity in units of vesc plotted vs. impact angle for all collisions in the regular scenarios with our collision model. Encounters exclusively among embryos and protoplanets (i.e., excluding planetesimals) are highlighted as larger symbols.

Fig. 12 Impact velocity in units of vesc over time for all regular scenarios with our collision model. Larger symbols mark collisions exclusively among embryos and protoplanets.

Fig. 13 Water loss in individual collisions over time for all regular scenarios with our model. Color coding is identical to Figs. 11 and 12. Water loss is defined as the fraction of the (combined) pre-collision water inventory that is not incorporated in either of the post-collision bodies. Larger symbols mark collisions exclusively among embryos and protoplanets.

6 Results II: water transport to potentially habitable planets

A central aspect of our work, and one that is crucial for assessing habitability, is to constrain initial water inventories of potentially habitable planets. We therefore analyzed (water) accretion of this subset in more detail. Due to the inherent stochasticity of long-term dynamical simulations and the uncertainties related to the extent of habitable zones themselves (see, e.g., Kasting et al. 1993; Kopparapu et al. 2013), we use a semi-major axis range of between 0.75 and 1.5 au to define “potentially habitable planets” for our purposes – broadly the space between Venus and Mars in the Solar System. This region is also highlighted in Figs. 2–4, showing the final systems and major properties of formed planets.

Fig. 14 All potentially habitable planets formed in the regular scenarios in the water-inventory vs. mass space. We note the log y-axis.

6.1 Main characteristics of potentially habitable planets

In all simulations, between zero and two planets form between 0.75 and 1.5 au, with one planet in most scenarios (Table 2). A complete list of these planets and their properties is provided in Table 3, and specific data on collisions related to our model are additionally summarized in Table 4. The overview in Fig. 14 reveals a distinctive twofold behavior of cJS and eJS simulations. For cJS, our model leads to significantly lower mass planets (around ~1 M_⊕) compared to PM, akin to the system-wide trend. Average water contents are ~30 O_⊕, compared toroughly doubled values for PM (Table 3). For eJS the situation is more complicated, and while masses are around 1 M_⊕ for both ourmodel and PM, final water inventories are highly scattered. The reason is again the high stochasticity of water transport with eJS, where often only a few water-rich bodies survive an initial removal phase of water-rich material due to forcing by the giant planets, as already detailed in Sect. 5.3. Water delivery to potentially habitable planets is even more stochastic than the global inventory, because it additionally depends on inward transport of water-rich material (instead of, e.g., remaining as/on stranded embryos in the outer disk). One manifestation of this is seen in the highly different (average) water contents of the PM-eJS-pp0 and PM-eJS-pp1 planets (Table 3), where closer inspection revealed that this is simply a result of the high stochasticity (as found also system wide, see Sect. 5.3). Therefore, it is justified and more reasonable to average over all PM-eJS planets, which leads to similar values with our model and PM of around 10 O_⊕ with eJS. The highly stochastic water transport with eJS thus disguises the actual differences between our model and simple PM in those scenarios. Stochastic variations aside, potentially habitable planets formed with cJS are likely more water-rich than with eJS (by a factor of two or more), while average masses are roughly equal and around 1 M_⊕.

6.2 Accretion and water delivery to potentially habitable planets

In this section, we discuss the accretion tracks of potentially habitable planets. Two particularly insightful examples are analyzed in Sect. 6.3. This is unfortunately only possible for a small, selected sample, but we provide comprehensive data on accretion histories of all planets to interested researchers⁶, in order to serve as a possible basis for in-depth modeling of further aspects and physical processes. Figure 15 illustrates growth of water content and mass for all potentially habitable planets formed in the regular scenarios. Finding the growth path of a planet is no longer entirely trivial once hit-and-run collisions come into play. Our results are based on the following algorithm: starting from the final planet, the evolution is traced backwards, where in collisions with a single survivor the more massive pre-collision body is followed. For hit-and-run, where the larger (smaller) pre-collision body practically always ends up as the larger (smaller) post-collision object, this characteristic is used to trace the respective pre-collision body further. A possible alternative would be to continue always with the more massive pre-collision body also for hit-and-run, which would not reliably trace the evolution passed these encounters however. Figure 15 makes it immediately apparent that our collision model (solid lines) frequently results in slower growth (see Sect. 6.5) and considerably less massive planets than simple PM (lower panel), and also a more complex evolution of their water inventories (upper panel), including frequent collisional erosion as well. In addition, the substantial scatter of final water contents and also the generally higher values for cJS (black + dark-green) than eJS (red + orange), are well visible. Final planets all end up between 0.5 and 1.5 M_⊕ with our model, while PM runs can result in considerably higher masses.

When one follows the evolutionary tracks for water content in Fig. 15 it may be surprising that often very few large steps – corresponding to single collisions – seem to decisively determine the final inventory. For the growth histories of all 18 potentially habitable planets formed with our model in the regular runs, this is indeed practically always the case, which means that the great majority of their final water content is delivered by a few – frequently only one or two – collisions. Furthermore, the water-delivering impactor is almost always an embryo or larger protoplanet, and is only very rarely a small (planetesimal) body (also in the pp0 scenarios). For the 18 planets, only twice did such a body deliver any appreciable amount of water (relative to the total accreted amount). This indicates that water-rich planetesimals are either removed from the system or incorporated into growing protoplanets before those can deliver water, thus showing a different behavior in terms of inward transport of these two types of source materials/bodies (Morbidelli et al. 2000). Therefore, the picture of inwardly scattered small bodies that directly deliver water to growing planets is rather not supported, at least not if the majority of mass has already accreted into embryos throughout the whole disk initially. We note that these considerations include all possibilities for the origin or formation of small bodies in our model – either pristine planetesimals from the outer disk, collisionally modified ones, or planetesimal-sized hit-and-run survivors.

This outlined mode of water delivery is basically the same in cJS and eJS scenarios. Exceptions are cases where the planet originates in a water-rich embryo from the outer disk beyond 2.5 au, which already start with 32–45 O_⊕, but have to be transported inwards over large distances. This happens for two of the eight planets with cJS (runs SPH03 and SPH08), and for none out of the ten planets with eJS. These planets are the most water-rich formed, with 64 and 50 O_⊕, respectively. The evolution of their water inventory can be easily followed in Fig. 15. All other 16 planets start formation inside 2.5 au, most from dry material within 2 au and some from embryos in the intermediate range between 2 and 2.5 au with WMF = 0.001 (but their initial water content is also modest ~0.5 O_⊕). The few water-rich impactors either originate directly from outer-disk material or have accreted considerable quantities from other bodies during their inward transport or diffusion. A significant difference between cJS and eJS is that the water-delivering embryos or protoplanets are often more water-rich with cJS, which explains the larger final water contents of those planets. A noticeable feature in several cases is that a first collision delivers a relatively small but significant amount (typically below 1 O_⊕), and a later one eventually delivers the majority of the final inventory. If the first collision does not happen the planet still ends up water-rich, but if the second impactor does not arrive, a relatively dry planet emerges. A general point to consider is that the few decisive events in the direct growth history do not take into account the impactors’ histories, which have their own, often complicated accretion paths, and are only rarely pristine bodies from the initial configuration. The decisive water-delivering collisions are sometimes accretionary and sometimes hit-and-run, with their typically low water-transfer efficiencies (see Sect. 6.4 and Fig. 16).

Fig. 15 Growth curves of all potentially habitable planets formed in the regular scenarios, with water inventories in the top panel (note the log y-axis) and mass in the bottom panel. The precise definition according to which evolutionary tracks are followed (especially for hit-and-run) is discussed in the text.

6.3 Two examples for accretion of potentially habitable planets

To illustrate the accretion of potentially habitable planets in detail, two typical examples are presented in Figs. 17 and 18, where growth of their mass as well as water content are plotted (thick solid lines). The planet in Fig. 17 forms in a cJS architecture (run SPH06) and ends up with ~24 O_⊕. It experiences four decisive collisions that deliver significant amounts of water, which is the highest number in all considered scenarios, but they are nevertheless representative and it is an interesting illustrative example. The first of these four collisions, a very grazing hit-and-run, happens after ~40 Myr and transfers less than 1 O_⊕ to the growing planet (note the log y-axis). After that, two further hit-and-run collisions deliver considerable amounts of water, where the impactor in the second one is actually the second hit-and-run survivor from the first one of these two events (~7 Myr before), which hits again and transfers some more of its water inventory. Finally, accretion of a water-rich protoplanet delivers ~24 O_⊕ at ~70 Myr. Snapshots of this latter graze-and-merge event are included in Fig. 16 (top row), illustrating water delivery by the impacting body (the target is color-coded as dry). Before the first of those events, practically only dry accretion occurs, and after the last one several collisions with dry bodies result in clearly visible losses of more than 5 O_⊕. The second example planet (Fig. 18) forms in an eJS architecture (run SPH15), and ends up with ~1.6 O_⊕, shaped by two decisive events. After 30 Myr of dry accretion, a hit-and-run encounter with a moderately water-rich protoplanet transfers ~0.35 O_⊕ to the growing planet, which is also illustrated as snapshots in Fig. 16 (bottom row), where the rather low water transfer efficiency can be clearly discerned. Only 6 Myr later, another (quite accretionary) hit-and-run with a more water-rich protoplanet delivers ~1.6 O_⊕.

For the planets in Figs. 17 and 18 we also include two additional growth curves, which are intended to provide some measure of the theoretical evolution without proper treatment of collisions. The first one gives the actual results with added cumulative collision losses, meaning that all material that is not included in either post-collision body is simply added (cumulatively). Final water and total masses are then the sum of the actual values (resulting from our model) plus all these losses during the accretion of the planet (also given in Table 4). In some sense this gives a lower limit for theoretical growth without proper collision treatment: while for events with a single post-collision body all material neglected as debris is now included (thus equivalent to PM), only a (typically small) fraction of the remaining material is included in hit-and-run.

The second theoretical growth curve, labeled “synthetic PM”, merges all material also in hit-and-run into the planet’s original growth path (again cumulative). In cases where the two survivors of a hit-and-run collision subsequently collide again, they are merged only once. Interestingly it turned out that defining such synthetic PM is not as straightforward as one might think, particularly due to hit-and-run in the original results, leading to convoluted accretion tracks. Eventually we decided to use the back-traced history of planets, that is, the regular results from our collision model since that istheir natural dynamical way to evolve, and to trace this path again forward in time while applying synthetic PM (instead of some direct and only forward-in-time approach). The obtained theoretical water contents are also included in Table 4. We note that even with synthetic PM the collision histories of the bodies impacting the growing planet are not also followed and replaced by PM, where complex accretion paths due to hit-and-run would complicate matters, for example avoiding that bodies or material are counted several times.

These two approaches yield growth curves that are significantly above the actual ones, clearly visible in Figs. 17 and 18. Due to the abundance of hit-and-run, it is furthermore not surprising that synthetic PM results in considerably larger values than counting only cumulative losses. A comparison of the average values for synthetic PM (Table 4) with the actual outcomes of our PM runs (Table 3) shows roughly similar results and indicates that synthetic PM provides quite reliable predictions for water contents if collisional effects were ignored. For the planet in Fig. 17 one obtains final water inventories of 35 and 75 O_⊕ for the two described approaches, respectively (vs. 24 O_⊕ as the actual result). For the planet in Fig. 18 these values are 3.5 and 5.4 O_⊕ (vs. 1.6 O_⊕ as the actual result).

Fig. 16 Snapshots of SPH collision simulations of two different scenarios (top and bottom row). The impactor (smaller body, in blue) approaches the target from the left. The color coding shows water on the impactor in blue and silicate material (basalt) in red. (Iron core material is hardly exposed and is light gray.) The target body was actually also simulated with a water shell, but is marked in red too in order to highlight water transfer from a water-rich impactor to a dry target. Top row: accretionary collision (more precisely “graze-and-merge”, where after a brief hit-and-run-like period the bodies fall back and merge), with vimp ∕vesc = 1.2, impact angleα = 28° and impactor-to-target mass ratio γ ~1:2. Bottom row: typical hit-and-run encounter with only relatively little water transfer, where vimp ∕vesc = 1.8, α = 47° and γ ~1:6. The collisionin the top panels is the last (at ~70 Myr) decisive water-delivering event for the planet in Fig. 17, the one in the bottom panels shows the first (at ~30 Myr) of two such decisive collisions for the planet in Fig. 18 (see Sect. 6.3 for details).

Fig. 17 Accretion history of the potentially habitable planet formed in run SPH06-cJS-pp0 (marked with superscript 1 in Tables 3 and 4, and also illustrated in Fig. 6). It forms at ~1.1 au, with 1.2 M⊕ and 24 O⊕ of water. Besides the actual results with our SPH collision model (solid lines), hypothetical growth curves with synthetic PM, and also with added cumulative losses are plotted – both defined and discussed in the text (Sect. 6.3). The crosses mark the respective property of the other (second) body in collisions, besides the traced planet.

6.4 Collision characteristics of potentially habitable planets

Table 3 provides overall collision numbers for all potentially habitable planets, while Table 4 summarizes more detailed collision information for the subset formed with our model. In the regular scenarios and with our model the average number of collisions in the growth history of a potentially habitable planet is ~20 with pp1 and ~30 with pp0. This is not surprising since with pp1 planetesimals can accrete mutually before hitting the growing planet, but they cannot with pp0. However, the total number of collisions is only of secondary importance anyway, since it depends strongly on the initial number of bodies, which is determined by computational limitations, not physics. More interesting is the number of impacts by large bodies (embryos and protoplanets), because the initial embryo population is accurately resolved according to our physical model. This number averages to approximately ten collisions, with pp1 and the same with pp0. In comparison, potentially habitable planets formed with PM experience roughly one-third less impacts by embryos or protoplanets, presumably due to the absence of hit-and-run. With our model, the prevalence of hit-and-run in the accretion histories of potentially habitable planets is similar to what was found and discussed in a system-wide context in Sect. 5.5. In the regular scenarios, about 50% of all collisions are hit-and-run, also if only impacts by embryos and protoplanets are considered, and even ~75% in embryo-only simulations.

Hit-and-run collisions where a water-rich target is hit by a considerably drier projectile can strip water from the target body, even thoughthis becomes only really transformative if hit by an approximately equally massive or larger impactor (e.g., Asphaug 2010; Burger et al. 2018). In the accretion histories of all 18 potentially habitable planets formed with our model in the regular scenarios, only once is the growing planet hit by a more massive impactor, and collisions with equally massive bodies are also relatively rare. Most impacts are by significantly lower mass objects (see Figs. 17 and 18). However, energetic individual events as well as a large number of smaller impactors can lead to considerable collisional water losses. Especially interesting for the delivery of water on the other hand are encounters where a drier target body is hit by a water-rich projectile, with potential transfer to the target. However, for hit-and-run this transfer is often not very efficient, especially compared to accretionary (e.g., more head-on) collisions (Burger et al. 2018, 2019). For commonly occurring collision parameters, typically less than ~25% of the projectile water is transferred in hit-and-run, but 50% or more is transferred in accretionary impacts (cf. Fig. 16). The rationale of past studies behind still assuming PM (in spite of usually recognizing the high frequency of hit-and-run) was often that the two hit-and-run survivors remain on similar orbits and re-collide soon after the first encounter. Our results support this assumption only to a low degree, as was also found in a recent, dedicated study by Emsenhuber & Asphaug (2019). Table 4 lists the total number of collisions, along with the number of hit-and-run and how many of those led to re-collision at a later time in brackets (and the same for purely embryo/protoplanet encounters). Re-collisions happen on average in only ~25% of cases, thus the majority of bodies that go into hit-and-run do not directly meet again. For comparison, Chambers (2013) found that 54% of hit-and-run events (considering only those among embryos and protoplanets) lead to re-collision of the same bodies later. A follow-up question that presents itself pertains to the nature of the re-collisions. Even though this analysis suffers from low sample sizes, our results indicate (not included in Table 4) that about half of all re-collisions after hit-and-run are again hit-and-run (Emsenhuber & Asphaug 2019, found a somewhat reduced probability for repeated hit-and-run). Therefore our data suggest no special trend in re-collision outcomes, and do not support the hypothesis that the majority of hit-and-run encounters are followed by subsequent accretion. The time between hit-and-run and re-collision varies over a very wide range (Table 4). For those planets where relatively good statistics (high hit-and-run and re-collision numbers) exist, median values are on the order of 1 to 10 Myr, and re-collision times can be as high as several 10 Myr (Chambers 2013, found a mean time of 1.3 Myr). Very immediate re-collisions seem to be rare – in our data (on potentially habitable planets) only one such event occurs after just 70 yr, where the next shortest time to re-collision is ~6000 yr. The high frequency of hit-and-run and the often low water transfer efficiency, combined with low re-collision rates, makes it particularly important to accurately model this type of encounter, especially if water delivery is dominated by few decisive collisions with large bodies (embryos or protoplanets), which naturally tend to be more similar sized (and thus more likely hit-and-run) than delivery by smaller bodies.

Figure 19 shows water losses over time for collisions in the growth histories of potentially habitable planets. This is a subset of Fig. 13, where all collisions in the regular scenarios are plotted. One apparent difference is the practical absence of very high water losses in Fig. 19, which are presumably often either very destructive events (thus not directly part of planetary growth), or happen in the innermost part of the disk (see Sect. 7.1). Otherwise there are no clear differences between either the giant planet architecture (cJS vs. eJS) or the planetesimal–planetesimal collision model (pp1 vs. pp0) in Fig. 19. Water losses are mostly between 1 and 10%, while very erosive events, and even losses significantly above 10% are rare. Knowing water losses in all collisions on the way to the final planet offers the possibility to compare their cumulative sum to the planet’s final water inventory. This is listed in Table 4 (third column) and plotted in Fig. 20. Cumulative water losses are typically somewhat less, but make up a sizeable fraction of the final inventory, despite considerable scatter especially with eJS. Considering also the histories of the impactors and associated losses (which are not included in Fig. 20 and Table 4), a reasonable first-order estimate might be that total losses in thewhole accretion tree of a planet are roughly equal to its final water content (the diagonal line in Fig. 20). For comparison with simple PM, the addition of cumulative losses likely still gives lower results, since in actual PM runs the whole impactor (including a potential second hit-and-run survivor) is always merged into the growing planet. This is also confirmed by the theoretical growth curves in Figs. 17 and 18, where synthetic PM (intended to emulate actual PM outcomes) gives considerably larger water content (and mass) than merely adding cumulative losses (see Sect. 6.3). We note that cumulative mass losses are only 10–20% in the regular scenarios (Table 4), while cumulative water losses are usually much higher, a clear indication that (surface) water is stripped more easily in collisions than other material.

In order to asses the actual efficiency of water delivery in individual collisions we use a modified version of the water accretion efficiency (Burger & Schäfer 2017), $${\Xi }_{\text{w}}=\{\begin{array}{cc}\frac{\alpha }{M_{\text{w},\text{p}}}& \hspace{1em}\hspace{1em}\text{if}\hspace{1em}\alpha >0\\ \to 0& \hspace{1em}\hspace{1em}\text{if}\hspace{1em}\alpha \to 0\\ \frac{\alpha }{M_{\text{w},\text{t}}}& \hspace{1em}\hspace{1em}\text{if}\hspace{1em}\alpha <0\end{array}$$(4)

where α = M_w,tf − M_w,t. M_w,t and M_w,p denote the (absolute) pre-collision water masses on target and projectile, where “target” always means the body that is followed (even if it is the smaller one). M_w,tf denotes the mass of water on the post-collision body that originated from the target (relevant only for hit-and-run, because otherwise there is only one surviving body). The interpretation of Ξ_w is the following: if the target body has increased its (absolute) water inventory during a collision (α > 0) then Ξ_w ∈ (0, 1] and the amount of increase in units of the (pre-collision) projectile inventory is given. Conversely, if the water inventory of the target decreases (α < 0) then Ξ_w ∈ [−1, 0) and the amount of decrease in units of the (pre-collision) target inventory is given. For practically unchanged water contents Ξ_w is close to zero. This definition provides an easy-to-interpret measure of the efficiency of water accretion and simultaneously also erosion, in all cases, independent of whether both bodies are water-rich or if either of them is dry. A plot of Ξ_w over time in Fig. 21 shows that in the majority of collisions it remains close to or below zero, which indicates losses of water inventory, albeit only by a few percent in most cases. Only a handful of collisions erode more than 10% of the target’s water content (Ξ_w below −0.1), and not a single one erodes more than 40%. This is likely a simple consequence of the fact that impactors are mostly smaller than the growing planets, especially in their later evolution (cf. Figs. 17and 18), and it is difficult for impactors to erode larger targets, even more so for hit-and-run (Asphaug 2010; Burger et al. 2018). A comparatively small fraction of collisions results in clearly positive Ξ_w, but these are distributed rather uniformly up to Ξ_w = 1. Results by Burger et al. (2019) indicate that the larger positive values are likely from accretionary collisions, while for hit-and-run they found values of ~0.3 at most, and frequently considerably less. Positive values of Ξ_w occur not only in the few decisive water-delivering events (see Sects. 6.2 and 6.3), but also whenever more water is delivered than removed. However, this becomes increasingly less likely for drier impactors, and it is the prevalence of significantly drier impactors in the growth histories of potentially habitable planets (clearly visible in Figs. 17 and 18) that cause Ξ_w to be at least slightly negative in the majority of collisions. We note that with simple PM, Ξ_w is always 1, and none of this diverse behavior is modeled at all.

Fig. 18 Accretion history of the potentially habitable planet formed in run SPH15-eJS-pp1 (marked with superscript 2 in Tables 3 and 4, and also illustrated in Fig. 6). It forms at ~0.9 au, with 1.03 M⊕ and 1.59 O⊕ of water. Lines and symbols have the same meanings as in Fig. 17.

Fig. 19 Water losses in individual collisions over time during accretion of all 18 potentially habitable planets in the regular scenarios. Water loss is defined as the fraction of the (combined) pre-collision water inventory that is not incorporated into either of the post-collision bodies.

Fig. 20 Cumulative water losses (from their individual collision histories) vs. final water inventory for all 18 potentially habitable planets formed with our model in the regular scenarios. The diagonal line indicates cumulative losses equal to the final water content.

Fig. 21 Water accretion efficiency Ξw in individual collisions over time during accretion of all 18 potentially habitable planets in the regular scenarios. Positive values indicate an increase in water inventory (in units of the projectile’s water content), while negative ones indicate a decrease (in units of the target’s water content). See the text for a precise definition of Ξw and further discussion.

6.5 Timing of accretion and water delivery to potentially habitable planets

Mass accretion timescales of potentially habitable planets are indicated by t₅₀, t₉₀ and t_last-emb-col in Table 3, and for their water inventories by t_25,w, t_50,w, and t_90,w. These represent the accretion of 25, 50, and 90% of their final mass and water content and the time of the last pure embryo/protoplanet collision (i.e., excluding planetesimals), and are all computed from the final simulation time backwards (not from t = 0 forwards). There appears to be no obvious strong difference between cJS and eJS in the regular scenarios, and they are also considerably lengthened by roughly a factor of two compared to PM, in concordance with accretion times of the entire systems (Sect. 5.4). However, the formation of individual (potentially habitable) planets generally proceeds faster, which is not surprising since only a single object is considered, and additionally accretion in the inner part of the disk is generally faster than for planets farther out. With our collision model it takes potentially habitable planets on the order of 10–50 Myr to accumulate 50% of their final mass, and 40–140 Myr to accrete 90%. In some cases t₉₀ and t_last-emb-col coincide, but frequently there are still further collisions with other embryos or protoplanets ahead even after a planet has accreted almost all of its mass, presumably often (not very accretionary) hit-and-run encounters. For comparison, PM simulations by Chambers (2001) found that accretion of 50% (90%) of Earth analogs is completed in about 20 (50) Myr, and O’Brien et al. (2006) found similar values, even though they considered all formed planets in the system. In very high-resolution simulations, Fischer & Ciesla (2014) found mean values of ~60–80 Myr for the last collision with an embryo or larger body for their habitable-zone planets. These results are consistent with our PM runs (albeit oftensomewhat larger), but are considerably lower than those of our model. There is a potential conflict especially with Hf-W measurements which seem to indicate shorter timescales for core formation on Earth. However, this is presumably connected rather to the last accretionary giant collision, and the influence of large-scale hit-and-run interaction (as happens regularly with our collision model) is unclear. O’Brien et al. (2006) also found somewhat shorter accretion times if strong (and thus probably more realistic) dynamical friction is included, which would indicate that our results have not yet reached resolution convergence.

The water accretion times t_25,w, t_50,w and t_90,w and particularly their frequent equality (Table 3), mainly show what has been described in Sects. 6.2 and 6.3 already, namely that the vast majority of the final planets’ water inventory is delivered in a few – often only one or two – decisive collisions, while all other impactors deliver only small amounts, or remove water. Radial mixing of water-rich material to the region of potentially habitable planets was found to take on the order of 10–30 Myr (see second panel in Fig. 8). This is confirmed by their accretion histories, where the decisive water-delivering collisions rarely happen prior to 10 Myr, and often also much later, especially with eJS (see also Fig. 5). Exceptions are the planets whose accretion starts from a water-rich embryo (those formed in runs SPH03and SPH08), which are water-rich from the beginning but nevertheless have to move inwards from their initial orbits. Once accreted, the chance of long-term water retention is certainly a strong function of the planet’s current mass, but also depends on its final mass, which determines how much collisional accretion (or erosion) is still ahead. Our results indicate that the bulk of final water inventories is often delivered late, between t₅₀ and t₉₀ or after t₉₀ (see Table 3). Deliveryafter t₉₀ seems more frequent in eJS scenarios than with cJS, and water retention is therefore more likely.

Anotheraspect of the timing of collisions is the intervals between those events for a growing planet, especially impacts by large bodies (embryos or protoplanets). The elapsed time since the last collision is plotted in Fig. 22 for the (direct) growth paths of potentially habitable planets. Giant collisions with embryos or protoplanets (large crosses) generally deposit enough energy to not only vaporize surface water and ice, but can also produce local or even global magma oceans. Therefore, they determine the physical state of water for typically a couple of Myr (until re-solidification of the surface at least), and thus not only the evolution until then but also how subsequent impactors encounter the growing planet, depending on the elapsed time since the last large event (see Sect. 7.6 for further discussion). The average time since the last collision naturally increases as planets grow and body numbers decrease. The diagonal lines in Fig. 22 indicate that intervals between collisions are almost always larger than 1/100 of the current system age (above lower line) and for the majority at least 1/10 of that age (above middle line). Typical intervals between collisions are therefore above 1 Myrat t = 10 Myr, and already above 10 Myr at t = 100 Myr. The relative scarcity of planetesimal impacts (small crosses not over-plotted by large ones in Fig. 22) is essentially a resolution effect, where a better-resolved planetesimal population would certainly lead to an increased number of planetesimal impacts between collisions involving only protoplanets and embryos. A closer look reveals that in the very end phase of accretion, considerably more planetesimal impacts occur with pp0 (dark-green and orange) than with pp1 (black and red), which indicates that a sizeable small body population survives longer with pp0 (see also Fig. 6).

Fig. 22 Time since the last collision (evaluated each time a collision happens) from the perspective of a growing planet as a function of time, for all 18 potentially habitable planets in the regular scenarios with our model. Small crosses count all collisions, large ones only impacts by embryos and protoplanets (i.e., excluding planetesimals). Each collision is assigned two data points in general – a small and a large cross – and if the last impacting body was an embryo or protoplanet then they coincide. The continuous and dashed lines indicate elapsed time (since the last collision): since simulation start (upper line), 1/10 of that time (middle line), and 1/100 of that time (lower line).

7 Discussion and conclusions

The direct combination of long-term N-body integrations and dedicated hydrodynamical simulations of all occurring collisions is basically a simple and obvious approach, but has nevertheless been applied only very few times to self-consistently model the collisional evolution in N-body planet formation simulations (Genda et al. 2011, 2017). Moreover, to our knowledge, this approach has never been used to also trace volatile constituents like water. The obviously high computational demands may have been one reason, the intricacies involved in the realization of such a method might have been another. Although it appears rather straightforward, we found it to be a formidable task to self-consistently implement the interface between the two worlds of N-body computations and SPH collision simulations, which resulted in more than 7000 lines of additional code, so far.

7.1 System-wide collisional evolution and water transport

The results of our perfect merging (PM) comparison simulations are generally in agreement with existing studies. In particular the influence of giant planets on eccentric orbits (eJS) leads to vastly reduced water contents compared to circular gas giants (cJS), where (system-wide) not more than a few tens of O_⊕ survive with eJS, but ~100–150 O_⊕ survive with cJS (Table 2). With eJS, strong dynamical forcing by the giant planets rapidly (~5–10 Myr) removes the majority of water-rich material from the system, leaving only a few water-rich bodies (cf. Fig. 6), which additionally renders final water content and distribution highly stochastic (with large scatter between scenarios). With cJS, an initial removal of water-rich material occurs too, but a substantial amount remains, which is the main reason for the much larger final inventories in PM runs with cJS and makes them also more robust and less stochastic. With our collision model on the other hand, final (system-wide) water contents with cJS are not much higher than with eJS (for both ~30–40 O_⊕ on average), where mainly two mechanisms appear to further reduce water amounts with cJS (see also Sect. 5.3): first, several tens of O_⊕ are lost to collisional erosion, but only a fraction of that with eJS, presumably simply because eJS systems contain much less water after the initial removal phase, while there are no large differences in collision characteristics (especially impact velocities) between cJS and eJS (Sect. 5.5). Second, our results indicate that the accretion phase is truly lengthened, by roughly a factor of two with realistic collision outcomes and especially accurate treatment of hit-and-run encounters (see Sect. 5.4, and also Sect. 6.5). The associated slower growth allows for further removal of water from the system, considerably more so with cJS, presumably again due to higher abundances after the initial removal phase.

Unlike with cJS, the resolution (number of bodies) in the regular scenarios appears to be still too low to obtain robust results on water transport with eJS, at least for initial configurations like ours, where embryos have already also formed in the outer disk, and with a rather sharp transition between the water-poor inner and water-rich outer regions. In contrast to that would be more smooth initial water distributions (Quintana & Lissauer 2014), possibly as a result of earlier radial mixing (e.g., Carter et al. 2015), which would likely also make results with eJS less stochastic. Eventually a higher number and/or fraction of small bodies (planetesimals) might be required to model dynamical friction more realistically and better resolve the water-rich outer disk. However, eJS simulations by O’Brien et al. (2006) with such higher resolution found that even less water-rich material remains to later diffuse inwards (or remains as “stranded embryos”) compared to our lower-resolution simulations. This would suggest that a better-resolved planetesimal population can robustly model their radial mixing and water delivery behavior, but inwards water transport by embryos (if they form in the outer disk) ultimately remains a stochastic process in eJS-like environments. We note that inward diffusion of water-rich material itself is likely a robust process dominated by interaction with the debris disk and not the giant planets, where water is reliably transported inwards provided it is not removed before. The same behavior has also been observed without any giant planets at all (Quintana & Lissauer 2014). In embryo-only scenarios the quick removal of water-rich material⁷ can have an even profounder influence, again more so with eJS, where the (chaotic) fate of individual water-rich embryos is often decisive. At least in one scenario (EO-PM17-eJS, see Fig. 3) all water-rich material is removed before delivering any water and both formed planets are composed exclusively of dry material from inside 2 au.

A central part of our model is the accurate treatment of material and water transfer and loss in hit-and-run encounters, which are prevalent in collisions among embryos and protoplanets with a roughly 50% occurrence rate (see Sect. 5.5 and Fig. 16). The often-stated claim that the two hit-and-run survivors mostly remain on similar orbits and re-collide (and accrete) soon after the first encounter is not supported by our results (Sect. 6.4), and therefore accurate treatment of hit-and-run appears to play a decisive role in water transport and late-stage accretion in general. Hit-and-run can be transformative especially for the smaller body of the colliding pair and strip it of mantle material and water (Burger et al. 2018), but still let the colliding bodies further evolve as individual objects. Therefore, it has been suggested that unaccreted second-largest bodies might often be shaped by hit-and-run events with the largest body in the region (Asphaug 2010). This is not supported by the accretion histories of the 18 potentially habitable planets formed with our model in the regular runs, where merely once the growing planet is hit by a larger body, and only rarely by an approximately equally massive one (cf. Figs. 17 and 18). However, these planets are often the largest formed (Fig. 2), so we also checked all other planets in those systems, with interesting results. In most scenarios there is no collisional interaction at all between the planets in their later growth phases, that is, they form independently, and also mostly by impactors that are smaller or equally massive at most. There are some instances where growing planets have hit-and-run encounters with larger bodies which are then removed from the system or accreted elsewhere, but very transformative hit-and-run collisions with larger impactors appear to be rare exceptions rather than the rule. In 5 out of 20 scenarios growing protoplanets do collide late in their evolution but still end up as separate planets, but in all except for a single case they have very grazing (and thus not particularly transformative) hit-and-run encounters. This happens for rather unsurprising configurations, like the two close inner planets in the SPH01 run (see Fig. 2), but also for planets that eventually end up on widely separated orbits, like in SPH02 or the middle and outer planets in SPH10. In only one instance, the SPH07 scenario, is a smaller planet decisively shaped by survived (hit-and-run) collisions with a larger one, where the outer planet has two late hit-and-run encounters with the roughly three times more massive inner planet. These collisions lead to considerable water losses on the more water-rich inner planet, but especially strip the outer planet of almost half its mass – preferentially rocky mantle material – and its CMF increases from 0.27 to 0.43. This behavior might be a simple consequence of the rather narrow allowed parameter range for such encounters: too grazing and they undergo no significant interaction, too central and they are more likely to merge and only a single planet forms instead of two (perhaps along with a smaller hit-and-run survivor). It is only in between, and with sufficient v_imp ∕v_esc, where highly transformative hit-and-run collisions are possible, which still leave two planet-sized objects. In addition, the protoplanets have to eventually achieve separated, long-time stable orbits again after those events. While the mostly larger mass (gravity) of growing planets aids in water accretion from the very few water-rich impactors (in accretionary as well as hit-and-run collisions), this factor also hinders erosion of existing water by smaller impactors in hit-and-run (Burger et al. 2018). In this sense, collisional water delivery is commonly more efficient than collisional erosion, especially during the later evolution when impactors are typically smaller than their targets.

A related observation is the high abundance of close-in (around 0.5 au) low-mass planets with high CMF, which form with our collision model but not with PM. Examples can be easily identified, for example in all scenarios in the first row of Fig. 2 (runs SPH01, SPH06, SPH11 and SPH16), and are formed in roughly half of all systems (see also first and fourth panels in Fig. 4). As evident from the discussion above, not individual mantle-stripping events (e.g., by neighboring planets), but rather a long sequence of energetic collisions gradually erode mantle material and increase their CMF, regularly up to 0.5 or more. The dynamical environment close to the star leads to high relative velocities, and planets that experience long chains of often erosive collisions with several times v_esc. The impactors frequently also have high CMFs, which suggests that they have undergone erosive encounters of their own, and that altogether mutual erosion towards ever more iron-rich bodies is typical for that region of the disk. This appears to be not only a possible natural explanation for Mercury’s large CMF – without the necessity of some individual chance event – but also suggests that similar planets might be common in extrasolar systems. We want to emphasize that our model was designed primarily to track collisional water transport and its limitations especially with respect to small-scale refractory debris (see Sect. 7.4) make it hard to assess the accuracy of these results. Nevertheless, the possibility that the outlined collisional evolution might be a regular and robust pathway to small, close-in, iron-rich planets certainly merits further investigation.

7.2 Water delivery to potentially habitable planets

A main focus of this study is collisional growth and water delivery to potentially habitable planets, which we define as lying between 0.75 and 1.5 au. In the regular scenarios and with cJS, potentially habitable planets with ~1 M_⊕ and average water contents of ~30 O_⊕ form with our model, while PM runs result in masses mostly above ~1.5 M_⊕ and roughlydoubled water inventories. Therefore, even with realistic collision treatment these planets end up still more water-rich than Earth, but are considerably closer than the very high final water contents in PM simulations. In the respective eJS scenarios, potentially habitable planets appear more Earth-like on average, with ~1 M_⊕ and water contents of ~10 O_⊕, but substantial scatter, typical for our eJS scenarios as discussed above (Sect. 7.1). Nevertheless, individual planets formed in the eJS runs may represent viable pathways to Earth-like planets with suitable water contents.

Our results indicate that the great majority of delivered water arrives by a few – often only one or two – decisive collisions, and that those are mostly encounters with other large bodies (embryos and protoplanets). The picture of (more or less direct) delivery by water-rich planetesimals scattered inwards from the outer disk (beyond the ice line) appears to be unsupported. Small bodies are rather either removed from the system or accreted by growing, inwards-diffusing protoplanets, which then in turn deliver water to growing planets. There are particularly two caveats, both with respect to our initial conditions: first, there is initially already more (water) mass in embryos than in planetesimals (~75 vs. 25%) – mainly due to computational constraints – and second, we assume that embryos have already formed uniformly throughout the whole disk at the start of simulations. With similar initial conditions like ours and an architecture similar to cJS (but PM), Raymond et al. (2007) found also that the majority (> 90%) of water is delivered by larger bodies to planets in the habitable zone.

Based on our results we propose the hypothesis that late-stage water delivery to forming planets, particularly in the habitable zone, is a stochastic process, dominated by a few decisive events instead of a more continuous accretion of water-rich material, at least if embryos have already formed throughout the disk. The water-delivering impactors are most often at least embryo-sized or larger, and are only rarely small bodies (planetesimals). Such similar-sized collisions are naturally sometimes in the accretionary and sometimes in the hit-and-run regime, where the typically low efficiency of water transfer in hit-and-run strongly limits the amount of delivered water (Burger et al. 2018, 2019), and makes accurate modeling thereof particularly important. The results of Burger et al. (2019) show that in relatively gentle head-on collisions water accretion efficiencies (Ξ_w, see Sect. 6.4) can be up to 0.8 (for v∕v_esc = 1.5), while for hit-and-run they are less than half of that at most, and often considerably less. We note that for simple PM, Ξ_w is always 1, and all of the impactor’s water inventory is delivered. After hit-and-run, the second survivor can in principle strike again and deliver more water, as seen for example during accretion of the planet in Fig. 17 (see Sect. 6.3), but our results indicate that this does not happen in the majority of cases (Sect. 6.4). If true, the above outlined hypothesis implies high scatter in water contents when formed under such conditions (cf. Table 3), where a few chance events can make the difference between practically dry, moderately water-rich (like Earth), or “ocean planets” (see also Morbidelli et al. 2000). The hypothesis can be tested, particularly with a better-resolved planetesimal population (increased numbers and larger initial fraction of mass and water in planetesimals). In PM simulations with such initial conditions (and ~cJS), O’Brien et al. (2006) found that still only a minority of accreted outer disk material comes from planetesimals, while Fischer & Ciesla (2014) find (with cJS and eJS) that more than half of the water-rich material accreted by habitable-zone planets originates in planetesimals (but they do not specify how much is accreted by other embryos or protoplanets before final delivery). Under the assumption that embryos only form in the inner disk, while the outer disk is initially populated exclusively by planetesimals, Raymond et al. (2007) found (with ~cJS) that a considerable amount of water is delivered by planetesimals, a less stochastic mechanism. Izidoro et al. (2013) on the other hand report that only a small fraction of planetesimals from such an outer disk are incorporated into final planets (the rest are presumably ejected or otherwise lost), but they do not specify how many of those are first accreted by other embryos or protoplanets. It would be interesting to test and compare all these results with realistic collision outcomes instead of simple PM.

Two of the 18 potentially habitable planets formed in the regular scenarios with our model start accretion as water-rich embryos from beyond 2.5 au with several tens of O_⊕, then migrate inwards during accretion, and eventually end up at 1.1 and 1.2 au (in runs SPH03 and SPH08). Both planets form with cJS, presumably because the chances of such an evolution are simply much lower if the majority of water-rich material is rapidly removed (with eJS) before it can set in. During the growth of both planets theyaccrete additional water-rich material, but also experience considerable collisional losses at later times (see Fig. 15). However, collisional erosion is difficult once protoplanets have grown to a certain mass and most impactors are significantly smaller (Asphaug 2010; Burger et al. 2018). Indeed these two planets end up as the most water-rich of the 18 potentially habitable planets – essentially “ocean planets” with 64 and 50 O_⊕. These two examples indicate that if sufficiently large bodies can form in the water-rich outer disk, and the environment is dynamically calm enough, then there is a chance that they will diffuse inwards and retain most of their water inventory despite collisional erosion, or even accrete more on the way – a viable pathway to very water-rich Earth-sized planets in the inner system. There is also an interesting connection to PM, where water-rich bodies are always entirely accreted, including their whole water inventory. In some sense this is also the case for these planets, where the majority of the water of their seed embryos is retained. Their final water contents are indeed relatively close to results of PM runs (Table 3). Planets that accrete dry until later water delivery on the other hand, receive only a (often small) fraction of the impactors’ water inventories with realistic collision treatment, much unlike PM.

Especially with respect to previous work on water delivery during late-stage accretion, questions pertaining to differences between ourmodel and PM and the amount of collisional water loss present themselves. For eJS scenarios, a direct comparison is not very meaningful due to the very stochastic evolution (but see Tables 2 and 3). For cJS, system-wide water contents end up (on average) between three and four times lower than with PM, and for individual potentially habitable planets the difference is still a factor of about two. Including realistic collision outcomes changes the process of late-stage accretion itself to some degree, which complicates comparison with PM, where frequent hit-and-run encounters lengthen growth times and can lead to complex, sometimes repeated, interaction among protoplanets. Moreover, high-energy and often erosive collisions can lead to results strongly deviating from PM, especially in the inner disk due to high relative velocities. To assess the deviations caused by not accounting for realistic collision outcomes we also computed “synthetic PM” accretion tracks from the actual results of our model (Sect. 6.3, and especially Figs. 17 and 18), which appear to be quite reliable indicators. In addition, we found that cumulative collisional water losses are mostly somewhat below final inventories (Sects. 6.3 and 6.4, and especially Fig. 20), but considering also the accretion histories of all impactors probably leads to typical losses on the order of the final water content. Collecting all the evidence together leads to the conclusion that collisional water losses likely reduce final inventories of potentially habitable planets by a factor of two or more in most situations.

7.3 Initial conditions

Our initial conditions are intended to represent the onset of late-stage accretion in the terrestrial planet region, once previous growth has led to a system of planetary embryos embedded in a large number of remaining smaller bodies. We assume that the gas component has already dispersed before larger protoplanets are able to form, which is still an open question also for the early Solar System. Our giant planet architectures (cJS/eJS) are in contrast to migrating gas giants, which can influence debris disks to an even greater degree (as in Grand Tack scenarios; Walsh et al. 2011; O’Brien et al. 2014), and giant planets that run into dynamical instabilities. Despite the great variety of known extrasolar systems, similar evolutionary paths might still be frequent. The initial amount and distribution of water, which strongly depends on the composition and conditions in the protoplanetary disk, is generally not well constrained, where meteorite data provide some information at least for the early Solar System. Our initial conditions generally follow what is often referred to as dry accretion, that is, no significant amounts of water are incorporated locally into building blocks in the inner disk, and we also assume that no earlier radial mixing has smoothed out the water distribution (see, e.g., Carter et al. 2015). Building blocks beyond the ice line on the other hand could be even more water-rich than our 5% (see O’Brien et al. 2018, and references therein), where for example 10% is not unreasonably high (for carbonaceous condrites, Morbidelli et al. 2000) and could have a large impact on final water inventories. In addition, the available budget of radionuclides (e.g., ²⁶Al) can alter initial water contents over a wide range (Lichtenberg et al. 2019). Even though we chose to limit this study to one particular initial water distribution, our collision model provides the straightforward possibility to apply arbitrary distributions as a simple post-processing step to whole simulations, for example models for local initial water absorption (Izidoro et al. 2013) or smoothly varying WMF (Quintana & Lissauer 2014). The mechanism for that is already applied in the frequently necessary pre-collision upscaling and post-collision correction of water contents in order to reasonably resolve them in SPH simulations (see Sect. 3.2).

Due to naturally limited computational resources our simulations also suffer the common problem of mostly orders-of-magnitude too large and massive planetesimals. With our model this additionally affects collisions, which typically involve bodies that are considerably too massive (in planetesimal–planetesimal collisions and planetesimal impacts onto larger objects), while their frequency is simultaneously too low. In some sense this can be considered as sampling only the top most part of the unknown planetesimal mass-distribution and neglecting (collisional) interaction with a large number of smaller bodies, except for the most top-heavy distributions. It is difficult to assess the consequences of an impactor flux composed additionally of a large number of small-to-medium-sized bodies instead of exclusively high-mass planetesimals. Particularly for water (and other volatiles), its location and physical state are likely important factors, especially when the bulk of a body’s water inventory is vaporized, like after a large collision (further discussed in Sect. 7.6). For impact erosion of (proto)planetary atmospheres, smaller bodies have been found to be much more efficient (per unit mass) than larger ones (see, e.g., Schlichting et al. 2015; Schlichting & Mukhopadhyay 2018; Lammer et al. 2019). This is mainly because most of their impact energy ejects atmosphere locally, while for larger impactors strong shocks have to traverse the target to eject atmosphere globally. In this study we performed a large number of more generic (e.g., no atmospheric component) collision simulations, where water is modeled as surface fluid, and evaluated water losses in all these events. This gives a quite robust sample, covering mass ratios down to ~ 10⁻³, even thoughthe accuracy towards the lowest mass ratios is likely limited by SPH resolution. The results are plotted in Fig. 23 and appear to support the notion that smaller impactors (per unit mass) lead to higher water losses, even though this is difficult to quantify precisely. We note that part of this effect could be due to hit-and-run, which becomes more likely in more similar-sized collisions but usually results in lower losses than accretionary or erosive events. Our results for water delivery to potentially habitable planets suggest that only a small number of decisive collisions with other large bodiesare responsible for delivery of the great majority of their final inventories, at least in their immediate collision history (Sect. 7.2). If this hypothesis proves to be correct, then at least for this aspect of water delivery, small bodies and their actual distribution would be of only little consequence. Ultimately a better resolution and modeling (e.g., via tracer particles; Levison et al. 2012) of planetesimals is required, which ideally also samples their mass distribution to some degree (instead of uniform masses). More realistic initial conditions will be a central aspect of future applications of our hybrid framework (see Sect. 7.7).

Besides the regular scenarios, we performed embryo-only simulations where all mass is assumed to have been initially accreted into approximately 30 planetary embryos in order to compare different resolutions and to study the influence of the planetesimal population. Similar setups have been used in the few existing direct hybrid simulations which are comparable to our approach (Genda et al. 2011, 2017). A central aspect of the dynamical evolution in these scenarios is missing dynamical friction by smaller bodies, which leads to generally more excited populations of embryos and protoplanets, and presumably more energetic collisions, eventually leading to planets on often more eccentric orbits. Additionally, water transport depends exclusively on the evolution of only a few water-rich embryos, which leads to more stochastic water delivery in general, and particularly with eJS (see Sect. 5.3). Altogether the embryo-only simulations show the same general trends for masses, water contents, collisional losses, and accretion timescales (cf. Tables 2–4), even though absolute values deviate sometimes considerably and outcomes have higher scatter. Therefore, results from such scenarios are mostly of limited accuracy, but could nevertheless be useful, for example in parameter studies in a less stochastic environment(like cJS), and with a focus on basic properties and fast computation with run-times of a few weeks instead of several months.

Fig. 23 Water losses in individual collisions as a function of the projectile-to-target mass ratio (where the target is always the larger body), for all regular scenarios with our model. Water loss is defined as the fraction of the (combined) pre-collisionwater inventory that is not incorporated into either of the post-collision bodies. Larger symbols indicate collisions among embryos and protoplanets alone (i.e., excluding planetesimals). Slightly negative values are merely numerical artifacts and can be assumed to be zero.

7.4 The collision model

Our physical model aims for realistic treatment of collisional water transport, which means simulation of the dynamical (long-term) evolution with high accuracy and collision of bodies only in actual physical encounters (including tidal collisions to some degree). Our approach of dedicated SPH simulations has the great advantage of accurate and reliable collision outcomes, and in principle allows the simulation of all desired physics for each event. Particularly for the application in this study, it is also well-suited (and currently without alternative) for modeling the rather subtle effects of collisional water losses and also transfer in hit-and-run encounters. The frequently used semi-analytical collision outcome model by Leinhardt & Stewart (2012); Leinhardt et al. (2015) for example includes scaling laws for both hit-and-run survivors, and also an approximate treatment of changes in (bulk) composition, but appears to be unable to accurately predict collisional water transfer and loss (Burger et al. 2018). It would certainly be a formidable task to derive reasonably precise scaling relations for these subtle processes. The drawbacks of our approach are increased computational expense and complexity. Computational limitations make it necessary to allow at most two post-collision bodies (or gravitationally bound aggregates in our model), and to ignore remaining debris. The error caused by this approach is however limited by the fact that gravity-dominated collisions among similar-sized objects typically result in either 0, 1, or 2 large post-collision bodies along with orders-of-magnitude smaller debris, but practically never in three or more. We are mostly interested in water transport and it may be safe to assume that the water content of these debris fragments will be frequently vaporized after a sufficiently energetic collision, and subsequently lost quickly from small parent bodies (e.g., Odert et al. 2018). Closer to the star (as around 1 au), exposed water on small collision fragments is also not long-time stable. The refractory component of the debris (silicates and iron) on the other hand, can be assumed to further participate in the evolution of the system, and neglecting this material is one of the shortcomings of our model in its current state. Such “losses” from all ~150–200 collisions in each regular scenario typically amount to ~1–1.5 M_⊕ – mostly mantle material with only a relatively small iron fraction of ~0.1 (see Table 2 and Sect. 5.2). Besides increasing the level of dynamical friction and thus eccentricity damping (Genda et al. 2017), it might be accreted at some later point and/or lead to further collisional water losses, even though these small bodies can be easily excited and may often be removed from the system.

We applied two distinct models for planetesimal–planetesimal collisions – either full collisional interaction (pp1), or none at all (pp0). The main motivation for performing also pp0 simulations was to investigate whether the large masses (and low numbers) of planetesimals influence especially the collisional evolution in a priori unknown ways. An additional reason was to increase computational speed, especially by omitting close encounters among planetesimals with pp0. However, it turned out that those gains are largely offset by slower decrease in the number of bodies due to lack of planetesimal–planetesimal accretion (cf. Fig. 9) or removal due to mutual scattering (Raymond et al. 2006). For the most part, our results indicate no obvious significant differences in the main outcome figures between pp1 and pp0 scenarios, even though it is often difficult to distinguish systematic from stochastic variations. One exception is the longer prevalence of a sizeable planetesimal population with pp0, while with pp1 accretion into embryo-sized or larger bodies proceeds faster (visible, e.g., in Fig. 22, Sect. 6.5). Collision statistics also indicate that events with the highest v_imp ∕v_esc (the most destructive) occur in pp1 scenarios (Fig. 12), presumably planetesimal–planetesimal collisions which are not possible with pp0. Those events also dominate the region of highest collisional water losses (Fig. 13), but could also be collisions of rather dry material in the inner disk, where orbital dynamics make them more likely. In this case the influence on final water budgets would actually be very low, but perhaps not on mantle stripping (see Sect. 7.1). In principle full planetesimal–planetesimal interaction (pp1) appears preferable, but the gain in realism from including their collisional evolution has to be considered in light of regularly unrealistically high masses. In addition, treatment with pp1 also includes dynamical planetesimal–planetesimal interaction (at least in close encounters), which was found to be an important part of general dynamical friction (Raymond et al. 2006). However, to identify the influence of these subtle effects a better resolved planetesimal population and better statistics (more simulations) would be required.

7.5 Accuracy of collision simulations

The question of how accurately our SPH simulations predict collision outcomes has two distinct aspects. First, there is the uncertainty of the actual physical state of water reservoirs immediately before the collision, which would require detailed modeling of the internal, atmospheric, and environmental evolution (and coupling between them) for all bodies in the intervals between collisions (see Sect. 7.6). Although a formidable task, we believe that including at least the most influential of these processes should be a future direction to extend our model. Due to the lack of this information, water is modeled as strengthless (surface) fluid, which we believe is the best compromise and is intended to represent roughly mean material properties between the limits of solid ice and hot steam. Burger & Schäfer (2017) and Burger et al. (2018) studied differences in collision outcomes between purely hydrodynamical and full solid-body material models, including icy bodies, from planetesimal- to planet-sized. For the mass range relevant here, these authors found considerable differences in water transfer and loss up to approximately Mars-sized bodies, beyond which strength effects become negligible. Final WMFs and losses differed by less than 20%, but deviations of up to 50% were found for the amount of transferred water from projectile to target and vice versa. Including material strength resulted in consistently lower transfer rates, presumably due to the greater amount of energy required to overcome tensile and shear forces. Differences in transfer and loss for vaporized water in the atmosphere versus liquid or solid reservoirs on or below the surface have not been studied so far to our knowledge. Without considering any surface volatiles in particular, Jutzi (2015) and Emsenhuber et al. (2018) also found that rheology can affect collision outcomes up to approximately Mars-sized bodies, where for this size scale however, only subtle effects (for our purposes) like differences in the post-impact temperature distribution were observed.

Second, accuracy of collision outcomes depends on the numerical resolution (see, e.g., Genda et al. 2015), where our approach naturally allows only relatively small SPH particle numbers (between 25 k and 75 k, see Sect. 3.2). Burger et al. (2018) studied resolution dependence in planet-sized collisions of water-rich bodies, including changes in WMF of the largest and second-largest survivor and also water transfer from projectile to target and vice versa, and thus all important aspects relevant to this work. They found that results for all these quantities were within 25% of approximate resolution convergence, and below 10% for the post-collision WMF of the target body alone. One issue directly connected to limited resolution is the modeling of low WMF, which we accomplish by increasing the water layer thickness to a minimum of approximately three SPH particles during the collision simulation, before its outcome is then corrected for this artificial increase (Sect. 3.2). The error estimated for this procedure is also ~25% or less, including water transfer in hit-and-run (Burger et al. 2018). These combined uncertainties could in principle mean that individual SPH results for water transfer and loss are relatively rough estimates. However, we believe that in the majority of cases outcomes are fairly robust, and that typically large collision numbers combined with the stochastic nature of accretion tracks means that global results and final planets are too.

Part of our numerical scheme is to run SPH simulations for considerably longer than the immediate collision aftermaths in order for all post-collision fragments to form and clearly separate. While this works as intended for the vast majority, “graze-and-merge” encounters are particularly demanding, where the colliding bodies loose sufficient energy in a gentle hit-and-run to remain gravitationally bound, and usually immediately re-collide (after one orbit about each other). The time up to this re-collision depends on the specific geometry and is hard to predict. Therefore, in a few cases SPH simulations may stop before re-collision, or not continue long enough for the final body to settle. The latter happens to some degree in the collision illustrated in Fig. 16, where towards the end a tail of material forms (due to fast rotation induced by the re-collision) whose fate is not simulated anymore, but has certainly only minor consequences there. If the simulation stops while the bodies are still approaching re-collision, they are simply merged (because they are gravitationally bound). However, since the re-collision is necessarily rather gentle (below v_esc), the error can be expected to be small in general.

7.6 The physical state of water and further losses

In our model we currently combine long-term dynamics and realistic collision outcomes, but these are only two aspects of the evolution of water inventories. Energy from short-lived radioactive isotopes, gravitational energy released by differentiation, and large collisions can all melt outer silicate layers, leading to local or even global magma oceans. Their evolution and lifetime is strongly coupled to the overlying (steam) atmosphere which can in turn be heavily modified by volatile outgassing. In addition to direct impact erosion (Genda & Abe 2003, 2005; de Niem et al. 2012; Schlichting et al. 2015; Schlichting & Mukhopadhyay 2018), vaporized water in an atmosphere is susceptible to hydrodynamic loss (escape of hydrogen from dissociated H₂ O molecules) driven by high EUV radiation of the young star (strongly depending, e.g., on its rotation state, see Odert et al. 2018, and references therein). However, it was found that impact-induced magma oceans might cool and solidify again in ~1 Myr and that subsequent condensation of a dense (outgassed) steam atmosphere and ocean formation can also happen on similar timescales (Lebrun et al. 2013), mostly shorter than typical intervals between giant impacts. Hamano et al. (2013) also found rapid (a few Myr) solidification and ocean formation beyond a critical distance from the star – about 0.7 au for the Sun. This is fast enough for hydrodynamic escape to play only a minor role. However, their results also indicate that for (proto)planets closer in a very long magma ocean phase (~100 Myr) can lead to loss of the majority of water via hydrodynamic escape, independent of further collisions. Salvador et al. (2017) found that degassing and solidification of a post-impact magma ocean and subsequent water ocean formation also depend strongly on the initial H₂ O (and CO₂) contents, and that rapid ocean formation is possible on protoplanets similar to early Earth, and possibly also early Venus. While the bulk of such work was done for basically already fully grown planets (typically Earth and Venus), Odert et al. (2018) studied magma ocean evolution and atmospheric losses for Mars-sized embryos at different orbital distances, atmosphere masses, and (EUV) radiation environments. The results of these latter authors indicate that at the orbit of Mars the time for steam-atmosphere condensation (e.g., after a giant collision) is very small compared to the timescale of atmosphericescape, and almost all post-collision water content can be retained. At Earth’s orbit, typically a relatively small fraction of water is lost before it can condense into oceans, while at the orbit of Venus virtually all atmosphere is expected to be lost before condensation is able to stop it for these small bodies. For even lower masses (≲Moon-mass) – like the collisional debris in our SPH outcomes – not only EUV-driven hydrodynamic escape but also thermal (Jeans) escape becomes important, which plays only a minor role for Mars-sized bodies and is negligible for Earth-mass objects (e.g., Odert et al. 2018). Depending on their (post-collision) temperature such small bodies can loose their atmospheres and water very rapidly. In addition, the evolution after a large collision event is connected to the current (dynamical) environment, where continued energy influx by impacting bodies for example can also delay surface solidification and water condensation (Maindl et al. 2015).

Our results indicate that, at least for potentially habitable planets, water delivery (by the few decisive events) happens typically late in their accretion histories once they have grown to at least 50% and often to over 90% of their final mass (Sect. 6.5). This increases the chances of water retention due to their already high gravity and also because fewer further (giant) collisions lie ahead. We note however, that this means only the final water-delivering events, where the (mostly smaller) impactors themselves may have complicated collision histories of their own. The (absolute) times when the bulk of water is delivered can vary strongly, with most between 10 and 100 Myr (Table 3) when intervals between large collisions are typically already several Myr (Sect. 6.5 and especially Fig. 22), and in some cases also considerably later with typical intervals >10 Myr. Therefore, repeated ocean formation between (giant) collisions seems to be a plausible scenario, at least for growing planets in the (Sun’s) habitable zone and beyond. The prevalence of hit-and-run adds an additional aspect to these considerations. Depending mainly on the impact angle, energy deposition can be greatly reduced compared to central (head-on) collisions, with possibly profound consequences for water vaporization (Burger et al. 2018) and magma ocean formation (or the lack thereof) in ways that remain mostly unstudied.

7.7 Future work

Besides the various foreseeable difficulties with other possible methods, for example reliably predicting all the subtleties of hit-and-run with scaling laws or machine learning techniques, we chose this direct computational approach because we believe that it is a solid basis for including further details and physical processes in future work. The simulations presented here follow only the dynamical and (immediate) collisional evolution of water reservoirs, while at the end of the day it is the coupled dynamical, collisional, geophysical, and atmospherical evolution that shapes a planet’s water inventory. In that sense our current results are only upper limits for final water contents, with the modification roughly decreasing with increasing orbital distance, where collisions are less energetic (due to lower relative velocities) and magma ocean phases are shorter (due to lower stellar flux). One possibility to incorporate further modeling without having to repeat whole simulations is to use already existing accretion histories of planets of interest as a backbone for modeling of additional physics and/or higher SPH resolution. This could include only a planet’s direct accretion track, or those of all impactors (and bodies impacting them, etc.) as well. To that end, we also offer precise data on the dynamical and collisional evolution of all planets formed in our scenarios to interested researchers (see footnote ⁶, page 17).

A shortcoming of the simulations in this study is the relatively poorly resolved planetesimal population and we are already running higher-resolution simulations with half the initial mass in (significantly more) planetesimals. In addition, embryo formation happens on different timescales depending on location, and presumably proceeds considerably faster in the inner part of the disk (Kokubo & Ida 2002; McNeil et al. 2005). More realistic initial conditions beyond equal planetesimal masses and uniform embryo formation throughout the whole disk could be obtained from simulating the collisional evolution of a large number of small bodies over a relatively brief period with our hybrid framework (see, e.g., Raymond et al. 2006; Bonsor et al. 2015; Carter et al. 2015). Alternatively, more realistic initial conditions and mass distribution can be based on semi-analytical models for earlier growth phases (e.g., McNeil et al. 2005).

The first application in this study is in the classical Solar System formation environment, including giant planet formation without migration, but the diversity of detected extrasolar planets shows that this is only one possible path for planetary accretion. We believe that realistic collision modeling can provide a better understanding of terrestrial planet formation in the great diversity of planetary systems, not least to constrain water contents of detected exoplanets for future observations and assessment of their principal habitability.

8 Summary

We have developed a new hybrid framework to combine the dynamical and collisional evolution of a debris disk, comprising planetary embryos and planetesimals, over several hundred Myr of late-stage terrestrial planet formation. Each physical collision is evaluated by a dedicated 3D SPH simulation of differentiated, self-gravitating bodies, the results of which are used self-consistently on the fly in the further N-body computation. The main focus is on collisional water transport, especially to potentially habitable planets, with volatile material like water being particularly susceptible to transfer and loss processes. Collisions can result in either zero (very destructive), one (accretion or erosion), or two (hit-and-run) post-collision bodies, which allows us to fully track water transfer and loss in hit-and-run encounters as well.

In this first application we chose a Solar System-like architecture, including already formed nonmigrating giant planets and a debris disk between 0.5 and 4 au. This well-studied setting allows us to compare and validate our results against existing studies, where at least for water delivery perfect merging (PM) was assumed almost exclusively. However, our main purpose is not to reproduce the inner Solar System in detail, but to study the importance and implications of realistic modeling of collisional water transport for the early Solar System, and also for planet formation in extrasolar systems. Results of our PM comparison runs are in general agreement with existing studies, while realistic collision treatment typically leads to lower mass planets with considerably reduced water inventories (Fig. 2). To our knowledge this is the first time that collisional water transport has been modeled self-consistently in combination with the dynamical evolution up to final planets.

Our results have revealed several key implications for water transport, and also for terrestrial planet formation in general:

• We find that the accretion phase is truly lengthened by up to a factor of two with our collision model compared to PM runs, often to 100–200 Myr or more (Sects. 5.4 and 6.5). This is signified by slower mass growth (Fig. 15), considerably slower decrease in the number of bodies, also for embryo-sized or larger bodies alone (Figs. 9 and 10), and significantly later final giant collisions (involving embryos or larger bodies).

• As also found in previous PM studies, eccentric giant planets (eJS) are detrimental to water delivery due to fast removal of water-rich material from the outer disk by strong resonances, while with circular giant planets (cJS) much more water-rich material initially remains (see Fig. 6). With our model however, the combination of considerable collisional losses and more ejections of water-rich bodies (presumably due to longer accretion times) strongly reduces final water contents of cJS planets as well, and brings them much closer to those of eJS scenarios – and Earth-like water contents (Sect. 5.3).

• Due to the quick removal of most water-rich material from systems with eJS architectures, water delivery remains very stochastic even with our realistic collision model, where the (chaotic) fate of individual water-rich embryos often becomes decisive. With cJS water delivery is more robust, since usually a much higher number of water-rich bodies survives and slowly diffuses inwards.

• We studied water delivery and associated accretion histories of potentially habitable planets (0.75–1.5 au) in particular detail and found that only very few – often one or two – decisive collisions deliver the vast majority of a planet’s final inventory (Sects. 6.2, 6.3, 7.2 and Fig. 15), at least if embryos form initially throughout the whole disk and with a relatively small planetesimal fraction of ~0.3. These encounters are predominantly with embryo-sized or larger bodies and only rarely with water-rich planetesimals, which either get removed or incorporated into the delivering embryos and protoplanets first. Since only a few collisions (accretionary and hit-and-run alike) seem to clearly dominate final water delivery, it naturally becomes an intrinsically stochastic process.

• Sixteen of the 18 potentially habitable planets formed with our model accrete from dry material in the inner disk before significant water delivery happens, typically after the planet has reached at least 50%, often already 90%, of its final mass (Sect. 6.5). Already long average intervals between giant impacts (several up to tens of Myr) then potentially result in repeated ocean formation between all-vaporizing collisions. The two remaining planets start as water-rich embryos in the outer disk and end up as very water-rich “ocean planets” around 1 au (see Fig. 15 and Sect. 7.2).

• Water losses in individual collisions are relatively low on average, mostly below 20–30% (Fig. 13, Sect. 5.5). Encounters in the direct accretion histories of (potentially habitable) planets generally lead to lower values mostly below 10%, but can nevertheless result in considerable cumulative losses. Growing planets in our scenarios experience typically 20–30 collisions overall, and still around 10 encounters involving only embryo-sized or larger bodies (Sect. 6.4). Collisional erosion clearly preferentially strips water and upper mantle material.

• Final water inventories are typically reduced by a factor of two or more compared to simple PM. This conclusion is based not only on PM comparison runs, but also on cumulative collision losses (Fig. 20) and a “synthetic PM” approach, computed from results with our collision model (Sect. 6.3).

• Hit-and-run is a frequent outcome during late-stage accretion with occurrence rates around 50%, consistent with earlier studies (Sects. 5.5 and 6.4). This has strong implications for water transport to growing planets, since often only a small fraction of the impactor’s water budget (and mass in general) is transferred to the target body, especially compared to more head-on, accretionary collisions (cf. Fig. 16).

• Our results indicate that (hit-and-run) interactions among those protoplanets that eventually grow into planets are rare, and their final growth phases are mostly rather isolated (from each other). Even though hit-and-run encounters with larger bodies can be very transformative (for the smaller target), this also actually happens only rarely to growing planets in a realistic dynamical environment (see Sect. 7.1).

• The often claimed justification for PM that hit-and-run survivors mostly re-collide (and merge) soon after the initial encounter is not supported by our results and occurs only in the minority of cases (Sect. 6.4). Therefore, it is crucial to include proper treatment of water transfer and loss for hit-and-run, and to track the further independent evolution of both survivors.

• An interesting additional result is the high formation rate of small, close-in, iron-rich planets resembling Mercury (Fig. 4), which seem to robustly form in a long sequence of high-velocity mantle stripping collisions instead of one or a few chance events (see Sect. 7.1). However, our framework was primarily designed to model collisional water transport, where considerable refractory debris is often neglected.

Acknowledgements

We thank John Chambers for many helpful comments that improved our manuscript. Simulations in this paper made use of the REBOUND code which can be downloaded freely at http://github.com/hannorein/rebound. C.B. and A.B. acknowledge support by the FWF Austrian Science Fund projects S11603-N16 and S11608-N16. The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no. INST 37/935-1 FUGG.

References

All Tables

All Figures

Fig. 1 Collision geometry of two differentiated bodies in a frame centered on the target. The impact angle α ranges from 0° for head-on to 90° for just grazing. The impact velocity is |v0|. The common convention is to label the larger body as the target, but in some situations the smaller object is the one followed (e.g., on its growth path towards a planet) and is thus considered (and referred to as) the target.

In the text

Fig. 2 Final systems for all regular scenarios (cf. Table 1). Each column represents a particular giant planet architecture (cJS or eJS) and mode of planetesimal–planetesimal interaction (pp1 or pp0). Upper five rows: results based on the SPH collision model, lower three: PM comparison runs. The sizes of planets and their iron cores (black) are proportional to mass1∕3. The region corresponding to our definition of potentially habitable planets (see Sect. 6.1) is indicated as a gray bar, and the Solar System planets are plotted as open circles as a reference.

In the text

	Fig. 3 Final systems for all embryo-only scenarios (cf. Table 1). The columns represent cJS and eJS architectures, while planetesimal–planetesimal interactions are always included (pp1). Upper five rows: results based on the SPH collision model, lower three: from PM comparison runs. Colors, symbols, and their sizes have the same meanings as in Fig. 2.
In the text

	Fig. 4 Various properties of final planets in the regular scenarios, as a function of semi-major axis. Mass and core mass fraction (first and last panel) are plotted linearly, while water mass fraction and total water mass (middle panels) are log-plots. The region of potentially habitable planets (see Sect. 6.1) is indicated as a gray bar.
In the text

	Fig. 5 Time evolution of the total mass (upper panel) and total water mass (lower panel) for all regular scenarios.
In the text

Fig. 6 Simulation snapshots in the a–e plane for three different runs (columns). Typical evolution with eJS in the left column (SPH15-eJS-pp1) vs. cJS in the middle (SPH06-cJS-pp0) vs. cJS with PM on the right (PM06-cJS-pp0). Symbol sizes are proportional to mass1∕3. Each row represents an equal time in all scenarios: up to 1 Myr the evolution appears similar, except some more dynamical excitation with eJS (left). Within only ~5 Myr however, the great majority of water-rich material is removed with eJS, but remains in the system with cJS and starts to diffuse inwards. At 50 Myr, accretion is still actively ongoing with the SPH collision model (even though less so with eJS in this particular scenario), while with PM (right) it is almost over, except for some remaining planetesimals. We note that the cJS scenarios (middle and right columns) are pp0 runs, where some more leftover planetesimals are common. After ~100 Myr, body numbers are declining in the SPH scenarios (left and middle columns), but accretion and associated collisional evolution is still not over. We also note the damping effect on the eccentricity of the giant planet with eJS, consistent with similar studies (e.g., O’Brien et al. 2006). These scenarios represent typical behavior and were not specifically selected, except that the left and middle runs are those where the potentially habitable planets in Figs. 18 and 17 form.

In the text

	Fig. 7 Time evolution of the mass in four semi-major axis bins for all regular scenarios. The curves have been time-averaged (smoothed) over ~3 Myr.
In the text

	Fig. 8 Time evolution of the water mass (in O⊕) in four semi-major axis bins for all regular scenarios. We note the log y-axis in the upper two panels. The curves have been time-averaged (smoothed) over ~3 Myr.
In the text

	Fig. 9 Total number of bodies (embryos or protoplanets + planetesimals) over time for all regular scenarios.
In the text

	Fig. 10 Number of embryos + protoplanets (i.e., excluding planetesimals) over time for all regular scenarios, split into cJS (upper panel) and eJS (lower panel) for better visibility. The respective “PM” curves have a slight y-offset of −0.2 for the samereason.
In the text

	Fig. 11 Impact velocity in units of vesc plotted vs. impact angle for all collisions in the regular scenarios with our collision model. Encounters exclusively among embryos and protoplanets (i.e., excluding planetesimals) are highlighted as larger symbols.
In the text

	Fig. 12 Impact velocity in units of vesc over time for all regular scenarios with our collision model. Larger symbols mark collisions exclusively among embryos and protoplanets.
In the text

	Fig. 13 Water loss in individual collisions over time for all regular scenarios with our model. Color coding is identical to Figs. 11 and 12. Water loss is defined as the fraction of the (combined) pre-collision water inventory that is not incorporated in either of the post-collision bodies. Larger symbols mark collisions exclusively among embryos and protoplanets.
In the text

	Fig. 14 All potentially habitable planets formed in the regular scenarios in the water-inventory vs. mass space. We note the log y-axis.
In the text

	Fig. 15 Growth curves of all potentially habitable planets formed in the regular scenarios, with water inventories in the top panel (note the log y-axis) and mass in the bottom panel. The precise definition according to which evolutionary tracks are followed (especially for hit-and-run) is discussed in the text.
In the text

Fig. 16 Snapshots of SPH collision simulations of two different scenarios (top and bottom row). The impactor (smaller body, in blue) approaches the target from the left. The color coding shows water on the impactor in blue and silicate material (basalt) in red. (Iron core material is hardly exposed and is light gray.) The target body was actually also simulated with a water shell, but is marked in red too in order to highlight water transfer from a water-rich impactor to a dry target. Top row: accretionary collision (more precisely “graze-and-merge”, where after a brief hit-and-run-like period the bodies fall back and merge), with vimp ∕vesc = 1.2, impact angleα = 28° and impactor-to-target mass ratio γ ~1:2. Bottom row: typical hit-and-run encounter with only relatively little water transfer, where vimp ∕vesc = 1.8, α = 47° and γ ~1:6. The collisionin the top panels is the last (at ~70 Myr) decisive water-delivering event for the planet in Fig. 17, the one in the bottom panels shows the first (at ~30 Myr) of two such decisive collisions for the planet in Fig. 18 (see Sect. 6.3 for details).

In the text

Fig. 17 Accretion history of the potentially habitable planet formed in run SPH06-cJS-pp0 (marked with superscript 1 in Tables 3 and 4, and also illustrated in Fig. 6). It forms at ~1.1 au, with 1.2 M⊕ and 24 O⊕ of water. Besides the actual results with our SPH collision model (solid lines), hypothetical growth curves with synthetic PM, and also with added cumulative losses are plotted – both defined and discussed in the text (Sect. 6.3). The crosses mark the respective property of the other (second) body in collisions, besides the traced planet.

In the text

	Fig. 18 Accretion history of the potentially habitable planet formed in run SPH15-eJS-pp1 (marked with superscript 2 in Tables 3 and 4, and also illustrated in Fig. 6). It forms at ~0.9 au, with 1.03 M⊕ and 1.59 O⊕ of water. Lines and symbols have the same meanings as in Fig. 17.
In the text

	Fig. 19 Water losses in individual collisions over time during accretion of all 18 potentially habitable planets in the regular scenarios. Water loss is defined as the fraction of the (combined) pre-collision water inventory that is not incorporated into either of the post-collision bodies.
In the text

	Fig. 20 Cumulative water losses (from their individual collision histories) vs. final water inventory for all 18 potentially habitable planets formed with our model in the regular scenarios. The diagonal line indicates cumulative losses equal to the final water content.
In the text

Fig. 21 Water accretion efficiency Ξw in individual collisions over time during accretion of all 18 potentially habitable planets in the regular scenarios. Positive values indicate an increase in water inventory (in units of the projectile’s water content), while negative ones indicate a decrease (in units of the target’s water content). See the text for a precise definition of Ξw and further discussion.

In the text

Fig. 22 Time since the last collision (evaluated each time a collision happens) from the perspective of a growing planet as a function of time, for all 18 potentially habitable planets in the regular scenarios with our model. Small crosses count all collisions, large ones only impacts by embryos and protoplanets (i.e., excluding planetesimals). Each collision is assigned two data points in general – a small and a large cross – and if the last impacting body was an embryo or protoplanet then they coincide. The continuous and dashed lines indicate elapsed time (since the last collision): since simulation start (upper line), 1/10 of that time (middle line), and 1/100 of that time (lower line).

In the text

Fig. 23 Water losses in individual collisions as a function of the projectile-to-target mass ratio (where the target is always the larger body), for all regular scenarios with our model. Water loss is defined as the fraction of the (combined) pre-collisionwater inventory that is not incorporated into either of the post-collision bodies. Larger symbols indicate collisions among embryos and protoplanets alone (i.e., excluding planetesimals). Slightly negative values are merely numerical artifacts and can be assumed to be zero.

In the text