COMPOSITIONAL EVOLUTION DURING ROCKY PROTOPLANET ACCRETION

In this work we present simulations initiated with either 10,000 (low resolution) or 100,000 (high resolution) particles arranged in an annulus, extending from 0.5 AU to either 1.5 or 3 AU, around a 1 solar mass star.

The initial disk of planetesimals is distributed according to a surface density, Σ = Σ₁ a^−3/2, similar to that of the Minimum Mass Solar Nebula (MMSN; Weidenschilling 1977; Hayashi 1981), using a surface density at 1 AU of Σ₁ = 10 g cm⁻² to match previous work (e.g., Kokubo & Ida 2002; Leinhardt et al. 2009; Paper I).

We begin our simulations with planetesimals of 15 different sizes, evenly spaced in mass, m, and chosen according to the cumulative number distribution:

$$\begin{eqnarray}&&{dN}={m}^{-p}{dm},\end{eqnarray} \tag{ 1 }$$

with p = 2.5, a slope appropriate for the runaway growth phase (Kokubo & Ida 1996). This initial mass distribution (black line in Figure 1) is similar to that used in Paper I (in which the size distribution was generated by running PKDGRAV using perfect merging with a disk with five times as many equal-size planetesimals). The simulations begin with most of the planetesimals having radii of a few hundred kilometers in accordance with the mass distribution shown in Figure 1.

Zoom In Zoom Out Reset image size

In this work we study the intermediate stage of planet formation; sizeable planetesimals have already formed and will continue to grow into embryos via runaway accretion during the simulations. As it is known that planetesimals can differentiate in the first few million years of the solar system's evolution (e.g., Lugmair & Shukolyukov 1998; Srinivasan et al. 1999; Amelin 2008; Kruijer et al. 2014), we assume that all planetesimals are differentiated at the start of our simulation, and we assign the same core fraction to all particles, regardless of their location within the disk. This uniform initial differentiation is probably a simplification, but will aid understanding of the role of collisions on the compositional evolution of growing bodies.

The initial conditions used for the simulations presented in this paper are summarized in Tables 1 and 2.

Table 1.  Initial Conditions for the Calm Disk Simulations

Parameter	Value	Notes
Mass of planetesimal disk	2.5 M${}_{\oplus }$	⋯
Radial extent of disk	0.5–1.5 AU	⋯
Initial number of particles	100,000	High resolution
	10,000	Low resolution
Initial planetesimal mass	$3.2\times {10}^{-11}$ ${M}_{\odot }$–1.5 × 10−8 M⊙	High resolution
	3.2 × 10−10 M⊙–3.6 × 10−8 M⊙	Low resolution
Initial planetesimal radius	196–1530 km	High resolution (f = 1)
	423–2050 km	Low resolution (f = 1)
Initial core fraction	0.22, 0.35	⋯
Planetesimal density	2 g cm−3	⋯
Resolution limit	3 × 10−11 M⊙ (1 × 10−5 ${M}_{\oplus }$)	High resolution
	3 × 10−10 M⊙ (1 × 10−4 ${M}_{\oplus }$)	Low resolution
Time step	0.01 yr, 0.005 yr (with gas drag)	⋯
Duration	600,000 yr (∼20 Myr effective time)	⋯
Initial gas density, ${\rho }_{{\rm{g}},1}$	$1.4\times {10}^{-9}$ g cm−3	MMSN
Gas dissipation start time	56 kyr (2 Myr effective time)	⋯
Gas dissipation timescale	2.8 kyr (100 kyr effective time), $\infty$	⋯

Download table as:  ASCIITypeset image

Table 2.  Initial Conditions for the Grand Tack Simulations

Parameter	Value	Notes
Mass of planetesimal disk	4.85 M${}_{\oplus }$	⋯
Radial extent of disk	0.5–3.0 AU	⋯
Initial number of particles	100,000	High resolution
	10,000	Low resolution
Initial planetesimal mass	6.2 × 10−11 M⊙–2.9 × 10−8 M⊙	High resolution
	6.2 × 10−10 M⊙–7 × 10−8 M⊙	Low resolution
Initial planetesimal radius	245–1900 km	High resolution (f = 1)
	528–2550 km	Low resolution (f = 1)
Initial core fraction	0.22, 0.35	⋯
Planetesimal density	2 g cm−3	⋯
Resolution limit	6 × 10−11 M⊙ (2 × 10−5 ${M}_{\oplus }$)	High resolution
	6 × 10−10 M⊙ (2 × 10−4 ${M}_{\oplus }$)	Low resolution
Time step	0.01 yr, 0.005 yr (with gas drag)	⋯
Duration	600,000 yr (∼20 Myr effective time)	⋯
Jupiter initial location	3.5 AU	⋯
Jupiter minimum semimajor axis	1.5 AU	⋯
Jupiter final location	5.2 AU	⋯
Migration timescale	2.8 kyr (100 kyr effective time), 17 kyr (600 kyr)	Fast, slow migration
Migration start time	56 kyr (2 Myr effective time), 220 kyr (8 Myr)	Early, late migration
Initial gas density, ${\rho }_{{\rm{g}}}$	1.4 × 10−9 g cm−3	MMSN
	see Equations (6) and (7)	Based on hydrodynamic simulations

Download table as:  ASCIITypeset image

N-body simulations of the intermediate phase of planet formation with large N have often employed an expansion of the particle radii in order to increase the collision rate and hence speed up the evolution (e.g., Kokubo & Ida 1996, 2002; Leinhardt & Richardson 2005; Morishima et al. 2008). This expansion factor was discussed in detail by Kokubo & Ida (1996, 2002), in which they determine that it has very little effect on the evolution during this phase, merely altering the timescale for planet formation. In this work we are studying the combined consequences of many collisions of different types, including erosive collisions, and since the expansion factor prevents very close approaches from occurring (which would excite planetesimal velocities), it is expected that our simulations slightly underestimate mantle erosion.

While the number of bodies in the simulation is large, dynamical friction from the planetesimals controls the velocity dispersion, and the expansion factor does not alter their growth mode (Kokubo & Ida 1996). As the number of planetesimals decreases, and massive embryos begin to significantly alter velocities via gravitational scattering, the expansion factor likely causes the results to become less accurate. Once the giant impact stage of planet formation begins, this becomes a significant problem, and so we cannot integrate accurately through the late stages of formation using the expansion factor. Owing to limitations in computing power, it is still not practical to conduct large N simulations covering millions of years of planet formation without using this expansion factor.

In previous work we have used an expansion factor, f = 6 (Leinhardt & Richardson 2005; Leinhardt et al. 2009; Paper I; Leinhardt et al. 2015), and, for consistency, we adopt this same value here. In our simulations giant planets do not have their radii inflated (any planetesimals colliding with giant planets are forced to merge to avoid complications), and gas drag is calculated using the natural radius, not the expanded value.

During the intermediate stages of planet formation we are investigating, the expansion factor causes the evolution to speed up by a factor ∼f². We therefore quote both a simulation time and an effective time corresponding to realistic radii (f = 1).

EDACM (Leinhardt & Stewart 2012; Leinhardt et al. 2015) provides an analytical description of the outcome of a collision between gravity-dominated bodies, like the planetesimals modeled in this work. Using this collision model, the size and velocity distributions of the fragments of each collision are calculated based on the impact parameter and collision velocity, with each collision classified as having one of seven different outcomes: perfect merging, partial accretion, hit-and-run (with the projectile intact, disrupted, or supercatastrophically disrupted), erosion, or supercatastrophic disruption.

Since EDACM can produce many fragments from the original two impactors, it is necessary to impose a resolution limit, a mass below which no further resolved particles will be generated. Any remnants smaller than this limit are treated semi-analytically as "unresolved debris." This debris is distributed in a series of circular annuli each 0.1 AU wide, with the debris placed in the annulus corresponding to the location of the collision that produced it. This material is assumed to have circular Keplerian orbits. As well as being generated by collisions, this unresolved debris is reaccreted by the resolved particles according to their geometric cross section (calculated using inflated radii) and orbit eccentricity (see Leinhardt & Richardson 2005). As was shown by Leinhardt et al. (2015), there is never a significant amount of mass in this unresolved component (see also Section 3.2.2).

In this work we track both the total mass of debris in each annulus and the mass of core material (since this proved to be a problem in Paper I). Additionally, each particle and debris annulus has an "origin histogram," which tracks the fraction of that object's mass originating from within each 0.1 AU annulus (see Leinhardt & Richardson 2005 for details).

Differentiated planetesimals are modeled with core and mantle components, which are tracked independently. In a collision between two such planetesimals core and mantle material may be exchanged, or distributed among fragments. Marcus et al. (2010) analyzed SPH collisions between differentiated bodies that started with separate core and mantle components, tracking the two types of particles to determine the result of the collision. They found that the fractional mass of the iron core of the largest post-collision remnant scaled with impact energy (that is, the more energetic the impact, the higher the core fraction of the largest remnant) and lies between two limiting cases: model 1, in which cores always merge; and model 2, in which core material from the projectile is only accreted if the largest remnant is more massive than the original target (see Marcus et al. 2010; Paper I). As in Paper I, we assume that any core material not accreted by the largest remnant is distributed to the smaller collision remnants (ensuring that mass is conserved), starting by placing as much as possible in the second largest.

In Paper I, the two empirical models from Marcus et al. (2010) describing the mass of the iron core of the largest remnant of a collision were applied to each collision in the simulation in a post-processing step. This proved to be problematic as it could not account for the composition of the unresolved material;⁵  this had to be estimated in order to calculate embryo compositions. To avoid this, we have incorporated the models used in Paper I into our N-body code, which can now track the core mass of both resolved particles and debris annuli directly, applying the mantle stripping model as each collision occurs.

In Paper I very little difference was found between the final core fractions of embryos using the two models from Marcus et al. (2010). We investigate the difference again in this work, and we also define a third model as the mean of the original two, ${M}_{\mathrm{core},\mathrm{lr},3}=({M}_{\mathrm{core},\mathrm{lr},1}+{M}_{\mathrm{core},\mathrm{lr},2})/2$ (see Table 3).

Table 3.  Mantle Stripping Models from Marcus et al. (2010)

Collision Outcome	M${}_{\mathrm{core},\mathrm{lr}}$ Model 1	M${}_{\mathrm{core},\mathrm{lr}}$ Model 2
Perfect merging	M${}_{\mathrm{core},\mathrm{targ}}$ + ${M}_{\mathrm{core},\mathrm{proj}}$	M${}_{\mathrm{core},\mathrm{targ}}$ + ${M}_{\mathrm{core},\mathrm{proj}}$
Partial accretion	$\mathrm{min}$(Mlr, ${M}_{\mathrm{core},\mathrm{targ}}$ + ${M}_{\mathrm{core},\mathrm{proj}}$)	M${}_{\mathrm{core},\mathrm{targ}}$ + $\mathrm{min}$(M${}_{\mathrm{core},\mathrm{proj}},$ Mlr−${M}_{\mathrm{core},\mathrm{targ}}$)
Hit-and-run	M${}_{\mathrm{core},\mathrm{targ}}$	M${}_{\mathrm{core},\mathrm{targ}}$
Erosive, supercatastrophic disruption	$\mathrm{min}$(Mlr, ${M}_{\mathrm{core},\mathrm{targ}}$ + ${M}_{\mathrm{core},\mathrm{proj}}$)	$\mathrm{min}$(Mlr, ${M}_{\mathrm{core},\mathrm{targ}}$)

Note. Model 3 is defined as the average of the original two: ${M}_{\mathrm{core},\mathrm{lr},3}=({M}_{\mathrm{core},\mathrm{lr},1}+{M}_{\mathrm{core},\mathrm{lr},2})/2$.

Download table as:  ASCIITypeset image

Collisions are not the only process that can effect bulk composition and core–mantle ratio of planetesimals; the initial oxidation state may have varied across the terrestrial planet feeding zone, but we do not know what any such distribution would have been. For this reason, and so that we will see the effects due to collisions alone, we assume the same initial core fraction across the terrestrial planet region. We explore two cases for this initial core fraction (for both planetesimals and unresolved debris) matching the values used in Paper I: 0.35, representative of high iron content chondrites (EH), and 0.22, representative of low iron content chondrites (LL).

Core and mantle materials are always isolated in our model, and we do not allow any mixing or reequilibration to occur. Accreted core material is assumed to merge instantly with the core of the target since we expect planets to be (at least partially) molten during accretion. This is a significant simplification for some proposed scenarios (in which ongoing equilibration between metal and silicate alters the metal–silicate ratios; e.g., Rubie et al. 2015), but a reasonable one given that the timescales involved and the degree to which mixing would occur are still poorly understood (e.g., Dahl & Stevenson 2010).

Walsh et al. (2011) presented the Grand Tack model of the evolution of the solar system, in which the inner planetesimal disk is truncated by the inward-then-outward migration of Jupiter and Saturn. In this scenario, Jupiter forms first, creating a gap in the disk, allowing it to migrate inward via type II migration. Saturn grew more slowly and later began migrating inward more rapidly than Jupiter, allowing Saturn to "catch up" with Jupiter, becoming trapped in the 2:3 mean motion resonance. This occurs as Jupiter reaches 1.5 AU; the gaps in the gas disk opened by the gas giants then overlap, changing the balance of the torques, and causing the migration direction to reverse. Jupiter and Saturn then migrate outward together as the gas disk dissipates, leaving the gas giants close to their present-day orbits and the outer solar system in the appropriate configuration for a late instability (associated with the late heavy bombardment). Jupiter migrating inward to 1.5 AU before "tacking" truncates the inner disk at ∼1 AU, producing the much sought after small Mars (Walsh et al. 2011; O'Brien et al. 2014).

In the simulations that employ the Grand Tack we include Jupiter as an N-body particle and model Jupiter's migration by giving it an additional acceleration to match the velocity changes described by Walsh et al. (2011). Unlike previous Grand Tack simulations (e.g., Walsh et al. 2011; O'Brien et al. 2014; Rubie et al. 2015), which begin with a possibly unrealistic bimodal distribution of planetesimals and embryos and allow Jupiter to begin migrating at the start, we begin with a less evolved planetesimal disk and allow embryos to begin growing via collisions while Jupiter remains at 3.5 AU for several Myr before it begins migrating.

In this work we test the effect of Jupiter's inward-then-outward migration starting at two different times (since Jupiter's formation time and the gas disk lifetime are uncertain): the commonly quoted time of 3 Myr (56 Kyr simulation time), and a "late" Grand Tack in which the migration begins 9 Myr (220 kyr simulation time) after the birth of the solar nebula. We also test two timescales for migration, both the standard 100 kyr (2.8 kyr simulation time) value of Walsh et al. (2011) and a "slow" migration timescale of 600 kyr (17 kyr simulation time). Using the expansion factor method makes matching the timescales problematic; we have opted to accelerate the migration by the same f² factor by which the accretion timescale is reduced. Comparison of the slow and standard timescales will show the effect of slower migration, but should also reveal any effect of the accelerated evolution relative to the migration caused by the expansion factor (the slow timescale is equivalent to the nominal migration rate accelerated only by a factor f).

The acceleration, ${\dot{v}}_{\mathrm{in}(\mathrm{out})},$ given to Jupiter to model its inward (outward) migration is described by

$$\begin{eqnarray}&&{\dot{v}}_{\mathrm{in}}=\displaystyle \frac{{v}_{3.5\ \mathrm{AU}}-{v}_{1.5\ \mathrm{AU}}}{\tau },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\dot{v}}_{\mathrm{out}}=\displaystyle \frac{{v}_{1.5\ \mathrm{AU}}-{v}_{5.2\ \mathrm{AU}}}{\tau }\;{e}^{(t-{t}_{\mathrm{tack}})/\tau },\end{eqnarray} \tag{ 3 }$$

where ${v}_{1.5\ \mathrm{AU}}$ is Jupiter's Keplerian velocity at 1.5 AU—the location at which it "tacks," t is the time since the start of the simulation, t_tack is the time at which Jupiter's migration reverses direction, and τ is the migration timescale. Jupiter begins the simulation at 3.5 AU, reaches an innermost semimajor axis of 1.5 AU, and ends the simulation after the Grand Tack at ∼5.2 AU. Jupiter feels the gravitational influence of the planetesimals just as any planetesimal or embryo; however, any object that collides with Jupiter is forced to merge.

We do not include Saturn or the ice giants in this work as we expect them to have no significant effect on the terrestrial planet region, as discussed by Walsh et al. (2011). We also ignore the outer planetesimals originally situated beyond the orbit of Jupiter, which would be numerically very expensive to model at these resolutions, as it is not expected that a significant number would be scattered into the terrestrial region (∼3% of final planet mass; O'Brien et al. 2014; Rubie et al. 2015), and most would be accreted after the end time of our simulations.

In the calm disk scenario the initial planetesimal disk extends from 0.5 to 1.5 AU, with no giant planets included.

In this paper we present some simulations that include a gas density distribution within the protoplanetary disk (the effects of gas drag were not included in Paper I). In these simulations we calculate the drag force on the planetesimals following the prescription of Adachi et al. (1976) and Brasser et al. (2007).

The drag force is applied to each body in the simulation as an additional acceleration in the direction opposite to that of the body's relative motion through the gas,

$$\begin{eqnarray}&&{\dot{v}}_{{\rm{D}}}=-\displaystyle \frac{1}{2\;m}{C}_{{\rm{D}}}\pi {r}^{2}{\rho }_{{\rm{g}}}{({\boldsymbol{v}}-{{\boldsymbol{v}}}_{{\rm{g}}})}^{2}\end{eqnarray} \tag{ 4 }$$

where m and r are the planetesimal mass and radius, respectively, C_D is the dimensionless drag coefficient, ${\rho }_{{\rm{g}}}$ is the gas density, and ${\boldsymbol{v}}-{{\boldsymbol{v}}}_{{\rm{g}}}$ is the velocity of the planetesimal relative to the gas velocity. The drag acceleration is calculated using the uninflated planetesimal radius and is applied only to resolved particles; unresolved material is assumed to accrete before its orbit is significantly affected by drag, or to be replenished by inward drift from larger radii.

The gas profile we adopt has one of two forms, either a power law matching the surface density of the solar nebula, or a third-order polynomial obtained from a fit to the profile used by Walsh et al. (2011), which matches the results of hydrodynamic simulations from Morbidelli & Crida (2007). The MMSN gas density is given by (Brasser et al. 2007)

$$\begin{eqnarray}&&{\rho }_{{\rm{g}}}={\rho }_{{\rm{g}},1}\;{s}^{-11/4}{e}^{-{z}^{2}/{z}_{s}^{2}},\end{eqnarray} \tag{ 5 }$$

where ${\rho }_{{\rm{g}},1}$ is the initial gas density at 1 AU, s is the distance from the central star in the (x, y)-plane, and z_s = (0.047s^5/4) is the scale height of the gas disk, and the hydrodynamically derived gas density is given by

$$\begin{eqnarray}{\rho }_{{\rm{g}}} & = & \;\displaystyle \frac{5.2}{{s}_{{\rm{J}}}}\left(21880{\left(\displaystyle \frac{s}{{s}_{{\rm{J}}}}\right)}^{3}\;-\;48010{\left(\displaystyle \frac{s}{{s}_{{\rm{J}}}}\right)}^{2}\right.\\ & & \left.+\;25420\displaystyle \frac{s}{{s}_{{\rm{J}}}}\;+\;663.3\right)\displaystyle \frac{{e}^{-{z}^{2}/{z}_{s}^{2}}}{\sqrt{\pi }{z}_{s}},\end{eqnarray} \tag{ 6 }$$

if $s\lt =\;0.9473\;{s}_{{\rm{J}}},$ or,

$$\begin{eqnarray}&&{\rho }_{{\rm{g}}}=\displaystyle \frac{5.2}{{s}_{{\rm{J}}}}\;90\;\displaystyle \frac{{e}^{-{z}^{2}/{z}_{s}^{2}}}{\sqrt{\pi }{z}_{s}},\end{eqnarray} \tag{ 7 }$$

otherwise, where s_J is the instantaneous s of Jupiter. These equations provide a good approximation for the gas profile given in Walsh et al. (2011) interior to the orbit of Jupiter, the region with which we are concerned. We adopt a temperature profile $T=280\;{s}^{-1/2}$ K as in Brasser et al. (2007), and thus the pressure gradient is given by $\eta =2.13\times {10}^{-3}\;{s}^{1/2}.$ The gas velocity is then obtained from

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{{\rm{g}}}={{\boldsymbol{v}}}_{{\rm{K}}}\;\sqrt{1-2\eta },\end{eqnarray} \tag{ 8 }$$

where ${{\boldsymbol{v}}}_{{\rm{K}}}$ is the Keplerian velocity at a distance s from the star.

In our Grand Tack simulations the gas density begins to decay exponentially as Jupiter migrates outward through the disk, whereas in the calm disk simulations the gas density either decays over the same time period or remains constant throughout.

Figure 2 shows the semimajor axis versus eccentricity at five snapshots in time of a system of planetesimals evolving in the Grand Tack scenario. This broadly reflects the results of the Grand Tack as seen in other works (e.g., Walsh et al. 2011; O'Brien et al. 2014), namely, Jupiter truncates the inner disk and scatters some material outward, leaving some higher-eccentricity planetesimals between 2 and 3 AU to populate the asteroid belt. Note that the banded structure seen beyond ∼1.3 AU in the third panel of Figure 2 and the grouping of excited planetesimals seen at ∼1.9 AU in the bottom two panels are caused by orbital resonances with Jupiter during its (outward) migration and would likely be disrupted by the presence of Saturn.

Comparing the evolution under the Grand Tack model (Figure 2) to the evolution of a calm disk environment (Figure 3), it can be seen that planetesimal eccentricities become much higher in the Grand Tack scenario and larger embryos grow more rapidly. Examining the colors of the planetesimals and embryos, which indicate the initial location of material, it is clear that some blue or green material originating beyond 1.5 AU has been pushed into the inner disk by Jupiter's migration, and that the disk is more mixed compared to the calm disk scenario, with only orange, yellow, or green objects present in the inner disk at the end of the Grand Tack simulation, compared to the calm disk in which the final colors of objects approximately match the initial colors at those locations. It is worth noting that this mixing may be reduced compared to what would be seen with the nominal migration rate as the efficiency of resonant transport may have been changed by the f² increase in the migration rate.

The results of our simulations for calm disks (simulations 1–19) are broadly similar to previous investigations (e.g., Kokubo & Ida 2002; Leinhardt & Richardson 2005), with the same general differences noted in Paper I, namely, that the larger number of particles (smaller starting sizes) and the inclusion of a sophisticated collisional model with fragmentation increase the timescales for runaway and oligarchic growth, as would be expected. The main obvious difference, before we consider the planetesimal compositions, between the simulations presented here and those presented in Paper I is the larger number of resolved fragments. This is due to the improved implementation of fragmentation described by Leinhardt et al. (2015).

The mass distributions at the end of simulations 21 and 5 are shown in Figures 1 and 4. Again this shows that some embryos grow larger in the same time period under the Grand Tack scenario. Here we are seeing the effect of shepherding significant amounts of mass from beyond 1.5 AU into the inner region of the disk; that the Grand Tack simulations result in embryos that are already one Earth mass or more in only ∼20 Myr of evolution is the result of this higher mass in the inner disk. The distinct lines of dark blue colored planetesimals seen in the left panel of Figure 4 show objects that underwent very little evolution in the outer part of the modeled disk (beyond ∼2.5 AU) before being scattered onto highly eccentric orbits by Jupiter. They have subsequently undergone very little or no interaction with other planetesimals, but could be responsible for late delivery of volatile elements (e.g., H, C).

The right panel of Figure 4 shows a very similar result to previous work investigating calm protoplanetary disks (e.g., Leinhardt et al. 2015; Paper I); embryos grow and detach from the mass distribution of the planetesimal population. As noted by Leinhardt et al. (2015), the embryos clean out the planetesimal population, starting in the inner disk where the evolution is faster, leading to an "embryo front" that propagates outward through the disk, and there is no period where the total mass is split uniformly between planetesimals and embryos across the entire disk. This "inside-out" growth that results from the more primitive and realistic initial condition is in contrast to many recent N-body models (e.g., Chambers 2013; Jacobson & Morbidelli 2014). This is also true in the Grand Tack scenario, but with the added complication of the influx of material shepherded by Jupiter.

Figure 5 shows the fraction of each type of collision that occurred as a function of simulation time for the same Grand Tack simulation, 21, and calm disk simulation, 5 (both high-resolution simulations that ignore gas drag; the effects of gas drag will be discussed below). It is immediately apparent that there are more erosive collisions (yellow) and supercatastrophic disruptions (red) under the Grand Tack scenario compared to the calm disk case. The left panel of Figure 5 shows a spike in the number of collisions and fractions of supercatastrophic disruptions and erosive collisions during Jupiter's inward migration, as well as a higher fraction of these disruptive collision types than in the calm disk throughout the rest of the simulation.

Zoom In Zoom Out Reset image size

The collision history for the calm disk simulation shown in the right panel of Figure 5 is very similar to that shown in Paper I. As was seen in both Paper I and Leinhardt et al. (2015), the fraction of disruptive collisions (yellow and red) increases with time as embryos grow, increasing the scattering of planetesimals. Note that the number of collisions decreases strongly with time in both scenarios (white lines), as the number of resolved particles falls. Again the fraction of accretional collisions is approximately constant with time.

Unsurprisingly, the dynamical excitation of the Grand Tack leads to a much more violent environment. The significant number of planetesimals with very large black cores visible in Figure 2 immediately suggests that migration of Jupiter has had a major effect on the planetesimal core fractions. The influence of the migration on core fraction of the growing embryos is more clearly evident in Figure 6. In this case with no gas and the standard fast, early migration (simulation 21, left panel of Figure 6), there is a large variation in core fraction following Jupiter's migration (blue triangles), and some evolution toward lower core fractions is noticeable as the simulation proceeds (black circles). Even though the largest core fraction of an embryo falls during this later evolution, the maximum is still higher than is ever reached in the calm disk, and the range in embryo core fraction is substantially larger (0.15–0.25 compared to 0.20–0.23 for the calm disk).

Zoom In Zoom Out Reset image size

Figure 7 shows the final core fractions of all resolved objects against their final semimajor axis, with the size of the points representing the radii of the planetesimals/embryos. In both simulations there is a small scatter in final core fractions of the embryos around the initial value of 0.22 (as seen in Figure 6), with the largest bodies in the left panel having core fractions very close to the initial value. The smallest resolved planetesimals, however, show a huge variation in core fraction, with some planetesimals composed entirely of iron, and some entirely of silicate mantle material (again in the left panel there are many largely unevolved blue planetesimals with small semimajor axis on high-eccentricity orbits that still have approximately their initial core fraction).

3.2.1. The Mantle Stripping Models

3.2.2. Unresolved Debris

Figure 8 shows the mass of unresolved debris present in the inner disk (0.5–1.5 AU) as the simulations progress. There is a clear spike in the left panel of Figure 8 corresponding to Jupiter's migration, as well as a number of smaller spikes in both simulations caused by collisions.

Zoom In Zoom Out Reset image size

In both simulations shown in Figure 8, the initial core fraction of all planetesimals and the initial unresolved material is 0.22; the average core mass fraction of the debris after 100 kyr is 0.14 and 0.10 in the Grand Tack and calm disk simulations, respectively, and the average mass in unresolved debris after 100 kyr is 4 × 10⁻⁴ M${}_{\oplus }$ and 2 × 10⁻⁴ M${}_{\oplus }.$ The unresolved material is enriched in silicate mantle material, which is as we would expect since it is mantle that is lost first in fragmenting collisions. As was shown by Leinhardt et al. (2015), the total mass in the unresolved component is always a small fraction of the total system mass and never builds to a significant mass.

Figures 9 and 10 show the collision histories and embryo core fractions resulting from slower migration (simulation 23) and late migration (simulation 24).

Zoom In Zoom Out Reset image size

The left panel of Figure 9 shows that the slower migration leads to a larger spike in the fraction of disruptive collisions (yellow and red) and a decrease in the fraction of "perfect" hit-and-run collisions (light blue) compared to the standard Grand Tack scenario (left panel of Figure 5). The hit-and-run collisions in which the projectile remains intact (light blue) cannot affect the composition of either body in our model, whereas the partial accretion collisions (dark blue) that replace them in this slower migration scenario can. The slower migration leads to a noticeable increase in the maximum final core fraction (0.26–0.28 for slow migration compared to 0.24–0.25 for the fast scenario), as shown in the left panel of Figure 10 compared to the left panel of Figure 6.

The late migration (shown in the right panels of Figures 9 and 10) causes an increase in the fraction of disruptive collisions similar to the standard fast, early migration simulations, but the embryo core fractions in this much more evolved planetesimal system do not undergo the same changes. The embryos have grown large enough during their longer period of evolution that the migration has little effect on their composition. The variation in the final core fractions is significantly smaller than in the early Grand Tack simulations, and there is no large variation after Jupiter's migration (green diamonds) to be averaged out by the end of the simulation (black circles).

3.4.1. Calm Disk

We have tested the influence of gas drag from an MMSN gas disk on the evolution of a calm planetesimal disk under two scenarios. In the first, the gas has a density that is constant in time throughout the duration of the simulations; in the second, the gas begins to dissipate after 2.1 Myr (58 kyr simulation time)—approximately 3 Myr after the birth of the solar nebula—and the surface density then decays exponentially on a timescale of 100 kyr (2.8 kyr simulation time). These timescales were chosen to match those we use in the Grand Tack simulations (see Section 3.4.2), which match those from Walsh et al. (2011).

The gas lifetime in the first scenario is probably unrealistically long compared to observed gas disk lifetimes (∼6 Myr; Haisch et al. 2001), but it provides a useful insight into the effects of gas drag on disk and planetary embryo evolution.

This scenario produces smaller oligarchs compared to the case without gas. Such an effect of gas drag was discussed by Tanaka & Ida (1997); it slows the growth of embryos as planetesimals become trapped outside the embryo feeding zones. The left panel of Figure 11 shows that the inner region of the disk has been cleared to a greater degree than in simulations without gas drag (right panel of Figure 7). This reflects the more rapid evolution of the system due to gas drag, as seen by Wetherill & Stewart (1993) and Kokubo & Ida (2000).

Download figure:

Video Standard image High-resolution image

Embryos grow more quickly in the runaway growth phase owing to the eccentricity damping caused by the gas; however, during the oligarchic growth phase, the embryos do not grow to be as massive as they are unable to sweep up the planetesimals. We are thus left at the end of our simulations with a system that has features that are indicative of both more and less evolved systems than the case without gas.

Figure 12 shows the significant change in the fraction of collision types when gas is included; the left panel shows that gas drag results in significantly fewer erosive and disruptive collisions (yellow and red) and a greater fraction of perfect merging collisions (black) compared to the case without gas drag (Figure 5, right panel). Gas drag damps the eccentricities and inclinations of the planetesimals, resulting in lower average collision velocities. The lack of these erosive and disruptive collisions in the constant gas scenario results in very few collisions that can significantly affect the core–mantle balance of the growing embryos. The final distribution of core fractions (left panel of Figure 11) therefore shows very little change from the initial value and much less variation than the simulations without gas drag (Figure 7, right panel).

Zoom In Zoom Out Reset image size

The second scenario, in which the gas dissipates early in the simulation, unsurprisingly has similarities to both the constant gas scenario and the "no gas drag" scenario. The gas drag again facilitates the rapid growth of small embryos; however, once the gas density begins to decay, the embryo eccentricities can be excited, allowing them to collide, producing a smaller number of larger embryos that go on to continue clearing the planetesimal population.

The collision histogram (right panel of Figure 12) looks similar to the no gas drag case (Figure 5, right panel) once the gas begins to dissipate, but in the first 58 kyr (2.1 Myr effective time)—during which the majority of the collisions occur—there are fewer erosive events and more perfect merging collisions. This likely leads to the smaller variation in final core fraction shown in Figure 11 (right panel) for this dissipating gas scenario compared to the case without gas drag (Figure 7, right panel).

3.4.2. The Grand Tack Scenario

Including gas in the standard fast, early migration scenario leads to embryo core fractions that evolve in a similar way to the simulations with slower migration, with the larger variation and enhanced final core fractions compared to the case without gas drag (see Figure 13). While the collision histograms look most similar to the fast, early migration without gas, the inclusion of gas drag has prevented the enhanced core fractions from being reduced by the end of the simulation. The results are similar for both an MMSN gas density profile (simulation 25) and the hydrodynamic derived gas density profile (simulation 26) used by Walsh et al. (2011).

Zoom In Zoom Out Reset image size

Including gas drag in the Grand Tack scenario increases the maximum core fraction compared to the scenario without gas drag; this is the opposite effect to that in the calm disk scenario, in which the damping of eccentricities caused by the drag force leads to a smaller variation in core fraction. The excitation caused by Jupiter's migration has a substantially larger effect than the damping caused by the drag force in an MMSN gas disk.

The "origin histogram" (see Figure 14) of each particle in our simulations tracks the initial locations of the mass of which that particle is composed, so we can easily examine our final embryos and determine from which regions of the disk they have accreted planetesimals. The first panel of Figure 14 shows the composition in terms of radial location of embryos produced in calm disk simulations. This shows a very similar result to that found by Leinhardt et al. (2015), that there is some localized mixing, with most of the mass accreted by each embryo originating from the regions 0.1 AU either side of its position at the end of the simulation.

Zoom In Zoom Out Reset image size

Figure 14 (lower panel) shows the mixing in the same way for the Grand Tack model. It is clear that Jupiter's migration in this model causes a much greater degree of mixing, with almost all embryos containing some significant fraction of material originating beyond 2 AU, and no pair of adjacent 0.1 AU regions dominating the mass of any embryo that remains in the terrestrial planet region. Figure 14 also suggests a much more uniform radial composition due to this increased mixing than is seen in any of the calm disk simulations. Again we note that the accelerated migration may have reduced the efficiency of transport and thus caused reduced mixing.

In this work we have assumed a uniform initial core fraction across the planetesimal disk; this allows for a clearer investigation of the effects of collisional evolution, but in a real system we would expect some variation in oxidation state with heliocentric distance (e.g., Righter et al. 2006). This would be expected to have some effect on the final core fractions of the terrestrial planets, especially in the calm scenario, which shows little mixing. The much more mixed embryos produced in the Grand Tack simulations may not be affected significantly by this simplification; however, the effect on the resulting system of planets will depend on the nature of the initial variation.

In the Grand Tack model, all embryos acquire similar fractions of "blue" material (3%–8% of final mass from beyond 2.5 AU; see Figure 14), which is often assumed to contain more water than planetesimals formed closer to the Sun (e.g., Raymond et al. 2009). If planetesimals originating beyond 2.5 AU are 10% water by mass (as is commonly assumed in similar calculations; e.g., O'Brien et al. 2006), the embryos at the end of these simulations would be 0.4%–0.8% water, and the larger embryos would already have acquired enough to match the Earth's water content. Here we have used only planetesimals originating interior to the orbit of Jupiter; O'Brien et al. (2014) similarly estimate that any Earth-mass planet forming under the Grand Tack would acquire the required amount of water from planetesimals originating beyond the orbit of Jupiter (which we have neglected). However, it is uncertain what fraction of the volatiles from planetesimals accreted during this period would be retained by the final planet, and water from both populations of planetesimals may be required if a substantial fraction is lost owing to collisions.

Bottke et al. (2006) proposed that the iron meteorites, which originate today from the inner asteroid belt, are remnants of planetesimals formed in the terrestrial planet region rather than fragments of bodies formed in the main belt. One outstanding question from Bottke et al.'s work is whether differentiated planetesimals from the terrestrial planet region could actually produce these iron-rich bodies via disruption. The results of this work suggest that many iron-rich, and complimentary silicate-rich, fragments are indeed produced by the collisional evolution of the inner disk, as seen in Figure 7. Though we see very few of these remnant planetesimals scattered into the asteroid belt region, it is important to note that the initial distribution of planetesimals in our calm disk simulations does not extend beyond 1.5 AU, and that such material may be important in scattering fragments into the main belt region (and keeping them there).

It has recently been realized that many white dwarfs display extrasolar planetary debris in their spectra (e.g., Jura & Young 2014). This material shows a diversity in the iron contents of accreted planetesimal and small embryo sized bodies (e.g., Farihi et al. 2011; Gänsicke et al. 2012; Raddi et al. 2015) that could represent leftover collisionally processed planetesimals similar to those produced in our simulations (e.g., Figure 7). The presence of differentiated rocky debris in white dwarf atmospheres appears to be a common phenomenon, which may suggest that many extrasolar planetary systems undergo similar collisional evolution to that shown in this work.

In order to estimate the bulk Fe/Mg ratios of the final planetary embryos, we use a series of starting bulk compositions for the planetesimals that span the range represented by chondritic meteorites. Following the technique used in Paper I, for each bulk planetesimal composition, a metallic core is formed of the proscribed mass (either 0.22 or 0.35, see Section 2.4, with 10% of the core as light elements in the former or 5% light elements and 5% Si in the latter) and a residual mantle composition calculated. These compositions are used together with the modeled mass fractions of core and mantle in the final embryos (see Figure 13) to determine their bulk Fe/Mg. Finally, we apply the additional constraint that the embryo must have enough Fe to match the Earth's measured core mass fraction (the FeO/Fe ratio is adjusted in order to match the Earth's core fraction of 32%, ensuring that the FeO content of the mantle is zero or greater).

Figure 15 illustrates the differences between starting, chondritic compositions (dashed lines) and bulk Earth composition (dark gray bar), indicating the need for some processing of chondritic material in order to produce the Earth. The ranges in Fe/Mg that result from a selection of our simulations are also shown in Figure 15, and it is evident that collisional processing can produce embryos matching Earth's Fe/Mg ratio. Earth-like Fe/Mg can be produced in extreme cases starting with planetesimal compositions close to bulk Earth in calm scenarios and more readily in the Grand Tack models. The collisional stripping of mantle caused by the Grand Tack can be sufficient to explain the non-chondritic bulk Fe/Mg ratio of the Earth starting with a wider range of planetesimal compositions.

Zoom In Zoom Out Reset image size

It is, however, important to note that the accretion of terrestrial planets is not complete at the end point of these simulations, and that it is not necessarily the case that the embryos with enhanced core fractions represent the Earth. Collisions during the giant impact phase may cause any compositional differences observed at the end of oligarchic growth to be averaged out again, or they may serve to enhance this variation (for more detailed discussions of compositional changes during the giant impact phase see Stewart & Leinhardt 2012; Paper I).

The compositional variation that we see in these simulations also has strong implications for exoplanetary systems. The bulk composition of planets could be significantly different from the expected composition based on the stellar photosphere. The CI chondrites are a good match for our Sun, but as Figure 15 shows, the bulk composition of planetary embryos formed from CI chondrites may not be the same. Given the interest in the compositions of exoplanets and their atmospheres, it is important to consider the possible deviation from the composition implied by their host stars, especially if they may have undergone dynamical processes similar to the Grand Tack.