{'fraction_non_alphanumeric': 0.04518108101681212, 'fraction_numerical': 0.028701876203838014, 'mean_word_length': 4.690663058186739, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Timing of Martian Core Formation from Models of Hf–W Evolution Coupled with N-body Simulations

Determining how and when Mars formed has been a long-standing challenge for planetary scientists. The size and orbit of Mars are difficult to reproduce in classical simulations of planetary accretion, and this has inspired models of inner solar system evolution that are tuned to produce Mars-like planets. However, such models are not always coupled to geochemical constraints. Analyses of Martian meteorites using the extinct hafnium–tungsten (Hf–W) radioisotopic system, which is sensitive to the timing of core formation, have indicated that the Martian core formed within a few million years of the start of the solar system itself. This has been interpreted to suggest that, unlike Earth’s protracted accretion, Mars grew to its modern size very rapidly. These arguments, however, generally rely on simplified growth histories for Mars. Here, we combine likely accretionary histories from a large number of N-body simulations with calculations of metal–silicate partitioning and Hf–W isotopic evolution during core formation to constrain the range of conditions that could have produced Mars. We find that there is no strong correlation between the final masses or orbits of simulated Martian analogs and their 182W anomalies, and that it is readily possible to produce Mars-like Hf–W isotopic compositions for a variety of accretionary conditions. The Hf–W signature of Mars is very sensitive to the oxygen fugacity (fO2) of accreted material because the metal– silicate partitioning behavior of W is strongly dependent on redox conditions. The average fO2 of Martian building blocks must fall in the range of 1.3–1.6 log units below the iron–wüstite buffer to produce a Martian mantle with the observed Hf/W ratio. Other geochemical properties (such as sulfur content) also influence Martian 182W signatures, but the timing of accretion is a more important control. We find that while Mars must have accreted most of its mass within ~5 million years of solar system formation to reproduce the Hf–W isotopic constraints, it may have continued growing afterwards for over 50 million years. There is a high probability of simultaneously matching the orbit, mass, and Hf–W signature of Mars even in cases of prolonged accretion if giant impactor cores were poorly equilibrated and merged directly with the proto-Martian core. This manuscript has been accepted for publication in Geochimica et Cosmochimica Acta


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). Later suites of simulations showed that it is possible, though unlikely, to produce Mars with classical dynamics (Fischer & Ciesla, 2014), but many proposed solutions to the small Mars problem involve altering the dynamics of the protoplanetary disk. Such approaches include planetary growth by collisional fragmentation (Clement et al., 2019) or pileup of centimetersized bodies (Drążkowska et al., 2016). One of the most popular alternate dynamical narratives, 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 are poorly constrained, and therefore must be tuned to reproduce the observed solar system configuration (e.g., . 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. Much 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 rocks. 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 ejected 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, 1995;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 less than 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 our understanding of both BSM and 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), though not all simulations form Mars analogs (Section 3.1). 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 sulfur (S) content, and each starting body was differentiated into a core and mantle at the time of solar system formation (equated with CAI condensation at 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 of the differentiated analogs were evolved undisturbed for 2 Myr before the start of the N-body 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: 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 (e.g., Morishima et al., 2013). The rate of ε182W growth in a differentiated body 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). Equilibration occurred at a fixed temperature (ΔT) above or below 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). In real planets, these values were likely correlated; for example, equilibration with ΔT < 0 implies a partially crystallized mantle and thus a low equilibration pressure. However, we allow all parameters to vary independently over reasonable ranges (Table 1) to isolate the effect of each. The fO2 was 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. In most simulation suites, starting bodies were assigned a fixed initial fO2 (Section 3.3) which was constant for all equilibration steps. In cases with non-homogenous initial fO2 (Section 4.2), the fO2 of each equilibration was a mass-weighted average of all bodies that had accreted to the target and impactor. The fO2 of equilibration was used to calculate the corresponding Fe partition coefficient (DFe): Note that this approach does not account for the possibility of other elements dissolving into the core and is therefore not self-consistent regarding the number of O and Si atoms present in each protoplanet, but this is a negligible effect due to the lithophile character of those elements 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 and C were 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, and d are constants derived from metal-silicate partitioning experiments, n is the valence of W, 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), and W is the activity coefficient of W in the metallic phase, calculated after Ma (2001). The reported partitioning behavior of W varies between experimental studies (e.g., Cottrell et al., 2009;Seibert et al., 2011;Wade et al., 2013;Jennings et al., 2021), influencing the values of the constants in Equation 5 (Table S2). In this study, we use the values of Jennings et al. (2021), as it compiles and reanalyzes the results of the earlier studies. We also use that study's chemical activity parameters and consider W-W, W-C, 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. 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. In contrast to ε182W, the change of f Hf/W in each impact 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;Jennings et al., 2021). 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.
IW-1.7 to IW-1.2 IW-1.6 to IW-1.3 bulk S 1-5 wt% 1.6-3.5 wt% Table 1. Model parameters and ranges tested. "Complete range tested" is the total range investigated for each parameter (the first three parameters must fall between 0 and 1 by definition). "Constrained range" is a more realistic subset of parameter space that we used to match the Hf-W signature of Mars within the restrictions of previous geochemical studies. 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). With this definition, 33 of the 50 EJS simulations produce a Mars analog, along with 23 of 50 CJS and 13 of 16 GT simulations. This is a broader definition than is typical for analyses of N-body outputs (e.g., Fischer & Ciesla, 2014;Rubie et al., 2015;Zube et al., 2019), because 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 changes 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 terms of mass or Hf-W signature (Section 4.1). Supplementary Table S3). Symbol size is proportional to body mass. Actual solar system bodies shown for context (Ceres' mass is increased 10× for visibility). All analogs are included in the subsequent analysis; the cluster of twelve analogs with the most "Mars-like" orbits (indicated with stars) are examined in Section 4.1.  (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.4 Myr: 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 also tend to be 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., . Simulations begin 2 Myr after CAI formation. The analytical growth curve (solid black curve) and 95% confidence intervals (dashed black curves) are from Dauphas & Pourmand (2011). 182 W 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.47; see Section 3.3). The shaded bar indicates the observed ε182W of Mars (Kruijer et al., 2017;Supplementary Table S4). Each impact is associated with a drop in ε182W proportional to kcore and the target-to-impactor mass ratio.

Parameters influencing Hf-W evolution
One way to visualize the influence of various model parameters on the resulting Hf-W 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 is high valence (+4 or +6), its partitioning is particularly redox-sensitive (e.g., Cottrell et al., 2009) and 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 increasing ΔT, increasing Pfrac, 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 lower f Hf/W but a larger ε182W. This is due to the action of two competing effects. A low degree of accreted metal re-equilibration causes the final body to inherit more of the signature of its building blocks, and W is less siderophile at shallower depths in the Jennings et al. (2021) parameterization, decreasing f Hf/W . Simultaneously, a small kcore reduces the drawdown of radiogenic 182 W in each impact, allowing ε182W to reach higher values. The effect is qualitatively similar for kmantle but, as is the case for the Earth , is only significant at very low degrees of mantle equilibration; 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.  Table S4). 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.47, 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 4 points evenly spanning each parameter range (not shown). Figure S2 is an alternate version of this figure using the Siebert et al. (2011) W metal-silicate partitioning parameterization.
The presence of C can significantly alter W partitioning (e.g., Jennings et al., 2021), but the C content of the Martian core is thought to be very low due to the absence of isotopic fractionation signatures  and the decrease of C solubility in S-rich alloys (Tsuno et al., 2018). Studies with the most C-rich Martian cores (e.g., Steenstra & van Westrenen, 2018) tend to propose abundances of ~1 wt% (i.e., 0.2 wt % bulk planetary C). We performed one set of calculations with C at this level and found that its presence increased f Hf/W by 0.05 and ε182W by 0.02, a very small effect.

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.3 to IW-1.6) 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. Note that the large kcore value interpreted for Mars can equivalently imply efficient mixing of differentiated impactors or accretion of undifferentiated impactors (which necessarily have kcore = 1), and these two cases are indistinguishable in terms of their effect on Hf-W.
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), temperatures during planetary formation are not well constrained and would vary depending on the time and size of each analog's impacts, so we do not impose a ΔT.
With these restrictions, we define our "reference case" as a set of parameters that produces a close match between the median analog Hf-W signature and Martian values ( Figure  5a): Pfrac = 0.6, kcore = 0.85, kmantle = 0.4, S = 3.5 wt%, ΔT = 0, and initial fO2 = IW-1.47. With these values, 20% of EJS analogs, 31% of CJS analogs, and 14% of GT analogs fall within uncertainty of the Martian Shergottite source. If we instead reduce the S content to produce a core with ~8 wt% S (i.e., 1.6 wt% bulk S, which is 70% less than CI chondrites), the best match to Mars is obtained with slightly more reducing conditions (IW-1.50; Figure 5b). A tighter clustering around Martian f Hf/W is achieved at intermediate S abundance (2.7 wt%; Figure 5c). If most analogs have approximately Martian f Hf/W , matching ε182W depends almost entirely on accretionary history. Keeping all other parameters the same but allowing cores to merge directly (kcore = 0) during giant embryo-embryo impacts (Figure 5d) creates more matching analogs than in our reference case (31% EJS, 28% CJS, 9.5% GT), and some of these matches experience giant impacts tens of Myr into solar system evolution (Figure 6d Table S4). a. Reference case: fO2 = IW-1.47, 3.5 wt% S, kmantle = 0.4, kcore = 0.85, Pfrac = 0.6. b. A low-S case: fO2 = IW-1.50, 1.6 wt% S, kmantle = 0.4, kcore = 0.85, Pfrac = 0.6. c. A good match to f Hf/W only: fO2 = IW-1.5, 2.7 wt% S, kmantle = 0.4, kcore = 1.0, Pfrac = 0.4. d. Same as c but with kcore = 0 (direct core merging) for embryo-embryo impacts. 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., 2014. Given this low f Hf/W , Mars analogs that end up matching ε182W tend to be those that avoid having their ε182W values reset by significant accretion after ~2 182 W half-lives have elapsed (i.e., >20 Myr post-CAI). 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 Mars-like orbit and matches the f Hf/W and ε182W of Mars within uncertainty. If the cores of giant impactors equilibrate with the same kcore as planetesimal impactors, then Mars is unlikely to have experienced an embryo-embryo collision, especially after 20 Myr post-CAI. However, if giant impactor cores are poorly equilibrated, then even large Mars analogs that experienced giant impacts can match the Hf-W signature of Mars for prolonged (e.g., Marchi et al., 2020) accretionary histories (in Figure 5d, 30% of the matching analogs are ≥0.75 MMars). 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 Hf-W values, demonstrating that analogs can simultaneously match the orbit and the core formation signature of Mars in our model. Indeed, any analog can be made to match the Martian Hf-W signature if the accretionary parameters (especially fO2) are tuned appropriately.
It is important to note that the impacts modelled here are not part of a late veneer. Any impactor that mixes into the mantle without participating in core formation (whether because it is fully oxidized or arrives after core formation is complete) will add isotopes in chondritic ratios and therefore reduce f Hf/W and ε182W proportional to the impactor-to-mantle ratio (since the chondritic values are 0 for each parameter). The mass of the late veneer on Mars was constrained by Dale et al. (2012), which showed that Martian highly-siderophile element (HSE) abundances are approximately 2.5 orders of magnitude lower than those of CI chondrites. Assuming no HSEs remained in the mantle after core formation, the late veneer cannot account for more than 0.3% the mass of the mantle, or ~1.5 × 10 21 kg for a Mars-sized analog. This is four times smaller than the planetesimals in our N-body simulations, which implies that the ε182W reduction from late veneer would minimal, about 3 × 10 -3 ε units. If Mars did, in fact, accrete mass for tens of Myr, all those impactors participated in core formation, implying that the Martian magma ocean was either long-lived (e.g., Debaille et al., 2007) or it was regenerated by each impact.

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 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. As Rubie et al. (2015) concluded based on Martian FeO content, it is difficult to produce Mars by mixing reasonable inner solar system and outer solar system reservoirs (Figure 7). Since Hf-W is so sensitive to fO2, the only "successful" fO2 distribution (i.e., one that allows numerous analogs to match Mars) is one in which a large portion of primordial material closely matches the final fO2 of Mars. In any other case, very few analogs accrete exactly the right proportions of reduced and oxidized material. This could be taken as an indication that the bulk fO2 of Mars represents a single, local reservoir (e.g., Mah & Brasser, 2021) 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.  Table S4). The counts of analogs whose f Hf/W plot off-scale are shown in parentheses at the top of the figure. Unless fO2 distributions are contrived such that much of the disk has approximately Martian fO2 (i.e., distribution d), almost all analogs end up either too oxidized (very low f Hf/W ) or too reduced (very high f Hf/W ). Provenance is quantified as the mass-weighted semi-major axis of an analog's building blocks and indicated by the color of each symbol.

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). Studies disagree on whether "pebble accretion" promotes (Levison et al., 2015;Matsumura et al., 2017) or discourages (Voelkel et al., 2021) the formation of small terrestrial planets, but a pebble contribution speeding up pre-chaotic growth (Lichtenberg et al., 2021) could possibly have helped Mars form within the timeframe of Dauphas & Pourmand (2011). Additionally, since very small bodies will fully equilibrate in the magma ocean (i.e., have kcore = 1; Section 4.1), significant pebble accretion could help explain Mars' high kcore value.
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. As noted above, however, later impactors have larger negative effects on ε182W because of the ongoing decay of 182 W, with anomaly growth dramatically slowing after 20 Myr post-CAI. Simulations in which long-lived nebular gas postpones the start of chaotic accretion (e.g., Walsh & Levison, 2019;Clement et al., 2020) are likely to produce systematically low ε182W if there is late addition of material to Mars ( Figure S3d). A possible workaround could be if late accreting material avoids resetting the mantle anomaly because it does not sequester radiogenic 182 W to the core (kcore = 0). Simulations with long-lived nebulae 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.6 to IW-1.3) 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.
As in previous studies (e.g., Kleine et al., 2004;Righter & Chabot, 2011;Yang et al., 2015;Zube et al., 2019;Brennan et al., 2020), we find that Martian material could have been highly equilibrated, with the caveat that larger, later-accreting analogs best match the Hf-W signature of Mars if giant impactor cores directly merged with that of proto-Mars (e.g., Deguen et al., 2014).
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.