ANALYSIS OF TERRESTRIAL PLANET FORMATION BY THE GRAND TACK MODEL: SYSTEM ARCHITECTURE AND TACK LOCATION

R. Brasser; S. Matsumura; S. Ida; S. J. Mojzsis; S. C. Werner

doi:10.3847/0004-637X/821/2/75

1. INTRODUCTION

A successful physical model for the formation of the terrestrial planets is a long-standing problem (Tsiganis 2015). The first physically plausible idea came from Safronov (1969), who suggested that the earliest stage of the accumulation of dust into larger bodies was caused by gravitational instability in a thin dust layer. Safronov (1969) showed that the relative velocities between bodies are of the order of their escape velocity, and so the largest body's gravitational cross section is limited by the geometrical cross section, limiting growth. These findings were later used by Wetherill (1980), who showed that the terrestrial planets coagulated from planetesimals, and that the formation of the these planets was linked to the evolution of the asteroid belt. Wetherill & Stewart (1989) elaborated that in a disk of planetesimals, some would undergo runaway growth and form a sequence of planetary embryos. These embryos would then further collide to form the terrestrial planets.

These ideas were first rigorously tested by Kokubo & Ida (1996), who performed numerical simulations of a self-gravitating disk of planetesimals. They discovered that some objects in the disk underwent runaway growth, as was predicted, which resulted in a mixed population of protoplanets and planetesimals (Kokubo & Ida 1998). The protoplanets underwent so-called oligarchic growth: all would be roughly equally spaced and of similar mass as each vied for supremacy in accreting the last remaining planetesimals. The protoplanets (also dubbed "planetary embryos") would subsequently collide to form the terrestrial planets (Chambers 2001).

Early simulations of terrestrial planet formation yielded estimates for a growth timescale of several tens of millions of years and overall results which showed that the final terrestrial system would be assembled by 100 Myr (Chambers 2001). Most of these early simulated systems, however, were found to suffer from excess eccentricity and inclination for the final planets, but the inclusion of a large number of planetesimals to exert dynamical friction alleviated this concern (O'Brien et al. 2006). A further chronic and fundamental shortcoming of earlier simulations was that the output systematically yielded a far too massive Mars analog. This predicament led Raymond et al. (2009) to investigate how the mass of Mars might depend on the orbital configuration of the giant planets. It was found that only the current spacing of the gas giants led to the model being capable of producing a Mars analog much less massive than Earth, but under the special condition that the eccentricities of the gas giants were higher than their current values. As such, Raymond et al. (2009) highlighted the unrealistic nature of the initial conditions required to explain Mars' low mass, and left the problem as a lingering problem to be solved later.

A potential solution presented itself in the work of Hansen (2009), who studied terrestrial planet formation with planetary embryos situated in a narrow annulus between 0.7 and 1 au from the Sun. These initial conditions nicely reproduced the mass-semimajor axis relationship we have today, with two relatively large terrestrial-type planets book-ended by two much less massive planets. The main drawback of that study was that no mechanism was presented to gravitationally truncate the outer edge of the solid disk near 1 au. The same could be said for the inner edge, so that no mechanism existed to create such a high-density, narrow annulus with which to explain the terrestrial worlds.

This quandary led Walsh et al. (2011) to propose the so-called Grand Tack scenario, wherein the early, gas-driven coupled migration of Jupiter and Saturn sculpts the planetesimal disk and truncates it near 1 au. The Grand Tack at least partially explains the formation of a high-density region in the inner disk, although it cannot explain the existence of an inner cavity inside roughly 0.7 au. The inclusion of this scenario leads to a broad outline of how the early solar system evolved. First, Jupiter is assumed to form before Saturn, clearing the gas in an annulus with a width comparable to its Hill radius, and undergoing inward Type 2 migration (Lin & Papaloizou 1986). The inward migration of Jupiter shepherds material toward the inner portion of the disk while also scattering other material outward, creating an enhanced density region for terrestrial planet formation and mixing planetesimals from the innermost and outer portions of the protoplanetary disk. Once Saturn grows to about half of its current mass (Masset & Casoli 2009), it is assumed to partially clear the disk in its vicinity, migrate rapidly at first to catch up with Jupiter, and subsequently become trapped in a mean-motion resonance near Jupiter, presumably the 2:3 resonance (Masset & Snellgrove 2001; Pierens & Raymond 2011) but it may also have been the 2:1 (Pierens et al. 2014). In this process, these two giant planets clear the disk together. The torque from the interaction with the disk is stronger for a shorter separation between the planet and the disk edge. Since Jupiter and Saturn are interacting and Jupiter is more massive, it is reasonable to think that Saturn is pushed outward by Jupiter's perturbation, and the separation from the disk edge is smaller for Saturn than for Jupiter. Thus, the torque on Saturn can be larger despite the planet's lower mass. The interaction with Jupiter prevents Saturn from creating a cleared annulus in the disk and allows gas from the outer disk to flow past Saturn and into the inner disk. If the gap-crossing disk gas flow is large enough, then the Jupiter-Saturn pair can migrate outward (Masset & Snellgrove 2001; Pierens & Nelson 2008). Consequently, the planets reverse their migration: they tack as a sail boat would change its direction by steering into and through the wind. Once the giant planets have completed this early migration phase, have left the inner solar system, and settled in the vicinity of their present positions, terrestrial planet formation could proceed as before, but only (as advocated by Hansen 2009) from material in a narrow circumsolar annulus. In this manner, Walsh et al. (2011) successfully reproduced the mass-semimajor axis distribution of the inner planets if the reversal of Jupiter occurred at 1.5 au because they truncated the inner edge of their planetesimal disk at 0.5 au. A successful feature of their model is that it also accounts for the apparent compositional differences across the asteroid belt (DeMeo & Carry 2014).

It is worth noting, however, that this stepwise reconstruction of the early evolution of the planetary system has some pitfalls of varying severity. For example, the particular configuration and outward migration of the giant planets favored by the Grand Tack is only supported for a narrow set of initial conditions (D'Angelo & Marzari 2012) and is not universal (Zhang & Zhou 2010). There may also be other pathways to produce such a high-density region through a deficit of material near Mars (Izidoro et al. 2014), although this idea was recently undermined in a follow-up study (Izidoro et al. 2015). Finally, up to this point, the Grand Tack fails to reproduce the current mass and location of Mercury, most likely because dynamical models always truncate the disk near 0.5 au or beyond. Clearly, further study is needed of both the gas-driven evolution of the giant planets and the subsequent formation and evolution of the terrestrial planets in order to explain what we see in our own solar system.

With this in mind, we sought to scrutinize the Grand Tack model and its consequences over most of the age of the solar system by running a large number of Grand Tack simulations with a range of initial conditions and varying tack locations. Our work also includes several dynamical effects that have hitherto been ignored. We report on the various methods that we employ to quantify whether or not Grand Tack successfully reproduces the observed dynamical features of the inner solar system and whether one set of initial conditions and tack locations is more favorable than another.

This report is organized as follows. In Section 2, we introduce several additions to the original Grand Tack simulations of Walsh et al. (2011) and justify our choice of disk model and the inclusion of type 1 migration. Section 3 describes our initial conditions, while Section 4 describes our numerical methods. Section 5 explains our criteria for a set of simulations to successfully reproduce the currently observed dynamical properties of the inner solar system. This is followed by Sections 6 and 7 where we describe the results of our numerical simulations. Section 8 is reserved for a discussion, and we draw our conclusions in Section 9.

2. DEVIATIONS FROM THE ORIGINAL GRAND TACK MODEL

In addition to the simple reproduction of the Grand Tack scenario, for this study we also chose to employ a substantially different model for the protoplanetary disk than that used by Walsh et al. (2011). An explanation for this choice is provided below. We have also included the effect of type 1 embryo migration.

2.1. Protoplanetary Disk

Walsh et al. (2011) employed the protoplanetary disk model of Morbidelli & Crida (2007), which in turn was based on the work of Guillot & Hueso (2006). The surface density of their disk profile is of the form ${\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{0}\mathrm{exp}(-{r}^{2}/{R}^{2})$ , where $R\sim 200\;\;{\rm{au}}$ is a scaling constant. This Gaussian profile of the surface density is markedly different from the often employed power-law slopes found elsewhere in the literature, e.g., Hartmann et al. (1998). The scaling constant at 1 au is ${{\rm{\Sigma }}}_{0}=100$ g cm⁻², which is much lower than the usual value of 1700–2400 g cm⁻² (Hayashi 1981). Since the disk model of Morbidelli & Crida (2007) is not widely used and relies on a constant viscosity, ν, rather than a constant α viscosity, we decided to make use of the disk model of Bitsch et al. (2015), which is based on the study by Hartmann et al. (1998). Bitsch et al. (2015) give fitting formulae to compute the disk's surface density, temperature, and scale height as a function of both heliocentric distance and time. Initially, the disk's gas surface density at 1 au is 2272 g cm⁻² and the temperature is 576 K in the midplane, and so the scale height at 1 au is about 0.057 au and the metallicity is equal to the solar value. This is higher than most traditional models have assumed; the higher temperature is caused by viscous heating. These fitting formulae are valid, however, as long as any embedded planet is not massive enough to significantly alter the disk structure, such as by opening an annulus, and as long as the disk remains mostly unperturbed. This latter requirement may not be entirely true because of the proximity of Jupiter, whose presence imposes a change in the disk's surface density and temperature. For a radiative disk, this change in the temperature profile will result in a change in the disk surface density, but for a viscous disk, which is the region we work in, this effect is less severe. Besides, a much lower disk temperature, caused by the influence of Jupiter, would mean that the ice line is close to 1 au very early on, which is inconsistent with solar system formation (Matsumura et al. 2016). Here, we adopt the disk of Bitsch et al. (2015), set theα viscosity equal to 0.005, use a solar metallicity, and adopt a molecular weight of the gas of 2.3 amu (Bitsch et al. 2015). The gas surface density scales as ${\rm{\Sigma }}(r)\propto {r}^{-\alpha }$ and the temperature profile is $T(r)\propto {r}^{-\beta }$ , where $\alpha =1/2$ and $\beta =6/7$ . These power-law relations are accurate out to ∼5 au, which is the region we are interested in, and so we adopted these profiles throughout the disk. In Figure 1, we plot the evolution of the surface density (top) and the temperature bottom as a function of time and distance to the Sun. There is a very rapid decay for the first 1 Myr and a slower decay after that. After 5 Myr, we photo-evaporate the disk away over the next 100 kyr.

**Figure 1.** Contour plots of $\mathrm{log}({\rm{\Sigma }}/1\;{\rm{g}}\;{\mathrm{cm}}^{-2})$ (top) and $\mathrm{log}(T/1\;{\rm{K}})$ (bottom) as a function of distance to the Sun (horizontal axis) and age in Myr (vertical axis).
Download figure:
Standard image High-resolution image

**Figure 1.** Contour plots of $\mathrm{log}({\rm{\Sigma }}/1\;{\rm{g}}\;{\mathrm{cm}}^{-2})$ (top) and $\mathrm{log}(T/1\;{\rm{K}})$ (bottom) as a function of distance to the Sun (horizontal axis) and age in Myr (vertical axis).
Download figure:
Standard image High-resolution image

2.2. Embryo Migration

Although Walsh et al. (2011) decided not to include the effect of type 1 migration (Tanaka et al. 2002) on the planetary embryos, we decided to take it into account to determine whether (or not) it drastically affects the outcome of the simulations. For the migration prescription, we follow Cresswell & Nelson (2008), which is partially based on the work of Tanaka et al. (2002). We chose not to employ the non-isothermal approach of Paardekooper et al. (2011) because generally none of the terrestrial planets are massive enough to begin outward migration, except during the very late stages when the disk surface density is low. Thus, for simplicity, we adhere to the isothermal case. The specific decelerations experienced by the planetary embryos, due to the disk, in eccentricity, inclination, and semimajor axis are given by

$\begin{eqnarray}&&{{\boldsymbol{a}}}_{e}=-2\displaystyle \frac{({\boldsymbol{v}}\cdot {\boldsymbol{r}})}{{r}^{2}{t}_{e}}{\boldsymbol{r}},\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&{{\boldsymbol{a}}}_{i}=-\displaystyle \frac{{v}_{z}}{{t}_{i}}{\boldsymbol{k}},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&{{\boldsymbol{a}}}_{m}=-\displaystyle \frac{{\boldsymbol{v}}}{{t}_{m}},\end{eqnarray} \tag{ 3 }$

where ${\boldsymbol{r}}$ and ${\boldsymbol{v}}$ are the position and velocity vectors of the embryo, ${\boldsymbol{k}}$ is the unit vector in the z-direction, v_z is the z-component of the velocity, and t_e, t_i, and t_m are the timescales to damp the eccentricity, inclination, and semimajor axis. The latter quantities depend in a complicated manner on the semimajor axis, eccentricity, inclination, surface density, and temperature of the gas, and the mass of the embryo. We refer the interested reader to Cresswell & Nelson (2008) and Tanaka et al. (2002) for details and for the conversion to the cartesian frame. All of the above timescales are a function of the wave time, given by (Tanaka & Ward 2004)

$\begin{eqnarray}{t}_{{\rm{wav}}} & = & \left(\displaystyle \frac{{M}_{\odot }}{{m}_{{\rm{emb}}}}\right)\left(\displaystyle \frac{{M}_{\odot }}{{\rm{\Sigma }}{r}^{2}}\right){\left(\displaystyle \frac{c}{r{{\rm{\Omega }}}_{{\rm{K}}}}\right)}^{4}{{\rm{\Omega }}}_{{\rm{K}}}^{-1},\\ & = & \left(\displaystyle \frac{{M}_{\odot }}{{m}_{{\rm{emb}}}}\right)\left(\displaystyle \frac{{M}_{\odot }}{{\rm{\Sigma }}{r}^{2}}\right){\left(\displaystyle \frac{H}{r}\right)}^{4}\left(\displaystyle \frac{r}{{v}_{{\rm{K}}}}\right)\end{eqnarray} \tag{ , 4 }$

where we used $H/r=c/{v}_{{\rm{K}}}$ , ${{\rm{\Omega }}}_{{\rm{K}}}$ is the orbital frequency, ${c}^{2}=\gamma {k}_{{\rm{B}}}T/\mu {m}_{{\rm{p}}}$ is the sound speed, ${v}_{{\rm{K}}}$ is the orbital speed, H is the scale height, γ is the ratio of specific heats (taken as 7/5), ${k}_{{\rm{B}}}$ is the Boltzmann constant, ${m}_{{\rm{p}}}$ is the mass of the proton, and μ is the molecular weight of the gas, assumed to be 2.3 amu. The value of ${t}_{{\rm{wav}}}$ depends sensitively on the slopes of the surface density and temperature profiles, but the product ${m}_{{\rm{emb}}}{t}_{{\rm{wav}}}\propto {a}^{3/2-2\beta +\alpha }$ is much less sensitive and is what determines the migration rate of the embryos. We generally have ${m}_{{\rm{emb}}}{t}_{{\rm{wav}}}\sim {a}^{2/7}$ to a² for our disk profile and the minimum disk profile of Hayashi (1981; which has $\alpha =3/2$ and $\beta =1/2$ ).

Walsh et al. (2011) included the first two deceleration contributions in Equations (1)–(3) because these are tidal effects cause by the disk that damp the eccentricity and inclination, but they omitted the third (Equation (3)), which is in part responsible for the inward migration of the planetary embryos. For the disk that we have chosen, near 1 au we compute ${t}_{m}\sim 2(0.1\;{M}_{\oplus }/{m}_{{\rm{emb}}})(2000\;{\rm{g}}\;{{\rm{cm}}}^{-2}/{\rm{\Sigma }})$ Myr, so that any inward migration the embryos experience will likely be severely restricted because the migration timescale is longer than the disk lifetime and increases as the disk surface density decreases. However, ${t}_{e}\;\sim 21$ kyr for a Mars-sized body near 1 au, and so the damping effect is a lot stronger than the migration. For the disk employed by Walsh et al. (2011), the migration and damping timescales are both at least an order of magnitude longer, so that their effects of type 1 migration are very weak.

To compare our migration times with earlier results, we note that for an Earth-sized body at 1 au, our nominal disk parameters yield ${t}_{{\rm{wav}}}\sim 1700\;{\rm{years}}$ and ${t}_{m}\sim 180$ kyr, which are much longer than the values of Tanaka et al. (2002) and Tanaka & Ward (2004) because our disk is initially hotter. The wave time scales as ${(H/r)}^{4}$ , so that a factor of 1.5 in H/r will result in a factor of 5 in the migration timescale; meanwhile, ${t}_{{\rm{m}}}\propto {(H/r)}^{2}$ , and so the effect of the disk temperature on it is weaker. As the disk evolves, the temperature decreases, speeding up type 1 migration, but the surface density also decreases, slowing it down. In general, the effect of type 1 migration weakens with time.

Thus, the effect of type 1 migration should be relatively weak in our disk model compared to the traditional results of Tanaka et al. (2002) and Tanaka & Ward (2004). Apart from the migrating force, the eccentricity and inclination damping forces from the disk will also cause some inward migration of the embryos due to angular momentum loss, and we expect them to migrate inward on the order of 0.1 au over the lifetime of the disk. In Figure 2, we plot a contour map of the normalized torque on the planetary embryos as a function of their distance to the Sun and their mass at the zero-time age of the disk. The normalized torque is ${{\rm{\Gamma }}}_{{\rm{n}}}={{\rm{\Gamma }}}_{{\rm{tot}}}/{{\rm{\Gamma }}}_{0}$ , where ${{\rm{\Gamma }}}_{0}={({m}_{{\rm{emb}}}/{M}_{\odot })}^{2}{(r/H)}^{2}{(v/r)}^{2}{\rm{\Sigma }}$ . The migration rate is then $\dot{r}=-2r{{\rm{\Gamma }}}_{{\rm{tot}}}/L$ , where L is the orbital angular momentum (Tanaka et al. 2002). In all of the disk models, the definition of the torque is always inward.

**Figure 2.** Contour plots of the normalized torque on a planetary embryo as a function of semimajor axis and mass.
Download figure:
Standard image High-resolution image

3. INITIAL CONDITIONS

For this project, we have run a large sample of numerical simulations of the Grand Tack scenario. These simulations are categorized into four large sets, with further subdivisions therein.

All of our simulations start with a high number of small planetesimals, planetary embryos, and with the gas giants Jupiter and Saturn. We do not include Uranus and Neptune in any of the simulations because they do not have any immediate effects on the formation of the terrestrial planets (Walsh et al. 2011). In all of our simulations, we chose to not take into account the effect of Saturn's mass growth. According to Walsh et al. (2011) and Jacobson & Morbidelli (2014), other effects, such as the radial evolution of Saturn and the gas giant migration timescale, did not substantially change the final terrestrial planet systems. Given the high number of free parameters, we decided to follow Walsh et al. (2011) with cases where Jupiter is assumed to have its current mass and is initially placed on a near-circular orbit at 3.5 au. A fully grown Saturn is placed in the 2:3 resonance with Jupiter at 4.5 au. During the first 0.1 Myr, Jupiter and Saturn migrate from their initial locations (3.5 and 4.5 au) to the tack locations (either 1.5 or 2.0 au for Jupiter, respectively). For the next 5 Myr, Jupiter and Saturn migrate out to ∼5.4 au and ∼7.5 au, which are appropriate initial conditions for late giant planet migration models (Morbidelli et al. 2007). Walsh et al. (2011) demonstrated that the migration speed of the gas giants has almost no influence on the final architecture of the terrestrial system. We therefore used the same linear inward migration of the gas giants with a timescale of 0.1 Myr, followed by outward migration via an exponential prescription with an e-folding time of 0.5 Myr. These timescales are comparable to the typical rate of Type 2 migration (Lin & Papaloizou 1986).

To compare the effects of different formation models on the architecture of the terrestrial planets, we use various initial distributions of embryos and planetesimals.

3.1. Equal Mass Embryo Initial Conditions

For the first two large sets of simulations, we employ the initial conditions of the embryos and planetesimals from Jacobson & Morbidelli (2014), but we use our model for the protoplanetary disk. These simulations were run because we want to directly compare the results of our modified simulations—employing a different protoplanetary disk, including type 1 migration and a realistic mass–radius relationship—with theirs. All of the simulations in this set use equal-mass embryos initially situated between 0.7 and 3 au. The embryos were embedded in a disk of planetesimals. The surface densities of the embryos and planetesimals both scaled with heliocentric distance as ${r}^{-3/2}$ . Following Jacobson & Morbidelli (2014), the total mass ratio of embryos and planetesimals in this inner disk is 1:1, 4:1, or 8:1, with the individual embryo masses being 0.025, 0.05, or 0.08 ${M}_{\oplus }$ . The equal-mass embryo assumption appears to agree with a pebble-accretion scenario of embryo formation (Morbidelli et al. 2015; Levison et al. 2015) rather than the more traditional oligarchic growth scenario (Kokubo & Ida 1998), although which scenario is favored is still a matter of considerable debate. To mimic the coagulation evolution of the solids in the disk, we follow Chambers (2006) and calculate the age of the disk to be 0.1 Myr when embryos have a mass of 0.025 ${M}_{\oplus }$ , 0.5 Myr when the embryos have a mass of 0.05 ${M}_{\oplus }$ , and 1 Myr when the embryos have a mass of 0.08 ${M}_{\oplus }$ .

Following Jacobson & Morbidelli (2014) again, the total mass in solids in the inner disk (embryos and planetesimals) is 4.3 ${M}_{\oplus }$ when the total mass ratio between embryos and planetesimals is 1:1, 5.3 ${M}_{\oplus }$ when the mass ratio is 4:1, and 6.0 ${M}_{\oplus }$ when the mass ratio is 8:1. Jacobson & Morbidelli (2014) further argued that these different initial disk masses were necessary to maintain the post-migration mass in solids between 0.7 and 1 au close to 2 ${M}_{\oplus }$ . Therefore, the surface density in the solids between these different initial conditions increases with increasing total embryo to planetesimal mass ratio. The number of planetary embryos ranged from 29 to 213 depending on their initial mass and total mass ratio. We kept the number of planetesimals at 2000, regardless of their total mass. The initial densities of the planetesimals and embryos was 3 g cm⁻³ (Walsh et al. 2011). The permutations of these initial conditions result in nine individual sets of simulations.

Following Matsumura et al. (2016), for most of the simulations, we also added an outer disk of planetesimals. This disk consists of 500 planetesimals with a total mass of 0.06 ${M}_{\oplus }$ . These planetesimals are, to some degree, considered responsible for volatile delivery to the otherwise dry terrestrial planets (Matsumura et al. 2016). The outer planetesimals are distributed between 5 and 9 au. The eccentricities and inclinations of both embryos and planetesimals are randomly chosen from a uniform distribution between 0 and 0.01 and 0 to 0 fdg 5, respectively. The other angular orbital elements were chosen uniformly at random from 0° to 360°. We ran two sets of these nine permutations, one with a tack at 1.5 au as in Walsh et al. (2011), and one with a tack at 2 au as described in Matsumura et al. (2016).

3.2. Oligarchic Initial Conditions

In addition to simulating the formation of the terrestrial planets from a disk of equal-mass embryos and planetesimals, we also run a second set of simulations where the initial conditions are reminiscent of the traditional oligarchic growth scenario (Kokubo & Ida 1998). To set up our simulations, we used the semi-analytical oligarchic approach of Chambers (2006). In that work, the mass of embryos increase up to their isolation mass as

$\begin{eqnarray}&&{m}_{{\rm{p}}}(t)={m}_{{\rm{iso}}}{\mathrm{tanh}}^{3}\left(\displaystyle \frac{t}{\tau }\right),\end{eqnarray} \tag{ 5 }$

where ${m}_{{\rm{iso}}}=2\pi a{{\rm{\Sigma }}}_{s}b$ is the isolation mass at semimajor axis a for embryos spaced b au apart embedded in a disk with solid surface density ${{\rm{\Sigma }}}_{s}$ . Here, τ is the growth timescale, which is a complex function of the semimajor axis, embryo spacing, solid surface density, and the radii of planetesimals that accrete onto the embryos. The growth timescale is given by (Chambers 2006)

$\begin{eqnarray}&&\tau =\displaystyle \frac{2{e}_{{\rm{Hi,eq}}}^{2}}{A}\end{eqnarray} \tag{ , 6 }$

where

$\begin{eqnarray}{e}_{{\rm{Hi,eq}}}^{2} & = & {2.7}^{2}{\left(\displaystyle \frac{{r}_{{\rm{c}}}\rho }{{{bC}}_{D}a{\rho }_{{\rm{gas}}}}\right)}^{2/5}\\ & = & 17.5{\left(\displaystyle \frac{{r}_{{\rm{c}}}}{10{\rm{km}}}\right)}^{2/5}{\left(\displaystyle \frac{\rho }{2.5{\text{g cm}}^{-3}}\right)}^{2/5}{\left(\displaystyle \frac{b}{10{R}_{H}}\right)}^{-2/5}\\ & & \times {\left(\displaystyle \frac{{\rho }_{{\rm{gas,0}}}}{1.4\times {10}^{-9}{\text{g cm}}^{-3}}\right)}^{-2/5}{\left(\displaystyle \frac{a}{1{\rm{au}}}\right)}^{2/5\alpha +1/5-1/5\beta },\end{eqnarray} \tag{ 7 }$

assuming the drag coefficient C_D = 2 and

$\begin{eqnarray}&&\displaystyle \frac{1}{A}=\displaystyle \frac{{b}^{1/2}P{\rho }^{1/3}{M}_{\odot }^{1/6}}{31.7{{\rm{\Sigma }}}_{s}^{1/2}}\end{eqnarray} \tag{ 8 }$

$\begin{eqnarray} & = & 30.1\;{\rm{kyr}}\;{\left(\displaystyle \frac{b}{10{R}_{H}}\right)}^{1/2}{\left(\displaystyle \frac{\rho }{2.5{\text{g cm}}^{-3}}\right)}^{1/3}\\ & & \times {\left(\displaystyle \frac{{{\rm{\Sigma }}}_{s}}{10{\text{g cm}}^{-2}}\right)}^{-1/2}{\left(\displaystyle \frac{a}{1{\rm{au}}}\right)}^{3/2-1/2\alpha }{M}_{\odot }^{1/6}.\end{eqnarray} \tag{ 9 }$

Here, α is the slope of the solid surface density, β is the slope of the temperature, ${\rho }_{{\rm{gas,0}}}$ is the gas density at 1 au in the midplane, ${r}_{{\rm{c}}}$ is the radius of planetesimals, ρ is their density, and P is the orbital period (Chambers 2006). Here, we used the fact that the gas density profile depends on the scale height and the gas surface density, which yields ${\rho }_{{\rm{g}}}\propto {a}^{-\alpha -3/2+1/2\beta }$ . Combining all of the above gives

$\begin{eqnarray}\tau & = & 527\;{\rm{kyr}}{\left(\displaystyle \frac{b}{10{R}_{H}}\right)}^{1/10}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{s}}{10{\rm{g}}{\mathrm{cm}}^{-2}}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{\rho }{2.5{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{11/15}{\left(\displaystyle \frac{{r}_{{\rm{c}}}}{10{\rm{km}}}\right)}^{2/5}\\ & & \times {\left(\displaystyle \frac{{\rho }_{{\rm{gas,0}}}}{1.4\times {10}^{-9}{\text{g cm}}^{-3}}\right)}^{-2/5}\\ & & \times {\left(\displaystyle \frac{a}{1{\rm{au}}}\right)}^{17/10-1/10\alpha -1/5\beta }{M}_{\odot }^{1/6},\end{eqnarray} \tag{ 10 }$

where ${{\rm{\Sigma }}}_{s}$ is the solid disk surface density at 1 au and we assumed that ${{\rm{\Sigma }}}_{s}$ has the same radial dependence as the gas surface density. This derived value of τ and the subsequent growth of any embryo near 1 au agrees well with Figure 1 in Chambers (2006). Adopting $\alpha =3/2$ and $\beta =6/7$ , the above equation is a reasonably steep function of the semimajor axis at a fixed epoch because $\tau \propto {a}^{193/140}\sim {a}^{1.38}$ . When we use the typically assumed value $\beta =1/2$ , then we have $\tau \propto {a}^{3/2}$ . At the Earth's location, the growth time for nominal parameters is ∼600 kyr, while at Mars' current orbit the growth time is ∼1 Myr. This is somewhat shorter than that advocated by Dauphas & Pourmand (2011), but is within error margins and is easily increased to 1.8 Myr assuming that the accretion was caused by planetesimals of ∼50 km.

We constructed our initial disk of embryos and planetesimals as follows. First, we computed the total mass in solids between 0.7 au and 3 au assuming the surface density in solids to be ${{\rm{\Sigma }}}_{s}=7$ g cm ${}^{-2}{(a/1{\rm{au}})}^{-3/2}$ . This setup is just the minimum-mass solar nebula (Hayashi 1981). Second, following Ogihara & Ida (2009), we subsequently increased the solid density by a factor of three at the ice line, assumed to be static at 2.7 au. Third, based on the results of Kokubo & Ida (1998), we imposed a spacing of 10 mutual Hill radii for the embryos, with the spacing computed assuming that the embryos had their isolation masses. In other words, the semimajor axis of embryo n is ${a}_{n}={a}_{n-1}[1+b{(2{m}_{{\rm{iso}}}/3{M}_{\odot })}^{1/3}]$ , so that the embryo spacing nearly follows a geometric progression. Following Chambers (2006), we assumed a planetesimal size of 10 km in computing the growth timescale of the embryos. Finally, we entered the epoch at which Jupiter was assumed to have fully formed and began migrating. This was either 0.5 Myr, 1 Myr, 2 Myr, or 3 Myr. Most embryos have only reached a fraction of their isolation mass and the remaining mass within their feeding annulus of 10 Hill radii is taken up by planetesimals, with each planetesimal having a mass of 10⁻³ ${M}_{\oplus }$ . The eccentricities of the embryos and planetesimals followed a Rayleigh distribution with a scale parameter equal to ${({m}_{{\rm{p}}}/3{M}_{\odot })}^{1/3}$ . The inclinations also followed a Rayleigh distribution with a scale parameter equal to half that of the eccentricities. The other angles were chosen uniformly at random between 0° and 360°. All of the embryos and planetesimals had an initial density of 3 g cm⁻². The most distant embryo was typically at 2.6 au.

4. METHODOLOGY

The gas giants, planetary embryos, and planetesimals are simulated using the symplectic SyMBA integrator (Duncan et al. 1998) with a time step of 0.02 years for 150 Myr. The end time of the simulations corresponds closely to the end of the purported Late Veneer (see the next section). The migration of the gas giants was mimicked through the fictitious forces described in Walsh et al. (2011). We first employed the same gas profile as in Walsh et al. (2011) with two large dips around the gas giants, which migrated inward with these planets, but we modified the gas density profile during the computation so that it followed the ${\rm{\Sigma }}(r)\propto {r}^{-1/2}$ surface density law of Bitsch et al. (2015) rather than the Gaussian of Walsh et al. (2011). The initial total mass of the disk was approximately 0.05 ${M}_{\odot }$ . The initial disk profile is depicted in Figure 3. The blue line is the disk from Bitsch et al. (2015), whereas the red one is from Walsh et al. (2011). Jupiter is assumed to be at 3.5 au and the disk age is zero. The planetary embryos experienced tidal damping of their eccentricities and inclinations and a negative torque from the gas disk as described in Section 2, and the orbits of the planetesimals evolve due to gas drag using the methods of Brasser et al. (2007). Following Walsh et al. (2011) and Jacobson & Morbidelli (2014), for the purpose of the gas drag, we assumed that each planetesimal had a radius of 50 km. These planetesimals are larger than the 10 km size assumed for oligarchic growth because we wanted to compare our simulations directly with those of Jacobson & Morbidelli (2014). The gas drag routines of Brasser et al. (2007) were modified slightly to allow for a smooth transition between those regimes where the Knudsen number crossed 1. We now have the following scheme:

1.
when ${ \mathcal M }\gt 2.727$ , C_D = 2 for all ${ \mathcal K }$ and ${ \mathcal R };$
2.
when ${ \mathcal R }\gt 1000,{C}_{D}\;$ $=\;0.44+0.2098{{ \mathcal M }}^{2}$ for ${ \mathcal M }\lt 2.727;$ and
3.
when ${ \mathcal R }\lt 1000$

$\begin{eqnarray}&&{C}_{D}=0.2689{{ \mathcal M }}^{2}+\displaystyle \frac{24}{{ \mathcal R }}\left[1-{\left(\displaystyle \frac{{ \mathcal M }}{2.727}\right)}^{2}\right]\displaystyle \frac{1+0.15{{ \mathcal R }}^{0.687}}{\sqrt{{{ \mathcal K }}^{2}+1}}.\end{eqnarray} \tag{ 11 }$

Here ${ \mathcal M }$ , ${ \mathcal K }$ , and ${ \mathcal R }$ are the Mach, Knudsen, and Reynold numbers of the planetesimals, and C_D is the drag coefficient (Brasser et al. 2007). In addition to the drag coefficient routine, we made one further modification to the code.

**Figure 3.** Density profile of the gas disk with Jupiter at 3.5 au and the disk age at 0. The blue line is the disk we employed from Bitsch et al. (2015) and the red line is the disk profile of Walsh et al. (2011).
Download figure:
Standard image High-resolution image

SyMBA treats collisions between bodies as perfect mergers, preserving their density. This works well in most circumstances, but given that the mean density of Earth is 5.5 g cm⁻² and not 3 g cm⁻², we implemented the mass–radius relationship of Seager et al. (2007) to make sure that the final planets have radii comparable to the current terrestrial planets, so that their collisional cross sections are not artificially large. The relation we employed was

$\begin{eqnarray}\mathrm{log}\left(\displaystyle \frac{{R}_{p}}{3.3\;{R}_{\oplus }}\right) & = & -0.209+\displaystyle \frac{1}{3}\mathrm{log}\left(\displaystyle \frac{{M}_{p}}{5.5\;{M}_{\oplus }}\right)\\ & & -0.08{\left(\displaystyle \frac{{M}_{p}}{5.5{M}_{\oplus }}\right)}^{0.4},\end{eqnarray} \tag{ 12 }$

where R_p and M_p are the radius and mass of the planetary embryo. This relation fits Mars, Venus, and Earth well.

During the simulations, we computed the mutual gravity between gas giants and embryos, but the planetesimals were not able to affect each other. This approximation was used to keep the CPU time within reasonable limits, and is justified because Jupiter clears the disk beyond 1 au in 100 kyr. Planets and planetesimals were removed once they were farther than 100 au from the Sun (whether bound or unbound) or when they collided with a planet or ventured closer than 0.2 au from the Sun.

For each permutation of the equal-mass embryo initial conditions of Jacobson & Morbidelli (2014), we ran 16 simulations (144 for each tack location). Meanwhile, for the oligarchic initial conditions, we ran 16 simulations for each starting epoch (64 for each tack location). In total, we ran 416 simulations, categorized in Table 1. The oligarchic simulations were run at the Centre for Computational Astrophysics at the National Astronomical Observatory of Japan, while the others were run at the Earth Life Science Institute at the Tokyo Institute of Technology.

Table 1. Summary of the Individual Sets of Simulations

	Equal Mass Embryos
Embryo Mass ( ${M}_{\oplus })$	${M}_{{\rm{emb}}}\;:{M}_{{\rm{pl}}}$	Tack Location	Migration Epoch (Myr)
0.025	1:1, 4:1 or 8:1	1.5 au or 2 au	0.1
0.05	1:1, 4:1 or 8:1	1.5 au or 2 au	0.5
0.08	1:1, 4:1 or 8:1	1.5 au or 2 au	1
	Oligarchic
	Migration Epoch (Myr)	Tack Location
Oligarchic	0.5	1.5 au or 2 au
Oligarchic	1	1.5 au or 2 au
Oligarchic	2	1.5 au or 2 au
Oligarchic	3	1.5 au or 2 au

Note. For each set of initial conditions, we ran 16 simulations.

Download table as: ASCII Typeset image

5. A MEASURE OF SUCCESS

According to its founders, the Grand Tack model has recorded several successes. These include, but are not limited to, the following: the ability to reproduce the mass-orbit distribution of the terrestrial planets (Walsh et al. 2011; though mostly only for Venus, Earth, and Mars), the compositional gradient and total mass of the asteroid belt (Walsh et al. 2011), the growth timescale of the terrestrial planets (Jacobson & Walsh 2015) and in some cases the angular momentum deficit (AMD), the spacing and orbital concentration of the terrestrial planets (Jacobson & Morbidelli 2014), and possibly the timing of the moon-forming impact (Jacobson et al. 2014). Some of these deserve further discussion before we outline our criteria for assigning success or failure to the individual simulations and the model itself.

Jacobson et al. (2014) use simulations of terrestrial planet formation based on the Grand Tack model to place constraints on the time of the moon-forming impact (also known as the Giant Impact or GI). They do this by requiring that, after the GI, the Earth subsequently accreted a further 0.5% of its mass, which they claim is the best estimate compatible with the fraction of highly siderophile elements in the mantle used to account for the Late Veneer (Walker 2009). From their simulations, they arrive at a timing of 95 ± 32 Myr. This is in agreement with the preferred moon-forming time from hafnium-tungsten and samarium-neodymium geochronology, although it is on the higher end (Kleine et al. 2005; Touboul et al. 2007). It is also in agreement with recent dynamical modeling and ${}^{40-39}$ Ar age compilations for the HED meteorites that likely originated from asteroid 4 Vesta (Bottke et al. 2015). In most of the simulations presented by Jacobson & Morbidelli (2014), however, the last giant impact occurs much earlier than their preferred value. O'Brien et al. (2014) pointed out that such an early impact violates the constraint posed by the amount of highly siderophile elements in the Earth's mantle used to define the Late Veneer in the first place. That said, a late accretion of 1% is still entirely within reason, and it may have been even higher: Albaréde et al. (2013) argue that upward of 4% of Earth's mass was added to the planet by the Late Veneer. Their arguments are based on the timing and amount of water delivery, and the vaporization of volatiles during accretion, throughout which a high fraction of impactor material is lost. Such a substantial amount of post-GI accretion would naturally push the epoch of the moon-forming impact much further back in time, and leaves the timing issue once again wide open. A different approach or chronometer may be needed.

Marty (2012) argues for an Earth formation time shorter than 50 Myr. In a subsequent study, Avice & Marty (2014) use iodine, plutonium, and xenon isotopic data to suggest that the closure time of the Earth's atmosphere is ${40}_{-10}^{+20}$ Myr, which would naturally coincide with the moon-forming event. The formation time suggested by Avice & Marty (2014) is in excellent agreement with the hafnium-tungsten dates of Kleine et al. (2005; 40 ± 10 Myr), but is on the lower end of that advocated by Touboul et al. (2007; 62 ${}_{-30}^{+92}$ Myr) and Halliday (2008; 70–100 Myr). In summary, the timing of the moon-forming event is still a topic of ongoing debate. It appears that radiogenic dating results in an earlier GI time than is suggested by the simulations of Jacobson et al. (2014), although the simulation results depend sensitively on the assumed amount of subsequent accretion. Most of the reported ages agree within the error bars, but the range remains at tens of millions of years.

Since individual simulations are chaotic and show a great variety of outcomes (Jacobson & Morbidelli 2014), we impose criteria that the model must adhere to, which are listed below. In what follows, we define the Mercury, Venus, Earth, and Mars analogs to have masses and semimajor axes within the ranges $(0.025\;{M}_{\oplus }\lt m\lt 0.1\;{M}_{\oplus },0.27\;\;{\rm{au}}\lt a\lt 0.5\;\;{\rm{au}})$ , $(0.4\;{M}_{\oplus }\lt {m}_{p}\lt 1.2\;{M}_{\oplus },0.55\;\;{\rm{au}}\lt a\lt 0.85\;\;{\rm{au}})$ , $(0.5\;{M}_{\oplus }\lt {m}_{p}\lt 1.5\;{M}_{\oplus },0.85\;\;{\rm{au}}\lt a\lt 1.15\;\;{\rm{au}})$ , and $(0.05\;{M}_{\oplus }\lt {m}_{p}\lt 0.15\;{M}_{\oplus },1.3\;\;{\rm{au}}\lt a\lt 1.7\;\;{\rm{au}}$ ). We require that objects in the region of the asteroid belt have their perihelia at $q\gt 1.6\;{\rm{au}}$ and their aphelia at $Q\lt 4.5\;{\rm{au}}$ . We want to add a note of caution. Since we employ initial conditions very similar to those of Walsh et al. (2011) and Jacobson & Morbidelli (2014), we expect to have a very low probability of reproducing the mass and semimajor axis of Mercury. One may then argue that we should only base our analysis on the other three terrestrial planets, but we decided against doing so.

Chambers (2001) introduced several quantities which describe the general dynamical properties of a planetary system. These are the following. First is the Angular Momentum Deficit (AMD), given by

$\begin{eqnarray}&&{\rm{AMD}}=\displaystyle \frac{{\displaystyle \sum }_{k}{\mu }_{k}\sqrt{{a}_{k}}(1-\sqrt{1-{e}_{k}^{2}})\mathrm{cos}{i}_{k}}{{\mu }_{k}\sqrt{{a}_{k}}}\end{eqnarray} \tag{ , 13 }$

where ${\mu }_{k}={m}_{k}/{M}_{\odot }$ . Second is the fraction of mass in the most massive planet (S_m). Third is a concentration parameter (S_c), given by

$\begin{eqnarray}&&{S}_{c}={\rm{max}}\;\left(\displaystyle \frac{{\displaystyle \sum }_{k}{\mu }_{k}}{{\displaystyle \sum }_{k}{\mu }_{k}{[\mathrm{log}(a/{a}_{k})]}^{2}}\right),\end{eqnarray} \tag{ 14 }$

and last is a mean spacing parameter (S_H), which is

$\begin{eqnarray}&&{S}_{H}=2\displaystyle \sum _{k=1}^{N-1}\;\displaystyle \frac{{a}_{k+1}-{a}_{k}}{{a}_{k+1}+{a}_{k}}{\left(\displaystyle \frac{{\mu }_{k+1}+{\mu }_{k}}{3}\right)}^{-1/3}.\end{eqnarray} \tag{ 15 }$

Unlike Chambers (2001), we use the mutual Hill sphere as the spacing unit. Ultimately, we invoke the following criteria to determine the success or failure of the Grand Tack model: statistically, the resulting terrestrial systems must have a lower median AMD than the current value, a concentration parameter of 107 ± 67 (2σ), a mean spacing in Hill radii of 45 ± 12 (2σ), a mass parameter of $\gt 0.5$ , they must produce at least one Mars analog (or a probability in excess of 5% for a whole set), and the most massive planet must have a greater than 5% probability of residing at 1 au, i.e., the cumulative semimajor axis distribution of the most massive planet must be lower than 0.95 at 1 au. The ranges in the listed values of S_c and S_H are computed using a Monte Carlo method adopting the range of masses and semimajor axes of our terrestrial planet analogs stated above.

6. RESULTS: TACK AT $1.5\;\mathrm{au}$

In this section, we present the results of our numerical simulations.

6.1. Equal Mass Embryos

We have run the same simulations as in Jacobson & Morbidelli (2014) in order to compare our results directly to theirs, and to determine whether a different protoplanetary disk model, realistic mass–radius relationship, and the inclusion of type 1 migration will substantially change the properties of the resulting planets. Generally, we find that most of our results are in good agreement. We do not provide a full comparison, but a systematic overview is presented and we highlight some similarities and differences.

In Figure 4, we compare the mass of each terrestrial planet that formed in our simulation versus their semimajor axis. The actual terrestrial planets are depicted as red bullets. For the most part, our results agree with Figure 1 in Jacobson & Morbidelli (2014). In our simulations, however, the peak of the distribution is situated near Venus' current location, while the region near Earth is empty in comparison. This is not the case for Jacobson & Morbidelli (2014) where the peak is in between these planets, encompasses both, and is probably caused by early type 1 migration of the embryos. The mean semimajor axis and mass of the most massive planet in our simulations are $\langle {a}_{h}\rangle =0.769\pm 0.109\;{\rm{au}}$ and $\langle {m}_{h}\rangle =0.969\pm 0.189$ ${M}_{\oplus }$ , respectively, with almost no variation within error bars as a function of either the embryo seed mass or total embryo to planetesimal mass ratio. Thus, our Venus analog is almost always more massive than the Earth analog and, with more than 95% confidence, the position of the most massive planet is inconsistent with a location at 1 au. This result is inconsistent with the findings of Jacobson & Morbidelli (2014) and we attribute this difference to the following: our use of a distinctive model for the protoplanetary disk with generally higher surface density, which causes stronger tidal damping; the inclusion of type 1 migration; and smaller planetary radii. Indeed, any material that is shepherded inward by Jupiter will be at high eccentricity and will be damped by interaction with the gas disk, which in turn will cause further inward migration since our gas disk has a higher surface density and stronger damping than in Jacobson & Morbidelli (2014). We emphasize that the inward migration caused by the damping forces is generally stronger than the direct effect of the type 1 term, so that even using the non-isothermal prescription of Paardekooper et al. (2011) would not substantially change the outcome.

**Figure 4.** Final mass of the terrestrial planets $[{M}_{\oplus }]$ vs. their semimajor axis [au]. The text above the panels indicates the embryo mass in Earth masses and the ratio of the total embryo to planetesimal mass. The beige regions indicate the range of our terrestrial planet analogs. Equal-mass embryo initial conditions with a tack at 1.5 au.
Download figure:
Standard image High-resolution image

In summary, our setup causes a peak density in solids near Venus' current position rather than in between Earth and Venus, as in Jacobson & Morbidelli (2014). For this reason, we ran another set of simulations with a tack location at 2 au rather than the typical location of 1.5 au to determine whether or not that would produce more Earth analogs. We also report a low success in producing Mercury analogs, similar to Jacobson & Morbidelli (2014), which, like them, we attribute to our initial conditions.

How do the resulting planets fare otherwise? In Figure 5, we plot the spacing parameter S_H versus the mass parameter S_m in the top panel and the concentration parameter S_c versus the normalized AMD (normalized to the current value) in the bottom panel. The typical value of S_c decreases with AMD because the higher eccentricities force the planets to be wider apart if they are to remain stable. In this figure and the ones that follow, the red symbols correspond to simulations with an initial embryo mass of 0.025 ${M}_{\oplus }$ , green to an initial embryo mass of 0.05 ${M}_{\oplus }$ , and blue to an initial embryo mass of 0.08 ${M}_{\oplus }$ . Bullets represent simulations with a total embryo to planetesimal mass ratio of 1:1, squares are for a 4:1 mass ratio, and triangles are for an 8:1 mass ratio. The large black bullet denotes the current terrestrial system while the gray bullet is for the system consisting of only Venus, Earth, and Mars. The beige regions denote 2σ regions around the mean values of S_H, S_m, and S_c. The scaled AMD range was chosen not to exceed 1 because it will increase with time (Laskar 2008; Brasser et al. 2013); the lower limit was chosen somewhat arbitrarily. Roughly 40% of all of our simulations fall in either one of the beige regions, though only 10% fall in both regions simultaneously. The results are summarized in Table 2. We reproduce the trend of Jacobson & Morbidelli (2014) wherein the AMD increases with a decrease in total planetesimal mass, though not with initial embryo mass. Jacobson & Morbidelli (2014) favor the cases with a high embryo to planetesimal mass ratio, but the high resulting AMD is inconsistent with the dynamical evolution of the terrestrial planets (Laskar 2008; Brasser et al. 2013). Jacobson & Morbidelli (2014) acknowledge that the high AMD is a potential problem, but they suggest that fragmentation during embryo–embryo collisions could produce sufficient debris to damp the AMD through dynamical friction. It is not clear whether or not this can be sustained if collisional grinding is important. Further study is needed to support or deny this claim.

Table 2. Properties of the Terrestrial Systems with a Tack at 1.5 au and Equal-mass Embryos

Embryo Mass ( ${M}_{\oplus })$	$\langle n\rangle$	$\langle {S}_{c}\rangle$	$\langle \mathrm{AMD}\rangle$	$\langle {S}_{H}\rangle$
0.025	5.0 ± 1.1	66 ± 16	2.2 ± 2.2 (1.5)	40 ± 11
0.05	4.4 ± 1.3	71 ± 29	2.6 ± 3.3 (1.5)	40 ± 11
0.08	3.7 ± 1.0	97 ± 53	1.6 ± 3.1 (0.35)	34 ± 10
${M}_{{\rm{emb}}}\;:{M}_{{\rm{pl}}}$	$\langle n\rangle$	$\langle {S}_{c}\rangle$	$\langle \mathrm{AMD}\rangle$	$\langle {S}_{H}\rangle$
1:1	4.1 ± 1.1	101 ± 47	1.1 ± 1.5 (0.37)	37 ± 13
4:1	4.8 ± 1.2	72 ± 33	2.0 ± 3.0 (1.0)	37 ± 10
8:1	4.3 ± 1.4	62 ± 19	3.3 ± 3.4 (2.11)	40 ± 10

Note. We list the average number of planets, concentration parameter, AMD, and average spacing with their standard deviations. Since the AMD distribution usually has a long tail, we list the median value in parentheses.

Download table as: ASCII Typeset image

The terrestrial system consists of four planets. We find that the average final number of planets decreases with initial embryo mass and is mostly independent of the initial total mass ratio between the planets and embryos. The relatively high number of planets with low embryo seed mass is most likely skewed by stranded embryos in the asteroid belt or near Mars' current position. The concentration parameter S_c increases with more massive embryos but decreases for a lower planetesimal mass, most likely because the AMD is higher and the planets need to be spaced farther apart to remain dynamically stable.

We define the average probability of producing a terrestrial planet analog as the fraction of planets in the designated mass-semimajor axis bin divided by the total number of produced planets. For a Venus analog, this is $24\%\pm 7\%$ , for an Earth analog it is $10\%\pm 3\%$ , and for a Mars analog it is $10\%\pm 4\%$ . All of these values are above the 5% threshold, but we see an overabundance of Venus analogs consistent with the density pileup reported earlier. We produced a total of two Mercury analogs (out of 635 planets).

Thus far, it appears that the only difference between our results and those of Jacobson & Morbidelli (2014) is our peak mass distribution being closer to the Sun than theirs, most likely because of our different disk model. Our results also differ when we investigate the timing of the moon-forming impact. Following Jacobson et al. (2014) and Jacobson & Morbidelli (2014) again, we compute the total amount of mass accreted by each planet between its last giant impact and the end of the simulation. We then plot this with a high-order polynomial best fit. The approximate timing of the moon-forming impact is the intersection of the fit and some assumed Late Veneer mass (Jacobson et al. 2014), which is at most a few percent of an Earth mass (Albaréde et al. 2013). The lower the amount of the assumed late-accreted mass, the later the giant impact had to occur because then less mass would had to have been around to impact the Earth afterward.

We plot our results in Figure 6 and obtain a best-fit value of 64 Myr for the timing of the moon-forming event assuming 1% subsequent accretion. This is in good agreement with the Hf-W results from Kleine et al. (2005) and Touboul et al. (2007), but is also sooner than what was suggested in the simulations of Jacobson et al. (2014). The difference in timing is most likely due to the fact that we consider an accreted Late Veneer mass of 1% rather than 0.5%. However, the range of mass accreted after the giant impact is rather large. From the figure, it is clear that this impact could have occurred anywhere between 30 Myr and 120 Myr, given the range of late accretion masses and uncertainties in the fit. Therefore, we do not think that our simulations, or others such as those of Jacobson et al. (2014) for that matter, can confidently predict the timing of the moon-forming event with this method.

Another issue that requires attention is the growth of Mars. Jacobson & Morbidelli (2014) conclude that it is very difficult to reproduce the rapid growth of Mars as advocated by Dauphas & Pourmand (2011). Figure 7 shows the evolution of the mass of several Mars analogs produced in our simulations as a function of time. The beige region should be avoided because the growth rate in this region is inconsistent with the Hf-W chronometer of Mars' formation (Nimmo & Kleine 2007). The blue dashed curve shows a stretched exponential growth function $m\propto {m}_{{\rm{Mars}}}(1-\mathrm{exp}[-{(t/\tau )}^{\beta }]$ ). We fit a seed mass of 0.04 ${M}_{\oplus }$ , a stretching parameter of $\beta \sim 0.5$ , and an e-folding time of $\tau \sim 10\;{\rm{Myr}}$ . These values are nearly identical to the growth of Earth and Venus (Jacobson & Walsh 2015). Some Mars analogs experience early giant impacts with other embryos, substantially increasing their mass, but even then the final growth is slow and inconsistent with the rapid growth advocated by Dauphas & Pourmand (2011), though still within limits of the Hf-W chronology of Nimmo & Kleine (2007).

**Figure 7.** Evolution of mass with time for several Mars analogs. The beige region should be avoided because the formation time is inconsistent with the Hf-W chronometer (Nimmo & Kleine 2007). The blue curve is a Weibull cumulative distribution with e-folding time $\tau =10\;{\rm{Myr}}$ , stretching parameter $\beta =0.5$ , and embryo seed mass 0.04 ${M}_{\oplus }$ .
Download figure:
Standard image High-resolution image

**Figure 7.** Evolution of mass with time for several Mars analogs. The beige region should be avoided because the formation time is inconsistent with the Hf-W chronometer (Nimmo & Kleine 2007). The blue curve is a Weibull cumulative distribution with e-folding time $\tau =10\;{\rm{Myr}}$ , stretching parameter $\beta =0.5$ , and embryo seed mass 0.04 ${M}_{\oplus }$ .
Download figure:
Standard image High-resolution image

One last thing that was not actively reported by Walsh et al. (2011), Jacobson & Morbidelli (2014), or O'Brien et al. (2014) is the amount of remaining mass in planetesimals. At the end of our simulations, we typically have a remnant mass in planetesimals of 0.051 ${M}_{\oplus }\pm 0.027\;{M}_{\oplus }$ , which is comparable to the total mass required to reproduce the Late Veneer (Raymond et al. 2013). The remnant mass depends on the original total mass ratio between planetary embryos and planetesimals. The simulations with an initial 1:1 ratio have a typical final mass of 0.08 ${M}_{\oplus }$ while the 8:1 simulations typically have 0.01 ${M}_{\oplus }$ . The decay follows a stretched exponential with best fits of $\beta =0.43\pm 0.03$ and $\beta \mathrm{log}\tau =0.40\pm 0.04$ , which are comparable to the results of Jacobson & Walsh (2015). Thus, after 150 Myr of evolution, the terrestrial planets would subsequently accrete an amount comparable to the Late Veneer. This could be problematic if the total accreted mass on the Earth after lunar formation is of the order of 1% because the subsequent accretion would overshoot the accepted 1% value, but only by a small amount.

6.2. Oligarchic Embryos

In this subsection, we report the results of Grand Tack simulations with oligarchic initial conditions. Since we do not expect substantial differences between this model and the equal-mass embryo one of Jacobson & Morbidelli (2014), we will only report on the overall results.

Figures 8 and 9 are the oligarchic equivalents of Figures 4 and 5. It appears that the correspondence with the real terrestrial system worsens as the time of the onset of migration increases. In the second figure, the red dots represent those simulations where the disk age (time of the onset of migration) is 0.5 Myr. The orange dots represent a disk age of 1 Myr, green represent 2 Myr, and blue represent 3 Myr. There are a few visible trends. First, the final AMD value tends to increase with increasing disk age. This is unsurprising because the total mass in planetesimals decreases as the disk ages, and so there is less mass to exert dynamical friction on the forming planets. Half of all systems are within the AMD-S_c boundaries in the bottom panel, but only 11% in the ${S}_{H}\mbox{--}{S}_{m}$ plot at the top, implying that only 5% fall into both regions simultaneously, which is lower than in the equal-mass embryo case. The equal-mass simulations have nearly uniform spacing anywhere from 20 to over 50 Hill radii. Systems with older disk ages are more widely spaced, while systems with younger disk ages—and therefore lower embryo seed masses and more mass in planetesimals—tend to be compact, with a typical spacing of 20 Hill radii, reminiscent of extrasolar systems (Fang & Margot 2013). The older systems also tend to have fewer planets, and these planets all appear to be of similar, sub-Venus masses because we observe a trend of a decreasing number of planets with older disk ages. The mean number of planets as a function of disk age are listed in Table 3.

**Figure 8.** Final mass of the terrestrial planets vs. their semimajor axis. The text above the panels indicates the time at which Jupiter and Saturn formed and migrated into the inner solar system. Beige regions indicate the range of our terrestrial planet analogs. Oligarchic initial conditions with a tack at 1.5 au are assumed.
Download figure:
Standard image High-resolution image

**Figure 9.** Top panel: scatter plot of the spacing parameter S_H vs. the mass parameter S_m. Red symbols correspond to simulations with a disk age of 0.5 Myr, orange to an initial disk age of 1 Myr, green dots have an initial disk age of 2 Myr, and blue ones have 3 Myr. Bottom: scatter plot of concentration parameter S_c vs. normalized AMD. Oligarchic with tack at 1.5 au.
Download figure:
Standard image High-resolution image

Table 3. Same as Table 2 for the Oligarchic Initial Conditions and a Tack at 1.5 au

Disk Age (Myr)	$\langle n\rangle$	$\langle {S}_{c}\rangle$	$\langle \mathrm{AMD}\rangle$	$\langle {S}_{H}\rangle$
0.5	3.8 ± 0.7	133 ± 26	0.32 ± 0.33 (0.18)	26 ± 9
1	3.7 ± 0.7	120 ± 28	0.23 ± 0.17 (0.18)	29 ± 6
2	3.2 ± 0.6	127 ± 71	1.73 ± 2.64 (0.47)	38 ± 9
3	2.6 ± 0.6	110 ± 64	7.1 ± 5.3 (4.3)	52 ± 13

Download table as: ASCII Typeset image

Visually, the results from the oligarchic model appear to be different from the equal-mass embryo setup, but the models are statistically nearly identical (see Table 3). We report no Mercury analogs, a probability of $32\%\pm 3\%$ for Venus analogs, $10\%\pm 5\%$ for Earth analogs, and $9\%\pm 5\%$ for Mars analogs. The location and mass of the most massive planet are $\langle {a}_{h}\rangle =0.782\;{\rm{au}}\pm 0.089\;{\rm{au}}$ and $\langle {m}_{h}\rangle =0.805{M}_{\oplus }\pm 0.160{M}_{\oplus }$ , which are on the low side for both quantities. Once again, the semimajor axis of the most massive planet is statistically inconsistent with that of Earth. In addition, unlike the equal-mass embryo case, we did not change the disk mass from one set of of simulations to the next, which could account for the lower mass of the most massive planet.

We point out that we only tested the oligarchic initial conditions for disk mass with a surface density of ${{\rm{\Sigma }}}_{0}=7$ g cm⁻² at 1 au. It is possible that a higher initial surface density would lead to higher planetary masses at the end of the simulations. That being said, it is unclear whether or not a higher surface density would increase the typical spacing between the planets because of the weak dependence of the Hill radius on the mass.

When investigating the timing of the moon-forming impact, we arrive at a time of 100 Myr, but once again the range is large, from 20 to 120 Myr. The nominal value is a little later than favored by the geochronology (Kleine et al. 2005; Touboul et al. 2007) or the plutogenic-xenon arguments of Avice & Marty (2014). In any case, the mass left in planetesimals after 150 Myr of simulation is 0.045 ${M}_{\oplus }\pm 0.029\;{M}_{\oplus }$ , which is comparable to the equal-mass embryo case. Both the timing of the moon-forming impact and the remnant mass are statistically the same as in the equal-mass embryo case. Finally, the growth of Mars proceeds similarly to that depicted in Figure 7. In summary, both the oligarchic model and the equal-mass embryo model are statistically identical within the error margins, and we cannot favor one over the other. Further study is needed to distinguish between the two.

7. RESULTS: TACK AT $2\;\mathrm{au}$

The simulations with a tack at 2 au were performed to determine whether a more distant tack location would result in the Earth analog being generally more massive than the Venus analog and whether it improves the overall fit of the model with the current architecture of the terrestrial planets. Since much of the underlying dynamics are the same, we will only report the highlights.

7.1. Equal Mass Embryos

Figure 10 is a scatter plot of the final semimajor axis and masses of the planets produced in our simulations. The wider, more massive disk will naturally produce more massive planets over a wider range of heliocentric distances, and the plot should be compared with Figure 4. There are two visual differences between the two sets of outcomes. First, the peak of the distribution is now farther out than with a tack at 1.5 au. Indeed, the mean semimajor axis of the most massive planet is at $\langle {a}_{h}\rangle =0.91\;{\rm{au}}\pm 0.19\;{\rm{au}}$ , which is much closer to the current position of Earth than with a tack at 1.5 au. The most massive planet now has a mean mass of $\langle {m}_{h}\rangle =1.15{M}_{\oplus }\pm 0.26{M}_{\oplus }$ , which is more massive than Earth but well within the error margins. This increased mass is most likely caused by the accretion annulus being wider, having been truncated at 1.25 au rather than at 1 au. The outward tail is also caused by the same effect. We explore whether this also implies that a tack at 2 au produces a better overall outcome.

First, we report that the average probability of producing a Venus analog is $17\%\pm 4\%$ , an Earth analog is $13\%\pm 5\%$ , and a Mars analog is $5\%\pm 3\%$ . We also produce some Mercury analogs ( $0.55\%\pm 0.83\%$ ). Overall, the production of Venus and Earth analogs is similar, while the production of Venus analogs was more than twice as high as Earth analogs when the tack occurred at 1.5 au. The production of Mars analogs is low, at the threshold of acceptability, caused by the fact that we generally create more massive planets near Mars' current position than with a tack at 1.5 au.

Compared to the equal-mass embryo simulations with a tack at 1.5 au, the planetary systems generated here on average have more planets. This is especially true for the cases with embryos of 0.025 ${M}_{\oplus }$ and a high mass in planetesimals. We list the number of planets and the standard deviation in Table 4. In Figure 11, we once again plot the spacing parameter S_H versus the mass parameter S_m in the top panel and the concentration parameter S_c versus the normalized AMD in the bottom panel. A similar trend of decreasing S_c with increasing AMD is visible, but not as pronounced. What is clear is that all of the systems have a concentration lower than or equal to that of the current terrestrial planets; none of them are higher (bottom panel). Since most systems have a spacing that is more compact than the current terrestrial planets (top panel), the lower concentration implies a lower variation in mass between the planets, which is generally what is observed in Figure 10. Things take a turn for the worse when we try to match the ranges of S_c, S_m, S_H, and AMD simultaneously. We find that only 3/144 cases do so, which is lower than the 5% threshold that we have adopted, and much lower than the 10% reported earlier when the tack was at 1.5 au. This low probability argues against a tack location at 2 au being suitable to reproduce the current architecture of the terrestrial planets with the equal-mass embryo initial conditions, despite the visual accuracy of the fit. We generally find that S_c is low while S_H and AMD are comparable to the case with a tack at 1.5 au, suggesting that the mass distribution is narrower and the mass variations between planets are less extreme. Most of the concentration values are low but within the acceptable range.

**Figure 11.** Top panel: scatter plot of the spacing parameter S_H vs. the mass parameter S_m. Red symbols correspond to simulations with an initial embryo mass of $0.025{M}_{\oplus }$ , green to an initial embryo mass of $0.05{M}_{\oplus }$ , and blue to an initial embryo mass of 0.08 ${M}_{\oplus }$ . Bullets are for simulations with a total embryo to planetesimal mass ratio of 1:1, squares are for a 4:1 mass ratio, and triangles are for an 8:1 mass ratio. Bottom: scatter plot of concentration parameter S_c vs. normalized AMD. Equal-mass embryo conditions with a tack at 2 au.
Download figure:
Standard image High-resolution image

**Figure 11.** Top panel: scatter plot of the spacing parameter S_H vs. the mass parameter S_m. Red symbols correspond to simulations with an initial embryo mass of $0.025{M}_{\oplus }$ , green to an initial embryo mass of $0.05{M}_{\oplus }$ , and blue to an initial embryo mass of 0.08 ${M}_{\oplus }$ . Bullets are for simulations with a total embryo to planetesimal mass ratio of 1:1, squares are for a 4:1 mass ratio, and triangles are for an 8:1 mass ratio. Bottom: scatter plot of concentration parameter S_c vs. normalized AMD. Equal-mass embryo conditions with a tack at 2 au.
Download figure:
Standard image High-resolution image

Table 4. Same as Table 2, but for the Equal Mass Embryos Initial Conditions and a Tack at 2 au

Embryo Mass ( ${M}_{\oplus })$	$\langle n\rangle$	$\langle {S}_{c}\rangle$	$\langle \mathrm{AMD}\rangle$	$\langle {S}_{H}\rangle$
0.025	5.4 ± 1.5	46 ± 8	3.2 ± 2.8 (2.6)	33 ± 8
0.05	4.5 ± 1.0	53 ± 13	2.4 ± 3.9 (0.93)	33 ± 8
0.08	4.4 ± 1.0	52 ± 15	2.1 ± 3.1 (0.74)	32 ± 8
${M}_{{\rm{emb}}}:{M}_{{\rm{pl}}}$	$\langle n\rangle$	$\langle {S}_{c}\rangle$	$\langle \mathrm{AMD}\rangle$	$\langle {S}_{H}\rangle$
1:1	5.0 ± 1.0	61 ± 13	1.0 ± 1.4 (0.36)	28 ± 7
4:1	4.9 ± 1.4	45 ± 9	2.6 ± 3.3 (1.1)	34 ± 8
8:1	4.5 ± 1.4	45 ± 9	4.1 ± 4.0 (3.1)	36 ± 7

Download table as: ASCII Typeset image

The last two issues are the timing of the moon-forming impact, which is at 60 Myr but with the same range (30 Myr to 120 Myr) as reported earlier, and a left-over planetesimal mass of $0.053{M}_{\oplus }\pm 0.04{M}_{\oplus }$ , once again on the high end but in agreement with simulations having a tack at 1.5 au. The simulations with an initial 1:1 ratio have a typical final mass of 0.1 ${M}_{\oplus }$ , while the 8:1 simulations typically have 0.02 ${M}_{\oplus }$ . The decay follows a stretched exponential with best fits of $\beta =0.44\pm 0.03$ and $\beta \mathrm{log}\tau =0.52\pm 0.05$ , which results in a slower decay than with a tack at 1.5 au and explains the higher left-over mass.

7.2. Oligarchic Embryos

In the previous subsection, we investigated whether or not a tack location at 2 au would yield a better outcome for the overall architecture of the terrestrial planets than a tack at 1.5 au in the case of equal-mass embryo initial conditions. We concluded that it appears to be difficult for this combination of tack location and initial conditions to simultaneously reproduce the combined spacing, concentration, mass distribution, and AMD of the terrestrial planets, despite generating more Earth analogs and with the heaviest planet being closer to Earth's current location. It is now worth investigating whether the oligarchic system fares any better.

Figures 12 and 13 depict the relation between the mass and semimajor axis, as well as the spacing, concentration, mass, and AMD distributions as usual. Once again, we see a broader semimajor axis-mass distribution and an overall closer spacing and lower concentration. The concentration is, however, generally a little higher than in the equal-mass embryo case. Indeed, we find that 10% of the outcomes fall within both beige regions simultaneously, which is higher than for the equal-mass embryo case and comparable to the simulations with a tack at 1.5 au. The production of terrestrial planet analogs is similar to that in the equal-mass embryo case above, so that we tend to produce more planets on average than with a tack at 1.5 au, but somewhat fewer than the equal-mass embryo case with a tack at 2 au. The final results are listed in Table 5.

**Figure 12.** Final mass of the terrestrial planets vs. their semimajor axis. The text above the panels indicates the time at which Jupiter and Saturn formed and migrated into the inner solar system. Beige regions indicate the range of our terrestrial planet analogs. Oligarchic initial conditions with a tack at 2 au.
Download figure:
Standard image High-resolution image

**Figure 13.** Top panel: scatter plot of the spacing parameter S_H vs. the mass parameter S_m. Red symbols correspond to simulations with a disk age of 0.5 Myr, orange to an initial disk age of 1 Myr, green dots have an initial disk age of 2 Myr, and blue ones 3 Myr. Bottom: scatter plot of concentration parameter S_c vs. normalized AMD. Oligarchic with tack at 2 au.
Download figure:
Standard image High-resolution image

Table 5. Same as Table 2, but for the Oligarchic Initial Conditions and a Tack at 2 au

disk Age (Myr)	$\langle n\rangle$	$\langle {S}_{c}\rangle$	$\langle \mathrm{AMD}\rangle$	$\langle {S}_{H}\rangle$
0.5	4.8 ± 0.8	73 ± 7	0.31 ± 0.14 (0.27)	25 ± 7
1	4.9 ± 0.6	72 ± 8	0.24 ± 0.11 (0.24)	24 ± 3
2	3.9 ± 0.8	53 ± 16	2.9 ± 3.6 (1.0)	38 ± 7
3	3.6 ± 0.7	48 ± 13	9.7 ± 8.8 (6.7)	43 ± 10

Download table as: ASCII Typeset image

Similarly to the equal-mass embryo case with a tack at 2 au, we produce Venus analogs $20\%\pm 3\%$ of the time, Earth analogs $13\%\pm 3\%$ of the time, and Mars analogs with a probability of $6\%\pm 3\%$ . The wider mass annulus also results in the heaviest planet having a mean semimajor axis of $\langle {a}_{h}\rangle =0.96\;{\rm{au}}\pm 0.18\;{\rm{au}}$ , but its mean mass remains a little low at $\langle {m}_{h}\rangle =0.88{M}_{\oplus }\pm 0.18{M}_{\oplus }$ , probably because we did not enhance the disk mass per set of simulations as we did for the equal-mass embryo case. Once again, the concentration value is low and we see ${S}_{H}\sim 25$ Hill radii for early disk ages and an increase in AMD with disk age.

The moon-forming impact occurs at 90 Myr, ranging from 30 to 120 Myr depending on the amount of late accretion. The mass in left-over planetesimals is comparable to earlier simulations at $0.054{M}_{\oplus }\pm 0.038{M}_{\oplus }$ .

7.3. Summary

The summary of our results is displayed in Table 6. The criteria we consider important are whether or not the model can produce Mars with a probability higher than 5%, whether or not the architecture of the system in terms of spacing, concentration, mass and AMD is consistent with the current planets, and whether or not the mass and semimajor axis of the heaviest planet are consistent with those of Earth. The moon-forming impact is a less stringent criterion because of the inherent uncertainty in the late-accreted mass and the wide range of ages that appear as a result. All that matters is that it is within the error bars of the reported hafnium-tungsten ages, which it is, but it does not provide any further constraints on the model.

Table 6. Summary of the Results of the Different Sets of Simulations

Type	Tack (au)	Mars	Architecture	$\langle {a}_{h}\rangle$	$\langle {m}_{h}\rangle$	${t}_{{\rm{moon}}}$	Mass Left ${\rm{(}}{M}_{\oplus }$ )
Equal Mass	1.5 au	✓	✓	×	✓	✓	0.05 ± 0.03
Equal Mass	2 au	✓	×	✓	✓	✓	0.05 ± 0.04
Oligarchic	1.5 au	✓	✓	×	✓	✓	0.05 ± 0.03
Oligarchic	2 au	✓	✓	✓	✓	✓	0.05 ± 0.04

Note. It is clear that a tack location of 1.5 au has difficulty reproducing the current solar system; a tack at 2 au is preferred.

Download table as: ASCII Typeset image

It appears that most initial conditions work to some extent. A tack at 2 au is able to place the most massive planet near Earth's current location, and thus is the only model to satisfy the semimajor axis constraint of the heaviest planet. This statistically rules out a tack at 1.5 au considering the constraints that we impose and the disk parameters and migration prescription that we used. However, the equal-mass embryo case with a tack at 2 au is problematic because of its poor ability to reproduce the spacing, concentration, mass, and AMD ranges simultaneously. We are confident that this low probability is inherent in the model and is not caused by sampling, and therefore we also rule it out. This leaves us with the oligarchic model with a tack at 2 au. The remnant mass in planetesimals is an obvious concern, but the somewhat low mass of the most massive planet when the tack occurred at 2 au cannot be statistically rejected. The disk age in the oligarchic model that best matches all of the constraints is 2 Myr, which is a typical time range to form the gas giants.

We point out that the above conclusions, drawn from Table 6, are valid for the disk model that we have employed. We return to its implications in the next section. One thing we have not mentioned thus far is the total mass of material that we place in the asteroid belt. This matter is discussed in detail in the next section.

8. DISCUSSION

In this study, we ran a high number of terrestrial planet simulations in the framework of the Grand Tack model with a range of initial conditions and two different tack locations. In the previous subsection, we concluded that, despite using a different model for the protoplanetary disk, the inclusion of type 1 migration and a realistic mass–radius relationship, the outcomes of our simulations are broadly similar, though each setup has its own unique pros and cons.

We have decided to use a different disk model than the traditional setup of Walsh et al. (2011) and we have outlined our reasons for doing so. We find that a tack at 1.5 au leaves too much mass near Venus' location. We attribute this to a combination of type 1 migration and inward shepherding by Jupiter, which decreases the semimajor axis of material down to below 1 au but increases the eccentricity, so that further inward migration will ensue. This raises the question as to how sensitive our results are to the choice of disk model. This can only be answered by running more simulations with varying disk parameters, which is beyond the scope of the current study. We have used a relatively hot and puffy disk with a scale height that is higher than traditional values (Hayashi 1981; Tanaka et al. 2002; Tanaka & Ward 2004), which is caused by viscous heating in the inner disk. This higher scale height causes slower embryo migration. Thus, making use of a colder, thinner disk would likely have exacerbated the overproduction of Venus analogs with a tack at 1.5 au. One way to mitigate this problem is to use a much lower surface density, as was done by Walsh et al. (2011) and Jacobson & Morbidelli (2014) but their values seem artificial. A much lower surface density would decrease the migration rate of Jupiter and Saturn, and even though their migration rate does not appear to affect the final orbital architecture of the terrestrial system (Walsh et al. 2011), it does increase the difficulty for these planets to reach their final positions beyond 5 au (D'Angelo & Marzari 2012). Thus, even when using a colder, less-massive disk, we still expect to see an overproduction of Venus analogs when the tack occurred near 1.5 au.

A second topic concerns the timing of the moon-forming impact. We generally find agreement between our typical time of 60–90 Myr and geochronology. However, our uncertainties are typically 30 Myr or longer, suggesting that GI occurred anywhere from 30 Myr to 120 Myr, which is what we have claimed in the previous sections. It is debatable whether this range implies anything meaningful. It is consistent with the value 95 ± 32 Myr reported in Jacobson et al. (2014), even though they ran their simulations for a little longer and used a much smaller amount of late-accreted mass. In summary, we do not think that our simulations, or those of Jacobson et al. (2014), can say anything meaningful about the timing of the GI beyond what is known from geochronology.

Another issue that requires discussion is the mass left over in planetesimals after planet formation. This is typically 0.05 ${M}_{\oplus }$ but can be as high as 0.1 ${M}_{\oplus }$ . This left-over mass has implications for the cratering rates on Noachian Mars and the Pre-Nectarian moon. The most efficient way to eliminate this material is through collision with the terrestrial planets. Ejection by the giant planets or collisions with the Sun is much more difficult. Preliminary simulations of this population of planetesimals indicate that it decays slowly, following a stretched exponential with a stretching parameter of $\beta \sim 0.83$ and an e-folding time of τ ∼ 85 Myr. A slow decay is preferred by lunar cratering records (Werner et al. 2014), but a slower decay is necessary so as not to have late melting of the crust of the planets (Abramov et al. 2013, Abramov & Mojzsis 2016).

One solution may be for these planetesimals to grind themselves to dust and subsequently be lost through Poynting-Robertson drag or radiation pressure. The difficulty with this idea is that the high ratio of highly siderophile elements in the Earth and moon suggest that the Late Veneer impactors were large (around 2000 km; Bottke et al. 2010). If these impactors were large, then there is no reason to believe that the impactors after the Late Veneer were substantially smaller. A simple argument is that Ceres is the only 1000 km body in the asteroid belt, and with a typical implantation probability of 0.1% (Walsh et al. 2011), there should have been at least 1000 Ceres-sized bodies. With of the order of 5% of the total mass remaining after 150 Myr, we have 50 Ceres-sized bodies still present, with perhaps 10 bodies the size of 4000 km. Since it was probable that the size distribution of the remnant planetesimals was shallow (Bottke et al. 2010), most of the mass is in the large bodies, and thus we consider it very unlikely that this mass was ground down by collisional erosion.

A quick estimate of the collisional timescale can be made with an $n\sigma {vt}=1$ argument, where n is the number density of planetesimals, v is their typical encounter velocity, and $\sigma =\pi {r}^{2}$ is their collisional cross section. The number density $n=M/(2{\pi }^{2}{a}^{2}\;m{\rm{\Delta }}a\mathrm{sin}i)$ , where M is the total mass of remnant planetesimals, ${\rm{\Delta }}a$ is the width of the annulus in which the planetesimals are situated and m is their individual mass. We then have

$\begin{eqnarray}&&\displaystyle \frac{3{Mvt}}{8{\pi }^{2}{a}^{2}{\rm{\Delta }}a\mathrm{sin}i\rho r}=1.\end{eqnarray} \tag{ 16 }$

From our simulations we find that the typical semimajor axis of these planetesimals is 1.5 au and ${\rm{\Delta }}a\sim 0.5\;{\rm{au}}$ , $i\sim 20^\circ$ and using a planetesimal density of $\rho =3000$ kg m⁻³ and encounter velocity $v\sim {{ev}}_{K}\sim 12$ km s⁻¹ we have for a planetesimal of radius 500 km that $t\gt 1\;{\rm{Gyr}}$ . This is longer than the estimate of 56 Myr of Raymond et al. (2013) because they consider a smaller annulus and a planar problem. Thus for large planetesimals, collisional grinding is most likely unimportant, though future studies need to verify or deny this claim. In any case, the amount of left-over material warrants a separate investigation, in particular if the size distribution is shallow and most mass is in large planetesimals. Its results will be discussed in a companion paper.

Another topic that was not discussed in the previous sections was the amount of mass that is placed in the asteroid belt. With our definition of the asteroid belt region (perihelion $q\gt 1.6\;{\rm{au}}$ and aphelion $Q\lt 4.5$ au), we find that simulations with a tack at 1.5 au place roughly 0.6% of material in the asteroid belt, while this increases to 0.9% when the tack is at 2 au, which is comparable to but slightly in excess of the percentage reported by Walsh et al. (2011). The main reason our simulations with a tack at 1.5 au have a higher amount of mass in the asteroid belt than theirs is due to the stronger gas drag acting on planetesimals in the beginning of the simulations because the surface density of our disk is higher than theirs. The higher value with the tack at 2 au is clearly the result of the weaker sculpting of the disk by Jupiter. All of these simulations leave us with an asteroid belt whose mass is at least an order of magnitude higher than the current value, with the caveat that we only used planetesimals with a size of 50 km for the purpose of the gas drag. Larger planetesimals would have reduced the remnant mass in the asteroid belt while a smaller size would have increased it (Matsumura et al. 2016). It has been suggested that chaotic diffusion over 4 Gyr and giant planet evolution deplete the belt by approximately 75% of its mass, but then the remaining amount of mass is still inconsistent with what is observed today (Minton & Malhotra 2010). Then again, our simulations do not take collisional erosion into account, which could erode the belt even further (Bottke et al. 2005).

A final topic that requires discussion is the total number of planets. When taking each set of initial conditions as a whole, the total number of planets is 4.4 ± 1.3 for the equal-mass embryo case with a tack at 1.5 au, 3.3 ± 0.8 for the oligarchic case with a tack at 1.5 au, 4.8 ± 1.3 for the equal-mass embryo case with a tack at 2 au, and 4.3 ± 0.9 for the oligarchic case with a tack at 2 au. All of these are consistent with four terrestrial planets at the end. However, most of our simulations do not produce a Mercury analog. This is a result of the initial conditions that we have employed, and it is noteworthy that the formation of Mercury has been mostly ignored in earlier works (Walsh et al. 2011; Jacobson & Morbidelli 2014; Jacobson et al. 2014; O'Brien et al. 2014; Jacobson & Walsh 2015). If we are to ignore it too, then the total number of planets we must produce is three. All of the models as a whole are statistically consistent with three terrestrial planets, though some subsets within the four models are not. Since the formation and evolution of Mercury are currently unknown, further study is needed to rule out whether the oligarchic model with a tack at 2 au can be made consistent with the current inner solar system. It is possible that, by the end of our simulations, the final system is not stable in the long-term and a very late collision could remove another planet. This is inconsistent with the historical evolution of the Solar System, and so we have to take the number of planets at the end of the simulations as final.

9. CONCLUSIONS

We have investigated the dynamical formation of the terrestrial planets in the framework of the Grand Tack scenario. It has been claimed that the Grand Tack reproduces several observed features of the inner Solar System that previous models failed to do, such as the low mass of Mars and the compositional gradient in the asteroid belt. We examined this scenario in more detail here but applied a different disk profile, a realistic mass–radius relationship, and took into account type 1 migration. We have stated our reasons for doing so in Section 2, and performed sensitivity tests to determine whether any of these differences matter. The answer appears to be "yes": with the initial conditions and disk and migration model that we employed, we statistically ruled out a tack at 1.5 au because we produced an excess of Venus analogs and a deficit of Earth analogs. We attribute this excess of Venus analogs to our disk model because it has a higher surface density than that of Walsh et al. (2011). Additionally, with more than 95% confidence, the semimajor axis of the most massive planet is inconsistent with Earth's location, while upon visual inspection of their results the same is not true in the simulations of Jacobson & Morbidelli (2014). Thus, the inclusion of type 1 migration and a much higher initial disk surface density, together with smaller radii of the planets, serve to shift the mass distribution closer to the Sun. This calls for a more distant tack location. We find that the model that best matches the current architecture of the terrestrial planets has a tack at 2 au and oligarchic initial conditions.

We thank Kevin Walsh for making available to us his version of SyMBA that incorporates the migration of the gas giants and the gas profile, Hal Levison for indicating that we should include a realistic mass–radius relationship, and an anonymous reviewer for constructive comments. R.B. is grateful for financial support from the Astrobiology Center Project of the National Institutes of Natural Sciences (NINS) grant number AB271017, and to the Daiwa Anglo-Japanese Foundation for a Small Grant. R.B. and S.J.M. acknowledge the John Templeton Foundation—Ffame Origins program in support of CRiO. S.J.M. is grateful for support by the NASA Exobiology Program (NNH14ZDA001N-EXO). S.C.W. is supported by the Research Council of Norway (235058/F20 CRATER CLOCK) and through the Centres of Excellence funding scheme, project number 223272 (CEED). Numerical simulations were in part carried out on the PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.

ANALYSIS OF TERRESTRIAL PLANET FORMATION BY THE GRAND TACK MODEL: SYSTEM ARCHITECTURE AND TACK LOCATION

Article metrics

Permissions

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. DEVIATIONS FROM THE ORIGINAL GRAND TACK MODEL

2.1. Protoplanetary Disk

2.2. Embryo Migration

3. INITIAL CONDITIONS

3.1. Equal Mass Embryo Initial Conditions

3.2. Oligarchic Initial Conditions

4. METHODOLOGY

5. A MEASURE OF SUCCESS

6. RESULTS: TACK AT $1.5\;\mathrm{au}$

6.1. Equal Mass Embryos

6.2. Oligarchic Embryos

7. RESULTS: TACK AT $2\;\mathrm{au}$

7.1. Equal Mass Embryos

7.2. Oligarchic Embryos

7.3. Summary

8. DISCUSSION

9. CONCLUSIONS

ANALYSIS OF TERRESTRIAL PLANET FORMATION BY THE GRAND TACK MODEL: SYSTEM ARCHITECTURE AND TACK LOCATION

Article metrics

Permissions

Share this article

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. DEVIATIONS FROM THE ORIGINAL GRAND TACK MODEL

2.1. Protoplanetary Disk

2.2. Embryo Migration

3. INITIAL CONDITIONS

3.1. Equal Mass Embryo Initial Conditions

3.2. Oligarchic Initial Conditions

4. METHODOLOGY

5. A MEASURE OF SUCCESS

6. RESULTS: TACK AT 1.5\;\mathrm{au}

6.1. Equal Mass Embryos

6.2. Oligarchic Embryos

7. RESULTS: TACK AT 2\;\mathrm{au}

7.1. Equal Mass Embryos

7.2. Oligarchic Embryos

7.3. Summary

8. DISCUSSION

9. CONCLUSIONS

6. RESULTS: TACK AT $1.5\;\mathrm{au}$

7. RESULTS: TACK AT $2\;\mathrm{au}$