Timing of Martian Core Formation from Models of Hf – W Evolution Coupled with N-body 1 Simulations 2

11 Determining how and when Mars formed has been a long-standing challenge for 12 planetary scientists. The size and orbit of Mars are difficult to reproduce in classical simulations 13 of planetary accretion, and this has inspired models of inner solar system evolution that are tuned 14 to produce Mars-like planets. However, such models are typically not coupled to geochemical 15 constraints. Analyses of Martian meteorites using the extinct hafnium–tungsten (Hf–W) 16 radioisotopic system, which is sensitive to the timing of core formation, have indicated that the 17 Martian core formed within a few million years of the solar system itself. This has been 18 interpreted to suggest that, unlike Earth’s protracted accretion, Mars grew to its modern size very 19 rapidly. These arguments, however, generally rely on simplified growth histories for Mars. Here, 20 we combine realistic accretionary histories from a large number of N-body simulations with 21 calculations of metal–silicate partitioning and Hf–W isotopic evolution during core formation to 22 constrain the range of conditions that could have produced Mars. 23 We find that there is no strong correlation between the final sizes or orbits of simulated 24 Martian analogs and their 182W anomalies, and that it is readily possible to produce Mars-like 25 Hf–W isotopic compositions for a variety of accretionary conditions. The Hf–W signature of 26 Mars is very sensitive to the oxygen fugacity (fO2) of accreted material because the metal– 27 silicate partitioning behavior of W is strongly dependent on redox conditions. The average fO2 of 28 Martian building blocks must fall in the range of 1.10–1.35 log units below the iron–wüstite 29 buffer to produce a Martian mantle with the observed Hf/W ratio. Martian 182W isotopic 30 signatures are more often reproduced if the planet’s building blocks are sulfur-rich and exhibit a 31 high degree of impactor metal equilibration, but the timing of accretion is a more important 32 control. We find that while Mars must have accreted most of its mass within ~5 million years of 33 solar system formation to reproduce the Hf–W isotopic constraints, it may not have finished 34 accreting until >50 million years later. There is a high probability of simultaneously matching 35 the orbit, mass, and Hf–W signature of Mars even in cases of prolonged accretion if giant 36 impactor cores were poorly equilibrated and merged directly with the proto-Martian core. 37


Introduction
The final stage of terrestrial planet formation began when large protoplanets finished consuming their neighboring small bodies and started perturbing each other's orbits (e.g., Chambers, 2004).This "chaotic growth" period was characterized by dynamical stochasticity: the mutual gravitation of numerous bodies introduces an unavoidable element of randomness into calculations of orbital dynamics, with even miniscule changes in starting conditions significantly altering the simulated solar system evolution (e.g., Lissauer, 2007).This stochasticity has made N-body accretionary simulations (which calculate the gravitational interactions between protoplanets) valuable tools in evaluating the range of possible growth histories for the terrestrial planets.For example, N-body simulations very frequently produce planets with Earth-like masses and orbits, so statistically meaningful interpretations can be drawn about Earth's likely accretionary history (e.g., Kenyon & Bromley, 2006;Canup, 2008;Rubie et al., 2015).
Producing Mars-like planets in N-body simulations, on the other hand, has proven more difficult.Early studies using classical dynamical regimes (i.e., with the Jovian planets near their modern orbits) produced planets near the orbit of Mars that were far too massive, creating the socalled "small Mars problem" (e.g., Wetherill, 1991;Chambers, 2001;Raymond et al., 2009).Though later suites of simulations showed that it is possible, though unlikely, to produce Mars with classical dynamics (Fischer & Ciesla, 2014), many proposed solutions to the small Mars problem involve altering the dynamics of the protoplanetary disk.One of the most successful of these approaches, the Grand Tack model, relies on an inward-then-outward migration of Jupiter to truncate the distribution of disk material at ~1.5 AU (Hansen et al., 2009;Walsh et al., 2011).This migration scatters material that originally condensed outside of 1.5 AU into the inner solar system and reduces the mass available in the Martian feeding zone.Grand Tack N-body simulations produce appropriately small Mars analogs (e.g., Walsh et al., 2011;Jacobson & Morbidelli, 2014;O'Brien et al., 2014), but the model has been criticized from a dynamical perspective; the mechanism, extent, and timing of the required giant planet migration is poorly constrained, and therefore must be tuned to reproduce the observed solar system configuration (e.g., Raymond & Morbidelli, 2014).N-body simulations run under these various accretion scenarios can provide important insights into the accretion history of Mars, such as plausible mass evolution histories and provenance, but Martian formation cannot be understood with dynamical simulations alone.
The accretion history of Mars can also be constrained using geochemical data.The extinct hafnium-tungsten (Hf-W) radioisotopic system, which is sensitive to the timing and conditions of core formation, is a common proxy for determining the formation timescales of terrestrial bodies (e.g., Lee & Halliday, 1995;Kleine et al., 2002;Jacobsen, 2005).This sensitivity comes from the differing chemical affinities of the parent and daughter nuclides: Hf is highly lithophile, but 182 Hf decays (with a ~9 Myr half-life) into 182 W, an isotope of moderately siderophile W. Most of a planet's primordial W is sequestered in its core, but any 182 W produced from 182 Hf decay after the end of core-mantle equilibration (or produced earlier but not efficiently partitioned into the core) remains in the mantle, creating anomalous "extra" 182 W in the planet's mantle.Due to its short decay time, 182 Hf went extinct early in solar system history, fossilizing this signature of core formation.
The Martian Hf-W isotopic composition has been determined from the Shergottite-Nakhlite-Chassignite (SNC) meteorites, a small family of achondrites jettisoned from the Martian crust (Treiman et al., 2000).These meteorites have been widely used to interpret the timing of Martian core formation (e.g., Lee & Halliday, 1997;Righter & Shearer, 2003;Jacobsen, 2005;Dauphas & Pourmand, 2011;Krujier et al., 2017;Marchi et al., 2020), but this effort has several limitations.First, our knowledge of the Bulk Silicate Mars (BSM) composition, derived as it is from <70 kg of material, is poor.Studies of BSM (e.g., Morgan & Anders, 1979;Dreibus & Wänke, 1985;Lodders & Fegley, 1997;Bertka & Fei, 1998;Sanloup et al., 1999;Bouvier et al., 2009;Taylor, 2013;Yoshizaki & McDonough, 2020) tend to either not consider trace elements (like Hf and W) or to disagree on their abundances.Furthermore, the SNC meteorites appear to be derived from several distinct mantle sources that formed during the lifetime of 182 Hf (e.g., Foley et al., 2005); since Hf and W are not equally compatible upon mantle melting, the SNCs have inherited a range of Hf-W signatures.While it is doubtful that any meteorites are directly derived from BSM, the Shergottites imply an earlier (and thus less likely to have been overprinted) core formation age than the other SNCs, and evidence from the Sm-Nd system suggests that their Hf-W signature may be representative of BSM (Kleine et al., 2004;Dauphas & Pourmand, 2011;Krujier et al., 2017).
Using the modern Hf-W signature of Mars to date its formation is difficult because both our understanding of BSM and of planetary accretion are incomplete (Nimmo & Kleine, 2007).Studies that assume the Martian mantle evolved undisturbed following a single core formation event (e.g., Jacobsen, 2005), or that Mars grew from perfectly-equilibrated mass added in infinitesimally-small steps (Harper & Jacobsen, 1996;Dauphas & Pourmand, 2011), have concluded that Mars accreted very early, within 5 Myr of solar system formation.These are not realistic depictions of planetary accretion and differentiation, however.During the chaotic end stages of accretion, growth may have occurred in random intervals from impactors with a variety of masses and compositions.Some studies have approached this issue with more sophisticated parameterizations of Martian formation (Marchi et al., 2020;Zhang et al., 2021) or have utilized N-body simulations that can match proposed formation timescales (Morishima et al., 2013;Woo et al., 2021).Morishima et al. (2013) calculated Hf-W evolution of three Mars-like bodies during N-body simulations of oligarchic growth and found that the Martian Hf-W signature can be matched even for very long (>100 Myr) accretion timescales, depending on conditions like the equilibration fraction of impactor material.Here, we examine a much larger number of Mars analogs produced by N-body models of chaotic growth under both classical and Grand Tack dynamics and trace their Hf-W isotopic evolution histories under a variety of accretionary conditions.This approach allows us to determine which narratives can match the observed geochemistry of Mars, thus providing more realistic constraints on the conditions of its formation.

Methods
We examined the outputs of 116 previously published N-body simulations: 16 in the Grand Tack (GT) regime (O'Brien et al., 2014) and 50 each in the classical Eccentric Jupiter and Saturn (EJS) and Circular Jupiter and Saturn (CJS) regimes (Fischer & Ciesla, 2014).The starting state of these simulations approximates the protoplanetary disk at the transition from oligarchic to chaotic accretion, with mass bimodally distributed between several dozen larger planetary embryos (each with mass of order 10 -2 MEarth) and a few thousand smaller planetesimals (each with mass of order 10 -3 -10 -4 MEarth).We identified Martian analog bodies from the final solar system configuration of each simulation (see Section 3.1) and calculated the Hf-W isotopic evolution implied by the accretionary history of each analog.Bodies in each simulation were assigned an initial composition (Supplementary Table S1), oxygen fugacity (fO2), and S content, and each starting body was differentiated into a core and mantle at the time of solar system formation (equated with CAI condensation, 4.567 Ga;MacPherson, 2014).To account for Hf-W evolution between the formation of the solar system and the start of chaotic growth, mantle Hf-W signatures were evolved undisturbed for 2 Myr before the start of the Nbody simulation (a timescale consistent with the oligarchic-chaotic transition of Kenyon & Bromely, 2006).The final 182 W anomalies of most Mars analogs are relatively insensitive to the details of this early accretionary phase (Supplementary Figure S1).
We tracked the isotopic evolution of every initial body that would eventually accrete into a Mars analog.Between impacts, mantle 182 Hf decayed to 182 W, increasing the 182 W anomaly, which is defined as: where � W 182 W 184 � � is the molar ratio of radiogenic 182 W to the stable reference isotope 184 W, and CHUR is the chondritic uniform reservoir, which is assumed to have experienced no core formation and thus to represent a pristine bulk solar system value (Kleine et al., 2009).In general, a larger ε182W implies more 182 Hf decay after core formation, and thus an earlier equilibration time (Jacobsen, 2005).A larger ε182W may also indicate that more of the 182 W produced before and during core formation was left in the mantle due to a low degree of coremantle equilibration (Morishima et al., 2013).The rate of ε182W growth in a differentiated body also depends on the overall Hf/W ratio of the mantle, quantified as: where 180 Hf and 184 W are stable, non-radiogenic reference isotopes.A mantle with a higher f Hf/W has relatively more Hf (including 182 Hf), and thus produces more 182 W per unit time until the system's extinction.
In our model, each impact was accompanied by an episode of metal-silicate equilibration between the entire impactor mantle, a fraction of the impactor core (kcore), and a fraction of the target mantle (kmantle) (Table 1).Equilibration occurred at a fixed temperature (ΔT) above the chondritic mantle liquidus (Andrault et al., 2011) at a constant fraction (Pfrac) of the core-mantle boundary (CMB) pressure at the time of the impact, which was scaled proportionally to the combined target+impactor mass (assuming a final Martian CMB at 20 GPa; Rivoldini et al., 2011).The fO2 of each equilibration reaction was a mass-weighted average of the fO2 of all the material participating in the equilibration, defined relative to the iron-wüstite (IW) buffer as: where    is the activity of component i in phase j and    is the corresponding mole fraction.This allowed us to use the fO2 of equilibration to calculate the corresponding Fe partition coefficient (DFe): Note that this approach is not self-consistent regarding the number of O atoms present in each protoplanet, but this is a negligible effect due to the lithophile character of O in Mars-sized bodies (e.g., Rubie et al., 2004;Steenstra & van Westrenen, 2018;Brennan et al., 2020).Ni was partitioned identically to Fe ( Ni =  Fe ), S was approximated as perfectly siderophile, all other major elements (plus Hf) were assumed to be perfectly lithophile, and W partitioned between the core and mantle with its partition coefficient calculated as: where a, b, c, d are constants derived from metal-silicate partitioning experiments (Siebert et al., 2011), P is the equilibration pressure, T is the equilibration temperature,   is the number of non-bridging oxygen atoms per silicate tetrahedron (a proxy for the degree of silicate melt polymerization, fixed here at 2.5), and  W is the activity coefficient of W in the metallic phase, calculated after Ma (2001) with activity parameters taken from steelmaking literature (Japan Society for the Promotion of Science, 1988) and considering W-W and W-S interactions.
After calculating the partitioning of Fe and W in a core formation episode, the compositions of the post-impact core and mantle were updated.All isotopes of W partition identically, so for kcore > 0 the equilibrating material contains a portion of 182 W-depleted impactor core.This results in a post-impact mantle closer to the CHUR isotopic ratio, and thus a decrease in ε182W proportional to kcore and the impactor mass.W becomes less siderophile at higher P-T (e.g., Siebert et al., 2011), so re-equilibration of impactor core material in a larger body extracts some W into the mantle.This means that in contrast to ε182W, f Hf/W tends to increase with every impact, though this effect is small since the change in W partitioning is modest over the range of conditions in a Mars-mass planet (e.g., Cottrell et al., 2009).After the final impact, all remaining 182 Hf in the mantle of the fully-grown Mars analog was converted to 182 W, allowing us to compare the implied modern Hf-W signature (i.e., the final ε182W and f Hf/W ) of the analog to that of Mars itself.
1-5 wt% 1.6-3.5 wt% Table 1.Model parameters and ranges tested."Complete range" is the total range investigated for each parameter (the first three parameters must fall between 0 and 1 by definition)."Restricted range" is a more realistic subset of parameter space which we used to match the Hf-W signature of Mars.See Section 4.1 for more details.

Mars analog criteria
We define a Mars analog as a body that survives until the end of a simulation with a semi-major axis of 1-3 AU and mass of <0.2 MEarth (Figure 1).This is a more inclusive definition than is typical for analyses of N-body outputs (e.g., Fischer & Ciesla, 2014;Rubie et al., 2015;Zube et al., 2019), but we are interested in sampling the broadest range of possible accretionary histories.The lowest-mass survivors of each simulation are stranded embryos that accreted only a few planetesimals and therefore remained at approximately their initial masses: ~0.05 MEarth in the EJS/CJS simulations (Fischer & Ciesla, 2014) and ~0.03 MEarth in the GT simulations (O'Brien et al., 2014).These smallest analogs are 2-3× less massive than Mars, but larger than any non-planet in our solar system (for comparison, Ceres has a mass of ~0.0002 MEarth).Planetary mass influences Hf-W isotopic evolution primarily because DW decreases with equilibration depth, but this effect is small over the size range of the analogs, making them viable candidates for investigating possible timescales of accretion for Mars.We find that the properties of the Mars analogs (mass, orbital eccentricity, mass-weighted provenance, accretion time) are uncorrelated with their final semi-major axes in these simulations (Figure 2), implying that dynamical scattering is strong enough that any small planetary body could have ended up in a Mars-like orbit.Furthermore, bodies with Mars-like orbits do not necessarily resemble each other, or Mars, in any other way (Section 4.1).

Figure 1.
Orbital parameters of all Mars analogs (full descriptions in Supplementary Table S2).Symbol size is proportional to body mass.Actual solar system bodies shown for context (Ceres' mass is increased 10× for visibility).We consider the cluster of twelve analogs near Mars (indicated with stars) to have the most "Mars-like" orbits.(Fischer & Ciesla, 2014;O'Brien et al., 2014).For each distribution, the box shows the interquartile range, the line within the box is the median, whiskers extend to ±1.5× the interquartile range, and any outlier analogs beyond that range are shown as open circles."Orbit cutoff" indicates the maximum semi-major axis allowed for Mars analogs (i.e., an orbit cutoff of 2.5 includes all bodies with masses of <0.2 MEarth and semi-major axes of 1-2.5 AU as Mars analogs).Numbers in the top panels indicate the number of analogs found using each orbit cutoff.Dashed red horizontal lines indicate observed Martian values."Provenance" is the massweighted semi-major axis of an analog's building blocks.Note that most analogs start the simulation at >50% of their final mass and that no parameters appear to vary significantly with the orbit cutoff used for any of these accretion scenarios.

Analysis of simulation suites
Evolution of mass and ε182W for all Mars analogs are shown in Figure 3. Prolonged accretion is common in all the simulations; few analogs match the Chambers (2006) homogenous, exponential growth function as parameterized by Dauphas & Pourmand (2011).As those studies pointed out, a planet with a Mars-like f Hf/W must accrete rapidly to match the ε182W of Mars since late equilibration of impactor cores will reduce ε182W.Marchi et al. (2020) found that accretion timescales of up to 15 Myr can be consistent with Mars under certain conditions, but many Mars analogs in these N-body simulations form even more slowly.Most analogs start the simulation at >50% of their final mass, but accretion within the simulation more strongly controls the final Hf-W signature of most analogs (Supplementary Figure S1).This is consistent with the results of Morishima et al. (2013), which found that the contribution of the oligarchic growth period to the final ε182W of Mars was small if accretion continued afterwards.Regardless, Mars analogs in these N-body simulations tend to accrete most of their mass within the brief formation timescales deduced by earlier studies (e.g., 3.3 Myr: Jacobsen, 2005; 1.8 Myr: Dauphas & Pourmand, 2011;2.4Myr: Kleine & Walker, 2017;4.1 Myr: Kruijer et al., 2017), even if they do not reach their final masses until much later.Mars analogs with the longest 50% accretion timescales are also those with the lowest final ε182W.
As expected, GT simulations, in which Jupiter's migration scatters protoplanetary mass towards the Sun, tend to produce many Mars analogs (~1.3 per simulation), and those tend to have more Mars-like orbits than analogs formed in EJS or CJS simulations (Figure 1).GT analogs also take longer to reach their final mass (median 95% accretion times: 9.6 Myr for EJS, 24 Myr for CJS, 45 Myr for GT; Figure 2).EJS and CJS simulations both produce many analogs that orbit further from the sun and with greater eccentricities than Mars, though EJS produces analogs almost twice as often (~1 per simulation versus ~0.6 for CJS), implying that mass ends up either lost or concentrated in a few large bodies under CJS dynamics.CJS is also the only suite to show even a slight possible trend between a body's final semi-major axis and the provenance of its material (Figure 2), consistent with the lower degree of radial mixing in classical dynamical regimes (e.g., Fischer et al., 2018).Dauphas & Pourmand (2011). 182W anomalies are dependent on various accretionary parameters; these ε182W evolution curves were calculated at "reference case" conditions (Pfrac = 0.6, kcore = 0.85, kmantle = 0.4, S = 3.5 wt%, ΔT = 0, fO2 = IW-1.22;see Section 3.3).The shaded grey bar indicates the observed ε182W of Mars (Kruijer et al., 2017; Supplementary Table S3).Each impact is associated with a drop in ε182W proportional to kcore and the target-to-impactor mass ratio.

Parameters influencing Hf-W isotopic evolution
One way to visualize the influence of various model parameters on the resulting Hf-W isotopic signatures is to define a reference set of parameters, then isolate the effect of each parameter by varying them one at a time (Figure 4).An analog's oxidation state is the single most important factor in determining its 182 W signature.The fO2 of equilibration determines the fraction of W that is sequestered in the core, thus setting the mantle f Hf/W and the rate of ε182W increase.Since W partitioning is particularly redox-sensitive (e.g., Cottrell et al., 2009), the mantle of a reduced planet contains much less of the total W than that of a more oxidized planet, resulting in a larger f Hf/W and ultimately a larger final ε182W.The same effect (though smaller in magnitude) can be seen by decreasing Pfrac, increasing ΔT, or increasing bulk S content.These parameters increase both f Hf/W and ε182W approximately equally, but the degree of impactor core equilibration has a different effect: a lower value for kcore results in a slightly lower f Hf/W but a significantly larger ε182W.A low degree of accreted metal re-equilibration causes the final body to inherit more of the signature of its building blocks (W is modestly more siderophile at shallower depths, increasing f Hf/W ), but it also reduces the drawdown of radiogenic 182 W in each impact, allowing ε182W to reach higher values.The effect is qualitatively similar for kmantle but is only significant at very low degrees of mantle equilibration, as is the case for the Earth (Fischer and Nimmo, 2018); in Figure 4, the "reference case" point (kmantle = 0.4) is nearly indistinguishable from the maximum mantle equilibration point (kmantle = 1).We also considered the possibility that kcore may have been lower in giant (embryo-embryo) impacts since hydrodynamic experiments have indicated that direct core merging is likely in these cases (Deguen et al., 2014).Decreasing kcore for giant impacts removes their otherwise irreversiblylarge ε182W reductions (i.e., the large vertical lines Figure 3b), allowing some analogs that experienced embryo-embryo impacts to reach Martian ε182W values.An example of this effect can be seen in Figure 5.  1).The shaded grey region indicates the uncertainty range of measured Martian Shergotty-source values (Supplementary Table S3).The "reference case" model parameters are Pfrac = 0.6, kcore = 0.85, kmantle = 0.4, S = 3.5 wt%, ΔT = 0, and fO2 = IW-1.22,representing a close match between the median of the analogs and Mars (Section 4.1).Other points were calculated with these same values except for the single parameter being varied, which was changed to the value indicated next to each point.Symbols denote the median of all analogs and error bars indicate interquartile ranges.Trends between symbols were calculated by a degree 2 polynomial fit to 10 points evenly spanning each parameter range (not shown).

Reproducing Mars
Considering the similar effects of several model parameters (Figure 4) and the wide distribution of analog properties, it is easily possible to get a few analogs to match the Martian Hf-W signature for many parameter combinations.However, there is a relatively restricted subset of parameter space that is both geophysically and geochemically plausible and results in a significant fraction of analogs matching the observed signature of Mars.As noted above, analog ε182W is most sensitive to initial fO2.There is a limited range (approximately IW-1.35 to IW-1.10) in which any analogs match Mars, regardless of other parameter values; even modestly more reducing conditions result in a wide distribution of analog properties that extends to high f Hf/W and ε182W, overshooting Mars.Fortunately, the average fO2 of core formation on Mars is constrained by its mantle FeO content (i.e., Equation 3), and previous studies have shown that the FeO-derived fO2 of Mars agrees with the permissible range found here (e.g., Righter & Drake, 1996;Rai & van Westrenen, 2013;Rubie et al., 2015;Brennan et al., 2020).Brennan et al. (2020) also used mantle trace elements to constrain the values of Pfrac (0.4-0.6), kcore (0.84-1.0), and kmantle (0.4-1.0).These conditions (high degree of equilibration, intermediate equilibration depth) broadly agree with other investigations of Martian core formation (Kleine et al., 2004;Righter & Chabot, 2011;Yang et al., 2015;Zube et al., 2019), so we restrict our further exploration of parameter space to these ranges.The bulk inventory of volatile elements (especially S) in Mars is controversial, with some studies (e.g., Wang & Becker, 2017;Yoshizaki & McDonough, 2020) favoring much lower abundances than others (e.g., Sanloup et al., 1999;Khan & Connolly, 2008;Taylor, 2013;Steenstra & van Westrenen, 2018).Furthermore, S content cannot be constrained by the Hf-W signature because the effect of changing S is indistinguishable from that of other parameters (Figure 4).We use a maximum bulk S value of 3.5 wt% (35% less S than CI chondrites; Palme & O'Neill, 2014).This corresponds to ~18 wt% S in the Martian core, within the preferred range of most S-rich models and close to the value interpreted from the first seismic measurements of the Martian core (Stähler et al., 2021).While Martian differentiation could have been unusually hot due to 26 Al heating (Sahijpal & Bhatia, 2015), this effect is not well constrained and would vary depending on the time and size of each analog's impacts, so we do not impose a ΔT above the mantle liquidus.
Regardless of the precise parameter combination, we can draw some general conclusions about the types of analogs that best match Mars.First, the Hf-W signature can be matched by any of the dynamical suites.The f Hf/W of the Martian mantle is quite low; for comparison, Earth's value has been estimated as 12 (Jacobsen, 2005), 14 (Kleine et al., 2009), or 25 (Dauphas et al., 14 2014).Given this low f Hf/W , Mars analogs that end up matching ε182W are those that avoid having their ε182W values reset by significant later accretion.This constraint is, however, not as severe as implied by the parametrized accretion curve of Dauphas & Pourmand (2011).Our reference case, for example, includes a GT analog with a 95% accretion time of 63 Myr that has a Marslike orbit and matches the f Hf/W and ε182W of Mars within uncertainty.If the cores of giant impactors are poorly equilibrated, then even large Mars analogs that experienced giant impacts can match the Hf-W signature of Mars for similarly prolonged (e.g., Marchi et al., 2020) accretionary histories (in Figure 5d, 40% of the matching analogs are ≥0.75MMars).Finally, analogs with the most Mars-like orbits (Figure 1) do not necessarily have Mars-like Hf-W signatures, nor do they cluster together in ε182W-f Hf/W space (Figure 5).Despite this, a few analogs with Mars-like orbits (which ones in particular depend on model parameters) often match Martian values, demonstrating that an analog can simultaneously match the orbit and Hf-W signature of Mars in our model.

Disk conditions
In the preceding analysis, we imposed a uniform Mars-like fO2 and bulk S content on all the initial bodies in each N-body simulation.Under these conditions, since the distribution of accretion times does not vary with orbit (Figure 2, bottom row), final ε182W values and orbits are similarly uncorrelated (Supplementary Figure S3).This does not reflect the reality of the protoplanetary disk.At the start of chaotic accretion, there would have been a relationship between a body's composition and the nebular properties of its indigenous orbit.As evidenced by Earth's relatively low fO2 of core formation (<IW-2.0;Geßmann & Rubie, 2000;Li & Agee, 2001;Chabot et al., 2005) and bulk S content (<1 wt%;McDonough, 2003) compared to that of Mars, higher nebular temperatures close to the Sun probably inhibited the condensation of more volatile species.Previous analyses of N-body simulations have examined this effect by imposing variable fO2 on their initial bodies, such as a highly reduced "enstatite chondrite-like" inner solar system surrounded by a more oxidized "ordinary chondrite-like" region (e.g., Rubie et al., 2015).There may be evidence for such discrete reservoirs of material existing in the early solar system (e.g., Warren, 2011;Morbidelli et al., 2016;Lichtenberg et al., 2021), but their spatial and temporal boundaries, as well as their bulk chemistries, are poorly constrained.
Matching the narrow range of permissible Martian bulk fO2 requires the oxidation state of these reservoirs (and the location of the boundary between them) to be precisely tuned to a particular accretionary provenance (Figure 6).This could be taken as an indication that the bulk fO2 of Mars represents a single reservoir rather than a mixture, but such a simple primordial fO2 distribution is probably unrealistic (e.g., Ciesla & Cuzzi, 2006).The situation becomes even more complicated if the bulk S gradient does not coincide with variations in fO2 or if disk dynamics displace material far from its region of condensation by the time of chaotic growth (e.g., Williams et al., 2020).Therefore, we have chosen to impose a Mars-like composition and redox state on all Mars analogs, allowing us to focus on how their accretionary history influences their Hf-W evolution.It is worth noting that most analogs (including those with Mars-like orbits) have provenances of ≥2 AU (Figure 2), so it is possible to match Martian geochemistry and simultaneously form a reduced Earth from material originating closer to the Sun.S3).The counts of analogs whose f Hf/W plot off-scale are shown in parentheses at the top of the figure.Most analogs are either too oxidized (low f Hf/W , average provenance outside the step) or much too reduced (very high f Hf/W , average provenance inside the step).Provenance is quantified as the mass-weighted semi-major axis of an analog's building blocks and indicated by the color of each symbol.It is difficult to produce Mars by mixing reservoirs; very few analogs accrete exactly the right proportions of material unless one or both reservoirs closely match the fO2 of Mars.GT analogs have a relatively narrow distribution of provenances (Figure 2), so their distribution is more sensitive to radial variations in the disk.

Other N-body approaches
While it is impossible to perfectly simulate the complex physics of planetary accretion, there have been significant advances in N-body techniques since the creation of our simulation suites.For example, Woo et al. (2021) used improved computational power to run an N-body model with more and smaller initial bodies, thus allowing the simulation's start time to closely coincide with that of the solar system.In agreement with our results, that study found that most Mars analogs accrete more slowly than the exponential growth curve of Dauphas & Pourmand (2011) and proposed various dynamical methods to make them grow more quickly.These include the implementation of non-perfect merging between colliding protoplanets, an effect which slightly prolongs Earth's accretion (e.g., Chambers, 2013;Dwyer et al., 2015) but could potentially form Mars more quickly via fragmentation (Kobayashi & Dauphas, 2013;Dugaro et al., 2019).N-body studies of "pebble accretion" disagree on whether that regime promotes (Levison et al., 2015;Matsumura et al., 2017) or discourages (Voelkel et al., 2021) the formation of small terrestrial planets.Broadly, our approach could be extended to any number of N-body simulation types, including ones with different disk dynamics, pre-simulation periods, or impact outcomes, but it seems likely that the Hf-W signature of Mars can be reproduced in a variety of circumstances despite prolonged accretion.One possible exception to this is if long-lived nebular gas postpones the start of chaotic accretion until well after solar system formation (Walsh & Levison, 2019;Chambers et al., 2020).In that case, matching the 182 W anomaly of Mars would require either very low kcore or that Mars experienced essentially no accretion during the chaotic epoch.Simulations with a long-lived nebula have yet to successfully solve the small-Mars problem, but Hf-W evolution under this regime could be a target of future studies.

Conclusions
The Martian Hf-W signature (f Hf/W and ε182W) can be reproduced by modeling W partitioning for successive stages of core formation in N-body accretion simulations.As suggested by some recent studies (e.g., Marchi et al., 2020;Woo et al., 2021;Zhang et al., 2021), we find that many Mars analogs experience substantially protracted accretion, in contrast to the rapid exponential growth of Dauphas & Pourmand (2011).While proto-Mars likely reached 50% of its final size within 5 Myr of solar system formation, it may not have finished growing until >50 Myr later.Exactly which accretionary histories match Mars is dependent on model parameters.Hf-W evolution is particularly sensitive to the oxidation state of metal-silicate equilibration, constraining initial fO2 to a narrow range (IW-1.35 to IW-1.10) consistent with the FeO content of the Martian mantle.This sensitivity means that reproducing Mars by substantial accretion of material from two reservoirs of dramatically differing fO2 is a low-probability event.
While GT dynamics allow analogs to form with Mars-like orbits much more often than EJS or CJS scenarios, we do not find that a Mars-like orbit correlates with Mars-like chemistry, or that GT analogs have a higher probability of matching the Hf-W signature of Mars.Indeed, analogs formed by GT dynamics tend to accrete material from a narrower range of orbits and finish forming later, slightly reducing their ranges of acceptable model parameters.Nonetheless, there are reasonable parameter combinations by which any of the dynamical regimes investigated can match both the orbit and Hf-W signature of Mars simultaneously, even with substantially prolonged accretion.

Figure 2 .
Figure 2. Distributions of orbital and accretionary parameters of Mars analogs from N-body simulations(Fischer & Ciesla, 2014;O'Brien et al., 2014).For each distribution, the box shows the interquartile range, the line within the box is the median, whiskers extend to ±1.5× the interquartile range, and any outlier analogs beyond that range are shown as open circles."Orbit cutoff" indicates the maximum semi-major axis allowed for Mars analogs (i.e., an orbit cutoff of 2.5 includes all bodies with masses of <0.2 MEarth and semi-major axes of 1-2.5 AU as Mars analogs).Numbers in the top panels indicate the number of analogs found using each orbit cutoff.Dashed red horizontal lines indicate observed Martian values."Provenance" is the massweighted semi-major axis of an analog's building blocks.Note that most analogs start the simulation at >50% of their final mass and that no parameters appear to vary significantly with the orbit cutoff used for any of these accretion scenarios.

Figure 4 .
Figure 4. Sensitivity of the Hf-W isotopic signature of Mars to model parameters varying over their "complete ranges" (Table1).The shaded grey region indicates the uncertainty range of measured Martian Shergotty-source values (Supplementary TableS3).The "reference case" model parameters are Pfrac = 0.6, kcore = 0.85, kmantle = 0.4, S = 3.5 wt%, ΔT = 0, and fO2 = IW-1.22,representing a close match between the median of the analogs and Mars (Section 4.1).Other points were calculated with these same values except for the single parameter being varied, which was changed to the value indicated next to each point.Symbols denote the median

Figure 6 .
Figure 6.Mars analog f Hf/W values calculated with an initial step function in fO2 (from IW-4 in the inner disk to IW-1 in the outer disk) imposed on the initial bodies in each N-body simulation.The shaded grey bar shows the Martian f Hf/W value (Kleine & Walker, 2017; Supplementary TableS3).The counts of analogs whose f Hf/W plot off-scale are shown in parentheses at the top of the figure.Most analogs are either too oxidized (low f Hf/W , average provenance outside the step) or much too reduced (very high f Hf/W , average provenance inside the step).Provenance is quantified as the mass-weighted semi-major axis of an analog's building blocks and indicated by the color of each symbol.It is difficult to produce Mars by mixing reservoirs; very few analogs accrete exactly the right proportions of material unless one or both reservoirs closely match the fO2 of Mars.GT analogs have a relatively narrow distribution of provenances (Figure2), so their distribution is more sensitive to radial variations in the disk.