© ESO 2018

1 Introduction

Basic details of the physiochemical dynamics which led to the formation of the terrestrial planets have remained elusive until relatively recently. Moore’s Law improvements in computational performance have transformed the way dynamicists are able to study the physical and chemical processes of terrestrial planet accretion via increasingly sophisticated numerical simulations. Analyses of terrestrial rocks and meteorites are being made with ever-more powerful analytical tools that yield data at unprecedented precision and accuracy. We now, therefore, have a firmer understanding of the chemical and isotopic compositions of Earth, the Moon, Mars, Vesta, various asteroidal meteorites derived from some 100–150 parent bodies, and some comets. These data help to constrain both the absolute and relative chronology of planet formation that is subsequently fed into modern dynamical models of accretion (e.g. Morbidelli et al. 2012).

Apart from Earth, Mars is the only other planet for which we have direct samples. Geochemical and isotopic analysis of these samples allows us to narrow the uncertainties of the timing of Mars’ formation. For example, Dauphas & Pourmand (2011) analysed a suite of martian meteorites to argue that Mars reached about half of its current size in approximately 1.8 Myr. Mars likely finished its accretion within ~10 Myr of the start of the solar system. It is straight forward to think that Earth requires more time to form because of its much larger mass. Indeed, derivation of timescales for terrestrial accretion from the Hf-W chronometer supports this idea by showing that Earth finished its core formation at least 30 Myr after the calcium–aluminium-rich inclusion (CAI) formation (e.g. Kleine et al. 2009 and references therein). These observations can be summarised as a general time line for the formation of the terrestrial planets, and combining such geochemical data with numerical simulations helps in efforts to sort out the mystery of the formation of the sampled terrestrial planets and by extension the whole of the inner solar system.

2 Background

2.1 The final stage of planet formation

According to traditional dynamical models of the late stage of planet formation, coagulation of planetesimals into planetary embryos and subsequent giant impacts between embryos gave rise to the terrestrial planets. Several different models have been proposed forthis scenario. The foundational one is the Classical model, in which all the terrestrial planets and the two gas giants formed near their current locations. Chambers (2001) performed simulations with different masses of planetesimals and embryos to study thegrowth of terrestrial planets with Jupiter and Saturn on their current orbits. Simulations showed that it is typical to produce three to four terrestrial planets with orbital characteristics similar to the current system. A recurring problem with the simulations, however, is that the masses of computed planets which form close to 1.5 AU are consistently greater than the current mass of Mars. This was recapitulated in simulations by Raymond et al. (2009), especially when both Jupiter and Saturn wereprovided circular orbits. To avoid forming a massive Mars, an alternative solution was proposed by Agnor et al. (1999) and expanded upon by Hansen (2009); this work initially confined all solid material to the area in between the current locations of Venus and Earth at 0.7 and 1 AU in an annulus model. With this setting, Hansen (2009) successfully reproduced the mass–distance relation of the terrestrial planets. If all the material making up the terrestrial planets was confined within a small annulus, what is the underlying mechanism that leads to this configuration? The Grand Tack model proposed by Walsh et al. (2011) provides a reasonable explanation.

In the Grand Tack, Jupiter first formed in a position closer to the Sun than its current orbit and opened a gap in the gas disk. Owing to the unbalanced tidal torque acting on the planet, Jupiter migrated inwards through type-II migration (Lin & Papaloizou 1986). At the time Jupiter was migrating, Saturn formed further away and slowly accreted its gaseous envelope. Saturn then began migrating inwards once it reached about 50 Earth masses, caught up with Jupiter, and the two gas giants were trapped in the 2:3 mean motion resonance (Masset & Snellgrove 2001). With this specific orbital spacing and mass ratio between the two gas giants, the total torque on both planets reversed and therefore the migration direction also reversed outward from the Sun. The tack location (the location where Jupiter reversed its direction of migration) is set to be 1.5 AU in order to confine the growth of Mars (Walsh et al. 2011). This mechanism shows how the migration of Jupiter reshaped the structure of the inner solid disk. Two possible scenarios emerge from this model: either Jupiter scattered the planetesimals and embryos in its path away, or pushed them into the inner region by trapping them in mean motion resonance. In either case, when Jupiter reversed its migration, the region in between 1 and 1.5 AU was cleared. On the other hand, solid material piled up within 1 AU, thus duplicating the outer edge configuration suggested by the annulus model. Therefore, the mass–distance relation of the terrestrial planets and the low mass of Mars are explicable by the Grand Tack model.

Another possible mechanism creating an annulus near the current orbits of the terrestrial planets is proposed by Drążkowska et al. (2016). In this model, pebbles pile up in radial distribution due to a combination of growth and radial drift of dust and therefore a narrow annulus of planetesimals close to 1 AU can be formed through the streaming instability (Johansen et al. 2007; Ida & Guillot 2016). This model implies that all the terrestrial planets form from the same building blocks and therefore that (see below) their bulk and isotopic composition should be similar to each other. However, the isotopic compositions suggest that this may not be the case.

2.2 Bulk compositional differences between Earth and Mars

The bulk isotopic compositions of Earth and Mars are different. The strongest evidence for this is the three oxygen isotope system (expressed as Δ¹⁷O_VSMOW = δ¹⁷O_VSMOW – 0.52δ¹⁸O_VSMOW, where δ-notation denotes deviations in parts-per-thousand where isotopic ratios are normalised to the standard mean ocean water and Δ expresses the deviation of these coupled ratios from the terrestrial mass-fractionation line of slope m = 0.52). Bulk Mars is enriched in the minoroxygen isotope (¹⁷O) with respectto terrestrial and lunar values (e.g. Franchi et al. 1999; Rubin et al. 2000; Mittlefehldt et al. 2008; Agee et al. 2013; Wittmann et al. 2015). Other isotopic systems that trace nucleosynthetic anomalies, such as Titanium (ε⁵⁰ Ti), Chromium (ε⁵⁴ Cr), and Nickel (ε⁶⁴Ni), likewise show clear differences between average terrestrial and martian values (Brasser et al. 2017). No known mass-dependent process can account for these deviations (Qin & Carlson 2016); they are instead attributable to implantation of neutron-rich isotopes by nearby supernovae at the time of solar system coalescence (Qin et al. 2010). This nucleosynthetic contamination is widely assumed to have been more influential in the composition of outer regions of the disk. Although under debate, a number of studies have proposed that the heliocentric composition gradient of the solid disk can be mapped by the isotopic differences of nucleosynthetic tracers that were implanted in different locations within the planet-forming disk. Taken together, it has been proposed that isotopic compositional differences between Earth and Mars hint that they were formed with different mixtures of source components (Wänke & Dreibus 1988, 1994; Lodders 2000; Warren 2011), with the further implication that source locations for the two worlds were different (cf. Fitoussi et al. 2016).

The difference between the source location(s) for components that gave rise to Earth and Mars is lent further credence by several recent isotope studies in terrestrial samples and martian meteorites. Assuming that the feed-stock for accretion of both Earth and Mars was dominated by chondritic material, analysis suggests that Earth should accrete about 70% of enstatite chondrite, about 25% of ordinary chondrite and less than 5% of carbonaceous chondrite in order to reproduce the well-documented isotopic compositions of ¹⁷O, ⁴⁸C, ⁵⁰Ti, ⁶²Ni, or ⁶⁴Ni, ⁹²Mo, and ¹⁰⁰Ru (Dauphas et al. 2014; Dauphas 2017; Fischer-Gödde & Kleine 2017). Mars, on the other hand, has a larger uncertainty in its bulk composition. Its bulk composition is suggested to be 45% enstatite chondrite and 55% ordinary chondrite to account for the planet’s documented ¹⁷O, ⁵⁰Ti, ⁵⁴Cr, ⁶²Ni, and ⁹²Mo isotopic values as determined from the martian meteorites (Sanloup et al. 1999; Tang & Dauphas 2014). With a more rigorous analysis of uncertainties, Brasser et al. (2018) recorded 29–68% entatite chondrite and 32–67% ordinary chondrite in their analysis, which is in agreement with previous studies within uncertainties. The generally much higher portion of ordinary chondrite in Mars as compared to Earth underscores the differences in bulk composition between the two planets. Other work suggested that the highly reduced and volatile depleted enstatite chondrites were formed in the inner region of the solar system (closer than 1.5 AU), and more oxidised and volatile-rich ordinary chondrite sources originated in a region further away (probably >1.5 AU; Morbidelli et al. 2012; Rubie et al. 2015; Fischer-Gödde & Kleine 2017). Consequently, Brasser et al. (2017) argued that the formation region of Mars is likely to be further from the Sun than Earth owing at least to the fact that its composition has a much higher computed fraction of ordinary chondrite.

3 Framework

Assuming that Mars formed distantly from Earth’s feeding zone, we first investigate the possibility of producing differential Earth and Mars isotopic compositions with both the Classical and the Grand Tack model from the results of our N-body simulations. This is accomplished via comparison of the accretion zone of the Earth and Mars analogues generated in the simulation as described in the following section. Next, we estimate the average percentage of enstatite and ordinary chondrite incorporated into Earth and Mars by assuming the initial solid disk is composed of only enstatite in the inner region and ordinary chondrite in the outer region (see Sect. 5). Finally, we calculate the expected isotopic composition of ¹⁷O, ⁵⁰Ti, ⁵⁴Cr, ¹⁴²Nd, ⁶⁴Ni, and ⁹²Mo in the martian mantle (Sect. 6). Our purpose in this study is to examine how the orbital dynamics affects Mars’ isotopic composition.

4 Comparison of accretion zones between Earth and Mars

4.1 Simulation set up

To estimate the possibility of reproducing Mars with a composition that is different from Earth’s, we analyse the data from a high number of N-body simulations for the Grand Tack and the Classical model as described above. The simulations and data used to examine the Grand Tack model were previously discussed in Brasser et al. (2016), and we refer the reader to that work. We first briefly summarise their initial conditions and the organisation of the simulations.

In Brasser et al. (2016), the formation of terrestrial planets were studied in the Grand Tack scenario by focusing on the traditional dynamic perspectives, which are the final orbital architecture and the masses of the terrestrial planets. The tack location for Jupiter at both 1.5 and 2 AU was considered, where “tack” location is the position at which Jupiter reversed its migration. Two types of initial condition for planetestimals and embryos were adopted: the equal-mass initial condition, and the oligarchic initial condition (embryo resulting from the oligarchic growth of planetesimals). Brasser et al. (2016) employed the initial condition of embryos and planetesimals from Jacobson & Morbidelli (2014) for the equal-mass initial condition, in which embryos in the same simulation have the same initial masses. The initial embryo masses are either 0.025, 0.05, or 0.08 M_⊕, where M_⊕ is the mass of Earth. All the embryos within the disk were embedded in a disk of planetesimals. The total mass ratio of the embryos to planetesimals is either 1:1, 4:1, or 8:1. In total there are nine different subsets of initial conditions for the equal-mass initial condition. They performed 16 simulations in each subset; in total 144 simulations for each tack location. Table 1 shows the initial number of embryos in these nine subsets of initial conditions. The number of embryos increases with the mass ratio of embryos to planetesimals but decreases with the masses of embryos. They kept the number of planetesimals at 2000. The initial densities of embryos and planetesimals are 3 g cm ⁻³.

In addition to forming terrestrial planets with equal-mass embryos, Brasser et al. (2016) also ran simulations with the initial conditions that are computed from the traditional oligarchic growth of the planetesimals (Kokubo & Ida 1998). Kokubo & Ida (1998) found that the orbital repulsion between embryos keeps their mutual separation wider than about 5 R_Hill during the oligarchic growth of the planetesimals, where R_Hill is the mutual Hill radii of two adjacent embryos. The oligarchic initial condition is more realistic compared to the equal mass initial condition because it is a result from the actual study of planetary embryo formation. Brasser et al. (2016) applied the semi-analytical oligarchic approach of Chambers (2006). First, they calculated the total solid mass between 0.7 and 3 AU according to the minimum mass solar nebula Σ_s = 7 g cm⁻²(a/1 AU)^−3∕2 (Hayashi 1981). Second, they increased the solid density at the ice line by a factor of three following the study of Ogihara & Ida (2009). The ice line was assumed to be static at 2.7 AU. Third, they set the spacing between the embryos to 10 R_Hill based on the results of Kokubo & Ida (1998); although in our work we ran additional runs with 7 and 5 R_Hill. The spacing was computed assuming that embryos had their isolation masses, m_iso = 2πabΣ_s, where a is the semi-major axis of embryo and b is the spacing between adjacent embryos (Chambers 2006). Here m_iso represents the largest mass of an embryo if it accretes all the solid mass within the annulus 2πab. The semi-major axis of embryo n is a_n = a_n−1 [1 + b(2m_iso∕3 M_⊙)], which nearly follows a geometric sequence, where M_⊙ is the solar mass.

Other than calculating the spacing and the semi-major axis of embryos, the initial masses of embryos also needed to be computed. Following Chambers (2006), the mass of embryos increases up to their isolation mass as (1)

where τ is the growth timescale, which depends on the semi-major axis a of the embryo, embryo spacing b, solid surface density Σ_s and the radii of planetesimals that assembled into the embryo, and t is the age of the solid disk before the start of the migration of Jupiter. Chambers (2006) assumed a planetesimal size of 10 km in calculating τ, which was also adopted by Brasser et al. (2016). By substituting t = 0.5, 1, 2, or 3 Myr into Eq. (1), they obtained four different subsets of initial conditions for the oligarchic growth scenario. They performed 16 simulations for each t, 64 in total for each tack location of Jupiter with 10 R_Hill embryo spacing. We perform more simulations with 5 and 7 R_Hill embryo spacing (64 each), but only for tack location of Jupiter at 1.5 AU. Figure 1 shows the initial masses of the embryos as a function of their initial semi-major axes for the oligarchic initial condition. Different colors correspond to different ages of the solid disk. By comparing the squares with different colors, we observe that embryos with similar initial semi-major axes have larger masses if the solid disk is older. This is because embryos accrete more planetesimals within the disk if the solid disk exists for a longer time before the start of Jupiter’s migration. The initial number of embryos is larger in the 1 and 2 Myr disk (see Table 2). We also compare the initial number of embryos with different initial R_Hill separation in Fig. 1, but only for the 0.5 Myr disk. Comparing the red points with different shapes, the mass of embryos is lower but the initial number of embryos is larger if they have smaller mutual separation (see Table 2).

The system consisting of gas giants, planetary embryos, and planetesimals was then simulated with the symplectic integrator package SyMBA (Duncan et al. 1998) for 150 Myr using a time step of 7.3 days. The migration of Jupiter and Saturn was mimicked through the fictitious forces (Walsh et al. 2011). The gas disk developed by Bitsch et al. (2015) was adopted, which is different from the one adopted by Walsh et al. (2011) and has a higher surface gas density. The initial gas surface density profile scales as Σ(r) ∝ r^−1∕2 and the temperature profile is T(r) ∝ r^−6∕7, where r is the radial distance to the Sun. Both gas surface density and temperature decay until 5 Myr after the start of simulation, after which the disk is artificially photo-evaporated away with the e-folding time 100 kyr (Bitsch et al. 2014). The existing gas disk caused type-I migration and tidal damping of the eccentricities of the planetary embryos. The planetesimals are affected by the gas drag, for which Brasser et al. (2016) assumed each planetesimal had a radius of 50 km.

We also study the Classical model, but only for the oligarchic initial condition due to the limited computational resources. In other words, we used the same oligarchic initial conditions but placed Jupiter and Saturn on their current orbits. The time step and the total simulation time for each simulation in the Classical model is the same as the Grand Tack model. We test the Classical model with 10 and 5 R_Hill mutual spacing between embryos.

Fig. 1 Initial masses of embryos against their initial semi-major axes for the oligarchic initial condition. The red points are embryos in a 0.5 Myr disk, green points are in a 1 Myr disk, blue points are in a 2 Myr disk, and black points are in a 3 Myr disk. Different point shapes represent different separation between adjacent embryos. Squares depict embryos separated by 10 RHill from each other, triangles are embryos with spacing of 7 RHill and rhombuses are embryos with spacing of 5 RHill. The comparison of different spacing between embryos is only plotted for the 0.5 Myr disk.

4.2 Definition of accretion zone

We compare the composition of the planets that formed in the current region of the terrestrial planets. Assuming the disk follows a simple relation between isotopic ratios and distance to the Sun (see Sect. 2.2), the bulk composition of each planet can be roughly quantified. The compositional differences between planets are quantified by tracking the initial heliocentric distance of the material (either embryos or planetesimals) accreted into each planet. Therefore, the mass-weighted mean initial semi-major axis of incorporated material, (2)

was calculated for each terrestrial planet analogue, where m_i and a_i are the mass and initial semi-major axis of the accreted object i and N is the total number of embryos and planetesimals incorporated into the planet (Kaib & Cowan 2015; Brasser et al. 2017; Fischer et al. 2018). Planets with similar a_mwmi should have similar bulk composition. The weighted standard deviations of a_mwmi, (3)

are also calculated for each planet analogue. a_mwmi ± σ_w provides us the accretion zone of the planet analogue, which is the main region of the disk that the planet analogue sampled; σ_w is defined as zero if the planet did not accrete any planetesimal.

4.3 Results and discussion

4.3.1 Classical: Oligarchic initial condition

Figure 2 depicts the a_mwmi of Venus analogues (green), Earth analogues (blue), and Mars analogues (red) as a function of their final semi-major axes for the oligarchic initial condition of the Classical model with Jupiter and Saturn on their current orbits. The results of initial embryos with spacing of 10 R_Hill and 5 R_Hill spacing, as well as different disk ages t are combined here. Venus, Earth, and Mars analogues are defined as having masses and semi-major axes within the ranges (0.4 M_⊕< m_p < 1.2 M_⊕, 0.55 AU < a < 0.85 AU), (0.5 M_⊕< m_p < 1.5M_⊕, 0.85 AU < a < 1.15 AU), and (0.05 M_⊕ < m_p < 0.15M_⊕, 1.3 AU < a < 1.7 AU; Brasser et al. 2016). The error bars are ±1σ_w, which represent the accretion zone of each planet.

Comparing the data in Fig. 2, we find that most of the Earth (blue squares) and Venus (green squares) analogues have 1 AU < a_mwmi < 1.5 AU and their error bars share a similar range, which implies that the accretion zones, and therefore their bulk compositions, are likely to be similar. There is a weak trend of planets with larger final semi-major axes having larger a_mwmi. A similar but stronger trend has been observed in the simulations performed by Kaib & Cowan (2015). One possible reason is that our initial solid disk truncates at 0.7 AU, while the inner boundary of Kaib & Cowan’s disk is at 0.5 AU. If the initial solid disk is extended to a region closer to the Sun, the Venus analogues should accrete more material originated from the inner-most region of the disk. Therefore, most of the Venus analogues would sample the disk closer to the Sun than the Earth analogues do and generate a more obvious linear trend between a_mwmi and the final semi-major axis. Comparing Earth (blue squares) with Mars (red squares) in Fig. 2 leads us to a similar conclusion, since six out of eight of the Mars analogues have 1 AU < a_mwmi < 1.5 AU and the error bars of all the red and blue squares share a similar range. Therefore, the Classical model predicts a similar bulk composition between Venus, Earth, and Mars. That is to say, the Classical model inadequately explains the bulk composition difference between Earth and Mars if we are able to form Mars with its current mass. If this restriction is ignored, it is possible to form a distinct Earth and Mars (see Sect. 5).

Table 3 lists the statistics for all sets of simulations performed in this study. The number of Venus and Earth analogues we obtained from the Classical model is similar to the oligarchic simulation of the Grand Tack model. Although 102 planets are formed in the orbital parameter space we set for Mars analogues, only 8 fall within a satisfactory mass range for Mars within the 128 Classical’s simulations performed. This is far fewer than those obtained from the Grand Tack simulations. The reason is that most of the planets formed near Mars’ current position have masses larger than the value that we defined for Mars analogues (>0.15 M_⊕). The average mass of the planets within 1.3 < a < 1.7 AU in the simulations is 0.45 M_⊕, which is more than three times the mass of Mars. As previously mentioned, the biggest problem of the Classical model is forming a Mars that is too massive (e.g. Chambers 2001; Raymond et al. 2009).

Fig. 2 Mass-weighted mean initial semi-major axis (amwmi) of material accreted into each Venus analogue (green), Earth analogue (blue), and Mars analogue (red) as a function of their final semi-major axis for the oligarchic initial condition of the Classical model. We combine the results of 10 and 5 RHill and differentages of the solid disk (either t = 0.5, 1, 2, or 3 Myr) within this plot.

4.3.2 Grand Tack: Oligarchic initial conditions

Figure 3 resembles Fig. 2, but depicts the results from the Grand Tack simulations with oligarchic initial conditions. The results of the tack at 1.5 and 2 AU of Jupiter, as well as different adjacent embryo spacing and disk ages, are all plotted in Fig. 3. Similar to Fig. 2, most of the Venus and Earth analogues have a_mwmi between 1 and 1.5 AU. A few of them have a_mwmi smaller than 1 AU; these sample the innermost portion of the disk. Statistics of the Mars analogues suggest that their bulk composition is likely to be similar to both Earth and Venus since only ~13% of the Mars analogues within the oligarchic simulations have a_mwmi > 1.5 AU, and their error bars in Fig. 3 mostly share the same range with the Venus and the Earth analogues. Only two Mars analogues in our output have a_mwmi > 2 AU. We return to this type of Mars analogue below (Sect. 4.3.3). According to simulations, the bulk compositions of these three planets are likely to be similar.

It is worth pointing out that about one-fourth of the Mars analogues have a_mwmi ≤ 1 AU. The fraction is even higher than the fraction of the Earth analogues with a_mwmi ≤ 1 AU. Therefore, the oligarchic simulation results show that it is possible for Mars to be formed in the innermost region of the protoplanetary disk and then scattered to its current location. This was also suggested by Hansen (2009) and Fitoussi et al. (2016) and might be taken as a natural explanation for the formation of Mars, that is if we do not consider its isotopic compositional difference from Earth (Brasser et al. 2017).

Reducing the spacing of embryos (from 10 R_Hill to either 7 or 5 R_Hill) produces different sets of initial conditions with lower average individual embryo mass and a higher initial number of embryos (see Fig. 1 and Table 2). Table 3 shows the statistics of 5 and 7 R_Hill spacing for the oligarchic Grand Tack simulations. The statistics of both 5 and 7 R_Hill spacing are similar to that of 10R_Hill, which indicates that the final results are not strongly dependent on the initial conditions of the embryos.

Fig. 3 As in Fig. 2, but for the Grand Tack model with oligarchic initial conditions. The results of Jupiter tacked at 1.5 or 2 AU, as well as different initial spacing between adjacent embryos (either 10, 7, or 5 RHill) and differentages of the solid disk (either t = 0.5, 1, 2, or 3 Myr) are combined.

4.3.3 Grand Tack: Equal mass initial condition

Figure 4 is similar to Fig. 3, but shows the result of using the equal-mass initial conditions of the Grand Tack model where Jupiter tacked at either 1.5 or 2 AU. This plot resembles the previously published Figs. 2 and 8 of Brasser et al. (2017), but we highlight the terrestrial planet analogues in our plot.

Similar to the oligarchic initial condition, we find that nearly all of our Venus, Earth, and Mars analogues have error bars that share a similar range of uncertainty. As such, most of the terrestrial planet analogues can be seen to have sampled a similar region of the disk. The value a_mwmi needs to be >2 AU if the error bars of the Mars analogues are not in the same range as all of the Venus and Earth analogues. The Mars analogues with a_mwmi > 2 AU are the ones that we are searching for because those Mars analogues would have a very different accretion zone with the Venus and Earth analogues. They formed in the outer region of the disk and accreted material mainly from there. Then Jupiter scattered them close to the current location of Mars when it migrated outward. This did not occur for Earth since we do not find any Earth analogue with a_mwmi > 2 AU. There are seven of the Mars analogues with a_mwmi > 2 AU in total in the Grand Tack simulations (including all the simulations with equal mass and oligarchic initial condition). We find none of them in the Classical simulations.

We obtained 635 planets in total when the tack location was at 1.5 AU and 685 planets when the tack location was at 2 AU, which greatly exceed the results from the oligarchic initial conditions (see Table 3). The main reason for this is that we have more simulations per set when using the equal-mass initial conditions. We only performed 64 simulations per set for the oligarchic initial conditions, whereas we had 144 simulations for the equal-mass initial conditions. The oligarchic initial conditions, however, also yield fewer planets per simulation. This is because the initial number of embryos in the oligarchic case is usually fewer than the equal-mass case (comparing the numbers between Tables 1 and 2). We would expect that putting more embryos at the beginning of the simulation would help to form more planets at the end. Yet, this is not the case when we consider the formation of the Mars analogues. Instead, we find that increasing the number of embryos does not help to produce more Mars analogues. Figure 5 shows the total number of Mars analogues in each subset of simulations against the initial number of embryos per simulation (see Table 1 for the initial number of embryo adopted for each subset with equal-mass embryos). We have 9 subsets and therefore 9 different initial numbers of embryos: from 27 to 211 embryos. Each subset consists of 32 simulations, in which 16 are simulations with Jupiter’s tack location at 1.5 AU and the other 16 are simulations with Jupiter’s tack location at 2 AU. Even if we start the simulations with more than 200 embryos, the total number of Mars analogues generated were fewer than the simulations with only about 50 embryos. Therefore, there is no correlation between the initial number of embryos and the final number of Mars analogues, which violates our earlier prediction that more embryos at the beginning of simulations would help to form more Mars analogues.

Fig. 4 As in Fig. 3, but for the equal-mass initial conditions. This plot is similar to Figs. 2 and 8 in Brasser et al. (2017), but we only show the Venus (green), Earth (blue), and Mars (red) analogues.

4.4 Summary

Here we examined the possibility of forming Earth and Mars with different accretion zones and therefore different bulk compositions in the Classical and the Grand Tack models. We find that simulations in the Classical model fail to form any of the Mars analogues with an accretion zone distinct from the Earth analogues (a_mwmi > 2 AU) when we only consider Mars analogues close to its current mass. Therefore, the Classical model has trouble reproducing the compositional difference between Earth and Mars. It may not be valid, however, to make such a conclusion since we only have eight Mars analogues in the Classical model due to fewer numbers of simulations and its own drawback in explaining the small mass of Mars. Therefore, we carry out another study comparing the bulk composition of Earth and Mars, presented in the following section (see Sect. 5). The Grand Tack model has a slightly higher success rate in producing compositionally different Earth and Mars because we obtain seven Mars analogues with accretion zones that are wholly distinct from that of the Earth analogues, although statistically the Grand Tack model has a much higher probability to form Earth and Mars with a similar accretion zone. Both the Classical and the Grand Tack models havedifficulty forming Earth and Mars with quite different accretion zones that then lead to different bulk compositions.

5 Chondrite fractions to make Earth and Mars

We now calculate the bulk composition of Earth and Mars from the output of our dynamical simulations. We assume that the solid disk was initially composed of only enstatite and ordinary chondrite. As mentioned in Sect. 2.2, it has been suggested – but not universally accepted – that the solid disk originally had a heliocentric composition gradient consisting of dry, reduced enstatite chondrite close to the Sun, moderately volatile-rich and oxidised ordinary chondrite in the region of the asteroid belt, and highly oxidised and volatile-rich carbonaceous chondrites beyond (Morbidelli et al. 2012; Rubie et al. 2015; Fischer-Gödde & Kleine 2017), mirroring the heliocentric distribution of asteroid groups (Gradie & Tedesco 1982; DeMeo & Carry 2014). Due to their different water fractions, a plausible reason for the transition from enstatite chondrite to ordinary chondrite is snow line migration. However, if this were true then we would also expect the ordinary chondrites to have as much water as carbonaceous chondrites. This is not the case. In truth, we are not aware of any physical explanation for this proposed heliocentric compositional gradient, and we only make use of previously published results.

Although recent studies provide an estimate for the initial location boundary between these two categories of chondrite (Morbidelli et al. 2012; Fischer-Gödde & Kleine 2017; O’Brien et al. 2018), there is still debate on the exact location for the disk to change its composition from enstatite to ordinary. We define this transition’s location as the break location of the disk. Asteroid 21 Lutetia, which is considered to be a candidate source of the enstatite chondrites (Vernazza et al. 2009; Coradini et al. 2011) is currently located at ~2.4 AU. However, it is possible – and consistent with our Grand Tack simulations – that it formed in the inner solar-system and was subsequently scattered by either emerging protoplanets or by the migration of Jupiter (Vernazza et al. 2011). Therefore, the break location of the disk may be closer to the Sun.

We therefore set this break location as a free parameter and assume that it ranges from 1 to 2 AU, with 0.1 AU intervals (i.e. if the break location of the disk is at 1.5 AU, all the planetesimals and embryos initially located within 1.5 AU are made of enstatite chondrite, and those initially located further than 1.5 AU are made of ordinary chondrite). We then make use of the simulation results in Sect. 4 to calculate the average percentage of enstatite and ordinary chondrite incorporated into the Earth and Mars analogues for each break location.

Figure 6 shows the average percentage of enstatite chondrite (red) and ordinary chondrite (blue) incorporated into the Earth and Mars analogues with different break locations of the disk for both models. We only examine data from the simulations with the equal-mass initial condition for the Grand Tack model because the oligarchic initial condition provides similar results in the accretion zone analysis (see Sect. 4). We relax the constraints on the Mars analogues by not considering their masses for the Classical model in this study to ensure that enough data are generated. The total percentage of enstatite and ordinary chondrite at each break location sums up to 100%. With a more distant break location, both Earth and Mars consist mainly of enstatite chondrite since the disk consists of more enstatite chondrite. The break location, however, cannot exceed 1.8 AU, otherwise the bulk composition of both Earth and Mars would consist of about 90% enstatite chondrite, which would make them too similar to each other and would also be inconsistent with Mars’ suggested bulk composition (Sanloup et al. 1999). Therefore, the key here is whether or not a particular break location exists at which we can reproduce the measured isotopic composition for both Earth and Mars simultaneously. Brasser et al. (2018) reported a best-fit composition of enstatite chondrite plus ordinary chondrite for Mars, whereas Mars is suggested to be 45% enstatite chondrite and 55% ordinary chondrite according to Sanloup et al. (1999) and Tang & Dauphas (2014). The break location in between 1.1 and 1.5 AU for both models can satisfy the average bulk composition of Mars in Fig. 6. Combining this with the most recent results for Earth (~70% enstatite chondrite and ~25% ordinary chondrite, Dauphas 2017) leads us to conclude that the break location of the disk should be at 1.3 to 1.4 AU in both the Classical and the Grand Tack models.

It is important to note, however, that even if the average mean composition of Mars analogues is close to 50:50 when the break location is at 1.3 AU, this result does not simply imply that most of the Mars analogues have such a composition. The error bars of Mars in Fig. 6 represent deviations of about ±20% to 30% from the mean values, which suggest that we are dealing with averaging over two extreme cases that happen to meet in the middle. To understand the origin of the large scatter, we plot the cumulative distribution functions (CDFs) of the fraction of enstatite chondrite in the Mars analogues and determine the fraction of them that fall close to the 50:50 mark. Figure 7 shows the result of this analysis, where the break location is at 1.3 AU. The left panel shows Mars analogues from the Grand Tack model with equal-mass initial conditions, and the right panel plots the CDF for Mars analogues from the Classical model with the oligarchic initial conditions. According to Fig. 6, the average mean composition of Mars in both models is close to 50:50. In the Grand Tack, however, only ~20% of them fall in the range wherein they possess 29% to 68% enstatite chondrite by mass (Brasser et al. 2018), whereas ~65% of them do so in the Classical model. Therefore, the Grand Tack model does not work better than the Classical model in producing Mars’ documented composition if we neglect the massive Mars problem of the Classical model.

To better understand the origin of Mars’ roughly 50:50 composition, we tracked its accretion history. By randomly picking 6 Marsanalogues from the Grand Tack simulations and 5 from the Classical simulations, we find that in order to reproduce Mars’ suggested composition, 9 out of 11 involve collision between two embryos of different initial composition. Figure 8 shows an example of the mass and semi-major axis evolution of a Mars analogue from the Grand Tack simulation. The squares represent the initial semi-major axis of the planetesimals and embryos and the time they collided with the Mars analogue. Their sizes are scaled by their masses. In this case, the Mars analogue starts at ~1.5 AU and is therefore an embryo composed of ordinary chondrite because the break location is at 1.3 AU. It stays near its initial location throughout its whole evolution. Its mass contributed by planetesimals made of enstatite chondrite is comparable to those made of ordinary chondrite: a total of seven of the planetesimals originated from the region further than the break location (squares above the horizontal dotted line at 1.3 AU) and five of the planetesimals started in the inner disk (squares below the horizontal dotted line at 1.3 AU), but they are not enough to increase its mass to the current value. Instead, the major reason for its 50:50 composition is a collision in the simulation at 5.7 Myr with another embryo whose composition is enstatite chondrite. This collision, doubling the mass of the Mars analogue is represented by a large square at ~1 AU in the lower panel of Fig. 8. This particular case agrees with the prediction by Brasser et al. (2017) that Mars formed in a distant region and is therefore initially mainly composed of ordinary chondrite. As pointed out in that work, it was then scattered into the inner disk by Jupiter’s migration and possibly collided with an embryo mostly made of enstatite chondrite from the inner disk. This scenario naturally explains the isotopic compositional difference between Earth and Mars. Although collisions in N-body simulations occur at random, a colossal impact at 5.7 Myr is entirely consistent with Mars’ formation timescale.

Fig. 5 Total number of Mars analogues in each subset of simulations as a function of the initial number of embryos per simulation (see Table 1) of the Grand Tack model with equal-mass initial conditions. The results from simulations with Jupiter’s tack location at 1.5 and 2 AU are combined.

Fig. 6 Average percentage contribution of enstatite chondrite (red) and ordinary chondrite (blue) to the bulk compositions of Earth and Mars. Data are plotted as a function of the distance where the initial solid disk changes its composition from enstatite chondrite to ordinary chondrite, which is defined as the break location of the disk. The error bars are ± 1σ. The two upper panels depict the average composition of Earth and Mars in the Classical model simulations. The two lower panelsshow the results from Grand Tack simulations with equal-mass initial conditions, which is similar to Fig. 7 of Brasser et al. (2017).

Fig. 7 Cumulative distribution function of the Mars analogues’ final composition (percentage of enstatite chondrite incorporated into each planet) when the break location of the disk is 1.3 AU. Left panel: results from the Grand Tack simulations with initial equal-mass embryos. Right panel: results from the Classical simulations with initial oligarchic embryos.

6 Using dynamics to produce the isotopic composition of Mars

Although the Grand Tack model is not a better model to explain the composition of Mars, its success in explaining the mass-distance distribution of the inner solar system should not be overlooked. Besides, if Mars formed in a distant region, it had to rely on the migration of Jupiter to scatter to its current location.

In the previous section, we computed the bulk composition of Mars from a mixture of both enstatite and ordinary chondrite. Even though this mixture agrees with published results of Brasser et al. (2018), such compositions are extrapolated from the measurements of isotopic anomalies in Mars’ mantle. Therefore a comparison between our computed composition of Mars from the dynamical simulations with available martian isotopic anomalies is warranted. Here we are interested not only in the predicted average values of the isotopic composition of Mars that result from the dynamical simulations, but also their uncertainties, and, furthermore, whether these uncertainties are comparable to the measured uncertainties, which may allow for additional insight regarding Mars’ origin. Therefore, we proceed to calculate Mars’ expected isotopic composition as follows.

The isotopic anomaly is here defined as the deviation of the isotopic composition of a nuclide with respect to the terrestrial standard value. To obtain the correct isotopic composition for both Earth and Mars, the break location of the disk should be around 1.3 AU (see Fig. 6). We did not include carbonaceous chondrite in Sect. 5 because it only takes up a tiny portion of Earth and Mars according to a few isotopic studies (Sanloup et al. 1999; Warren 2011; Dauphas et al. 2014; Dauphas 2017). When calculating the isotopic anomaly, however, the tiny contribution of carbonaceous chondrite may be enough in some cases to cause significant alteration in the value. To address this, we re-estimate the composition of Mars’ analogues in this section by assuming the solids further than 3 AU are made of carbonaceous chondrite (CO and CV). We multiply the mass of each planetesimal further than 3 AU by ten since we have a relatively low-mass outer disk in the simulations. The isotopic composition of an element in Mars’ mantle ε_Mars can be calculated by (4)

where f_i is the percentage of chondrite type i in Mars (i is either enstatite, ordinary, or carbonaceous chondrite), c_i is the concentration of the element in chondrite type i and ε_i is the isotopic composition of the element in chondrite type i (Dauphas 2017). The mean composition (±1σ) of Mars is enstatite chondrite, ordinary chondrite and carbonaceous chondrite if the break location of the disk is at 1.3 AU (see Fig. 6). The uncertainties in the final composition are the result of the great potential for variety in Mars’ dynamical history, and in the ε_i for each chondritic source. Therefore, a Monte Carlo method is adopted to calculate ε_Mars for each element. We pick the f_i from our dynamcial simulations and ε_i from measured values listed in Dauphas (2017). We use the Box–Mueller transform to approximate a standard normal distribution and then calculate the corresponding ε_Mars from this distribution. We repeat this 60 000 times and obtain a list of values for ε_Mars of an element. The mean and the standard deviation of ε_Mars of each elements is then calculated with this list of values. Our work is different from Brasser et al. (2018) since we adopt the composition of Mars from our dynamical simulation, whereas the results of Brasser et al. (2018) are based on a Monte Carlo Mixing model.

Table 4 shows the average ε_Mars and the 2σ uncertainty of each element. We compare the values that we calculated and those measured in different studies and summarised by Dauphas (2017). The measured value for μ¹⁴²Nd is the mean value of enriched shergottites in Kruijer et al. (2017). We adopt the same computed value as Brasser et al. (2018) for ε⁶⁴ Ni, where they rely on existing ε⁶² Ni isotope data for Mars and the correlation between ε⁶⁴Ni and ε⁶²Ni (Tang & Dauphas 2014). The ε⁶⁴Ni of the martian mantle has not been measured precisely. We observe that all the elements fit within the 2σ uncertainties of the measured values. For most elements, the computed isotopic uncertainties resulting from dynamics are comparable to the measured values. Oxygen, Titanium, and Neodymium are lithophile elements and therefore would remain in the mantle since they were accreted by Earth or Mars (Dauphas 2017). Chromium is moderately siderophile (e.g. Righter & Chabot 2011). A portion of accreted Cr in the mantle sank into the core during the formation of Mars. This may alter the ε_Mars that we calculated from the full accretion of Mars. Nevertheless, the portion of Cr that sank into the core during Earth’s formation istiny (Dauphas 2017). For simplicity, we treat Cr behaviour in Mars in a similar manner to that in Earth.

We adopt a different method in calculating the ε_Mars for ⁶⁴Ni and ⁹²Mo. Both Ni and Mo are moderately siderophile elements and behave more like siderophiles than Cr. It is probable that a significant portion of accreted Ni and Mo sank into the core throughout the accretion of Mars. Dauphas (2017) showed that nearly all the Ni and Mo accreted in the first stage of Earth’s accretion should have sunk into the core. The first stage of Earth’s accretion is defined as the first 60% of Earth’s accretionin terms of mass because previous models suggested that material accreted by Earth changed from more reduced to more oxidised at about 60% accretion (Rubie et al. 2015). It is expected that Ni and Mo should behave similarly in Mars. For Ni and Mo, we ignore all accretion during the first ~60% of Mars’ accretion. This requires a higher resolution of the late-stage accretion data for each Mars analogue, which our simulations lack. Therefore, we can only pick the data from simulations with smaller initial embryo mass (0.025 M_⊕). The average percentage of enstatite, ordinary chondrite, and carbonaceous chondrite is calculated for 20 Mars analogues from their last ~40% of accretion. The ε_Mars value of ⁶⁴Ni and ⁹²Mo are calculated with the same Monte Carlo method described above, but with the f_i obtained from the last ~40% of accretion. The mean bulk composition (±1σ) from N-body simulations is enstatite chondrite, ordinary chondrite, and carbonaceous chondrite for the last ~40% of Mars’ accretion when the break location of the disk is at 1.3 AU. Table 4 shows that the calculated values match with the measured values for both Ni and Mo.

To summarise, it is possible to reproduce the isotopic composition of the martian mantle for six isotopes from our simulation with the Grand Tack model. The uncertainties of the martian isotopic composition calculated from the N-body simulations for 4 out of 6 elements are of the same order of magnitude as the measured data. This is caused mainly by the great variety of Mars’ dynamical evolution in our N-body simulations. Therefore, to better constrain the outcome, more high-resolution simulations are needed. We encourage improved measurements of the isotopic composition of martian meteorites to lower the measured uncertainties.

Fig. 8 Evolution of mass and semi-major axis, respectively, of a Mars analogue (red solid lines in both panels) with a final 50:50 composition in the Grand Tack simulation with initial equal mass embryos. The squares in the lower panel represent the initial semi-major axes and the collision time with the Mars analogue of the planetesimals (smaller squares) and embryos (larger square). The black dotted horizontal line at 1.3 AU in the lower panel is the break location of the initial disk.

7 Conclusions

The isotopic anomaly between Earth and Mars in several nuclides suggests that the bulk compositions of Earth and Mars are different. Assuming they both accreted mainly chondritic material, Earth should have inherited mainly enstatite chondrite (~70%) and Mars is likely to have accreted more ordinary chondrite (32–67%) than Earth. This provides further support to the conclusion that the region where Mars accreted most of its mass was different from that of Earth, and was likely to have been more distant. We first examine the possibility of forming Earth and Mars with different compositions in two different terrestrial planet formation models: the Grand Tack model and the Classical model. By tracing the initial positions of the planetesimals and embryos assembled into each planet, an accretion zone comparison can be made between the Earth and Mars analogues after performing a large number of N-body simulations. Our initial conditions are either embryos with equal masses or with masses that were computed from the traditional oligarchic growth of the planetesimals. We find that there are seven Mars analogues with a_mwmi > 2 AU in the Grand Tack simulations. These are Mars analogues that have a very different accretion zone from all the Earth analogues. They only contribute <5% of the totalMars analogues. We discovered no Mars analogues with a_mwmi > 2 AU in the Classical model when we considered the massive Mars problem generally observed in the Classical model’s simulations. Both the Classical and Grand Tack have difficulty in explaining the compositional difference between Earth and Mars. We also find that the final number of Mars analogues we obtain from the simulations is independent of the initial number of embryos, which violates our original expectation that more initial embryos yield more Mars analogues.

We then estimated the average percentage of enstatite chondrite and ordinary chondrite incorporated into Earth and Mars based on the N-body simulation results. At first we assumed that the solid disk is made entirely of enstatite and ordinary chondrite, with enstatite chondrite residing in the inner region and ordinary chondrite in the outer region. The break location of the disk, at which the disk changes from enstatite chondrite to ordinary chondrite, is set as a free parameter from 1 to 2 AU with 0.1 AU intervals. Our results showthat the break location must be close to 1.3 AU for both models to possess mean compositions of Earth and Mars close to the documented ones. At this break location, the Classical model yields more Mars analogues close to Mars’ documented composition than the Grand Tack model if we neglect the problem of excessive mass of Mars in the Classical model. This composition requires an early collision between two embryos of different compositionthat merge together to form Mars, which agrees with the distant Mars formation scenario (Brasser et al. 2017).

In the previous section, we calculated the isotopic composition from our simulation results for four elements that trace the full accretion history of Mars (¹⁷O, ⁵⁰Ti, ⁵⁴Cr, and ¹⁴²Nd) and two elements that trace only the late accretion of Mars (⁶⁴Ni and ⁹²Mo) in the case of the break location of the disk at 1.3 AU, in which Mars is ~50% enstatite chondrite and ~50% ordinary chondrite. The Grand Tack model matches the measured values within their uncertainties for all elements. Nevertheless, we cannot at the present time confine the dynamical pathway of Mars by matching the calculated isotopic composition to the measured values since their uncertainties are comparable to each other. More dynamical simulations on terrestrial planet formation with higher resolution are needed. We further advocate improved measurements of the isotopic composition of martian meteorites in order to solve the mystery of the formation of Mars and the other terrestrial planets.

Acknowledgements

RB is grateful for financial support from JSPS KAKENHI (JP16K17662). RB and SJM acknowledge the Collaborative for Research in Origins (CRiO), which is supported by The John Templeton Foundation – FfAME Origins program: the opinions expressed in this publication are those of the authors, and do not necessarily reflect the views of the John Templeton Foundation. The source codes for the model used in this study are archived at the Earth Life Science Institute of the Tokyo Institute of Technology. The data, input, and output files necessary to reproduce the figures are available from the authors upon request.

References

All Tables

All Figures

Fig. 1 Initial masses of embryos against their initial semi-major axes for the oligarchic initial condition. The red points are embryos in a 0.5 Myr disk, green points are in a 1 Myr disk, blue points are in a 2 Myr disk, and black points are in a 3 Myr disk. Different point shapes represent different separation between adjacent embryos. Squares depict embryos separated by 10 RHill from each other, triangles are embryos with spacing of 7 RHill and rhombuses are embryos with spacing of 5 RHill. The comparison of different spacing between embryos is only plotted for the 0.5 Myr disk.

In the text

	Fig. 2 Mass-weighted mean initial semi-major axis (amwmi) of material accreted into each Venus analogue (green), Earth analogue (blue), and Mars analogue (red) as a function of their final semi-major axis for the oligarchic initial condition of the Classical model. We combine the results of 10 and 5 RHill and differentages of the solid disk (either t = 0.5, 1, 2, or 3 Myr) within this plot.
In the text

	Fig. 3 As in Fig. 2, but for the Grand Tack model with oligarchic initial conditions. The results of Jupiter tacked at 1.5 or 2 AU, as well as different initial spacing between adjacent embryos (either 10, 7, or 5 RHill) and differentages of the solid disk (either t = 0.5, 1, 2, or 3 Myr) are combined.
In the text

	Fig. 4 As in Fig. 3, but for the equal-mass initial conditions. This plot is similar to Figs. 2 and 8 in Brasser et al. (2017), but we only show the Venus (green), Earth (blue), and Mars (red) analogues.
In the text

	Fig. 5 Total number of Mars analogues in each subset of simulations as a function of the initial number of embryos per simulation (see Table 1) of the Grand Tack model with equal-mass initial conditions. The results from simulations with Jupiter’s tack location at 1.5 and 2 AU are combined.
In the text

Fig. 6 Average percentage contribution of enstatite chondrite (red) and ordinary chondrite (blue) to the bulk compositions of Earth and Mars. Data are plotted as a function of the distance where the initial solid disk changes its composition from enstatite chondrite to ordinary chondrite, which is defined as the break location of the disk. The error bars are ± 1σ. The two upper panels depict the average composition of Earth and Mars in the Classical model simulations. The two lower panelsshow the results from Grand Tack simulations with equal-mass initial conditions, which is similar to Fig. 7 of Brasser et al. (2017).

In the text

	Fig. 7 Cumulative distribution function of the Mars analogues’ final composition (percentage of enstatite chondrite incorporated into each planet) when the break location of the disk is 1.3 AU. Left panel: results from the Grand Tack simulations with initial equal-mass embryos. Right panel: results from the Classical simulations with initial oligarchic embryos.
In the text

Fig. 8 Evolution of mass and semi-major axis, respectively, of a Mars analogue (red solid lines in both panels) with a final 50:50 composition in the Grand Tack simulation with initial equal mass embryos. The squares in the lower panel represent the initial semi-major axes and the collision time with the Mars analogue of the planetesimals (smaller squares) and embryos (larger square). The black dotted horizontal line at 1.3 AU in the lower panel is the break location of the initial disk.

In the text