I/Pu reveals Earth mainly accreted from volatile-poor differentiated planetesimals

The observation that mid-ocean ridge basalts had ~3× higher iodine/plutonium ratios (inferred from xenon isotopes) compared to ocean island basalts holds critical insights into Earth’s accretion. Understanding whether this difference stems from core formation alone or heterogeneous accretion is, however, hindered by the unknown geochemical behavior of plutonium during core formation. Here, we use first-principles molecular dynamics to quantify the metal-silicate partition coefficients of iodine and plutonium during core formation and find that both iodine and plutonium partly partition into metal liquid. Using multistage core formation modeling, we show that core formation alone is unlikely to explain the iodine/plutonium difference between mantle reservoirs. Instead, our results reveal a heterogeneous accretion history, whereby predominant accretion of volatile-poor differentiated planetesimals was followed by a secondary phase of accretion of volatile-rich undifferentiated meteorites. This implies that Earth inherited part of its volatiles, including its water, from late accretion of chondrites, with a notable carbonaceous chondrite contribution.


INTRODUCTION
Earth must have accreted from diverse materials, but the nature and temporal sequence of delivery of these potential building blocks remain heavily debated (1)(2)(3)(4)(5)(6). To investigate these questions, the isotopes of xenon (Xe), the heaviest stable noble gas, are particularly useful. Because 129 Xe* comes from radioactive beta decay of now extinct volatile 129 I (t 1/2 = 15.7 Ma) and 136 Xe* Pu comes from spontaneous fission of extinct refractory 244 Pu (t 1/2 = 80 Ma), the 129-Xe*/ 136 Xe* Pu ratio evolves as a function of both time and reservoirs compositions (i.e., I/Pu ratio) early in Earth's history. Hence, the study of the 129 Xe*/ 136 Xe* Pu in silicate reservoirs of Earth has the potential to place strong constraints on Earth's accretion and evolution (7)(8)(9)(10)(11)(12)(13). According to recent high-precision analyses of Xe isotopes, ocean island basalt (OIB) samples [plume mantle sources, originating from as deep as the core mantle boundary (CMB)] display a uniformly low 129 Xe*/ 136 Xe* Pu (by a factor of 2.8) compared to mid-ocean ridge basalt (MORB) samples (upper mantle sources) (7)(8)(9)13). Previous work has shown that these signatures cannot simply result from shallow atmospheric contamination, mixing between subducted atmospheric Xe and MORB Xe, and/or different closure ages of Xe loss between the shallow and deep mantle reservoirs (7,9,13). Instead, the low 129 Xe*/ 136 Xe* Pu in the plume reservoir indicates that a low I/Pu was established before 129 I extinction (i.e., first~80 to 100 Ma of the Solar System) and has been preserved thereafter, avoiding rehomogenization, for about 4.45 billion years (7,10,13).
Because iodine is thought to be retained since the earliest stages of accretion (14,15), two main competing mechanisms have been proposed to explain the I/Pu contrast between MORB and OIB mantle reservoirs: (i) a heterogeneous volatile accretion history for Earth (7,10) or (ii) a homogeneous volatile accretion history where partitioning of iodine into liquid metal during core formation was taken into account (11,16). Both models have important ramifications for our understanding of Earth's evolution. In the first case, volatile elements (including iodine) would be depleted in early accreted materials compared to later building blocks of Earth and inefficiently mixed into Earth's whole mantle. In the second scenario, a change in the nature of Earth's building blocks is not required because the iodine depletion of the deeper mantle could be achieved through episodes of high-pressure core formation. However, both models also suffer from notable shortcomings. The heterogeneous accretion models (7, 10) did not consider the impact of core formation processes, and neither model considered the geochemical behavior of plutonium at the high P-T (pressure-temperature) conditions relevant to core formation. This last parameter could exert a strong control on the evolution of I/Pu ratios within Earth's reservoirs during its accretion. In the near absence of experimental data on metal-silicate Pu partitioning (17)-due to its highly radioactive nature and lack of access to suitable amount to perform these experiments-Pu has been implicitly assumed to be a rock-loving (lithophile) element independent of temperature and pressure throughout the planetary accretion and core formation processes.
To remedy this situation and assess whether core formation alone without a change in volatile content of building blocks could explain the difference in I/Pu ratios between MORBs and OIBs, we used the two-phase first-principles molecular dynamics (FPMD) method (18)(19)(20) to determine the partition coefficient of I and Pu between liquid iron and silicate melt. This method has successfully predicted the partition coefficients of C, He, Mg, Si, and O (18,20). Our newly derived partition coefficients can then be used in multistage core formation model to assess whether I/Pu could be efficiently fractionated under different accretion scenarios. Enhanced by a thorough compilation and analysis of available meteorite xenon isotope data, we present a model for the accretion history of Earth based on I/Pu evidence.

RESULTS
The FPMD simulations were ran at pressures ranging from 25 to 85 GPa along the mantle liquidus (21). These values were chosen to cover the range of P-T conditions relevant to Earth's accretion according to the current prevailing literature (22,23). Each simulation started with 256 Mg, Si, O, and Fe atoms (in proportions reflecting the bulk Earth composition, see Materials and Methods) and 4 I or Pu atoms placed randomly in the simulation boxes. After 30,000 steps of equilibration, the atoms are segregated into two phases, a liquid iron phase and a silicate melt phase. A polyhedron "alpha shape" of Fe cluster is constructed to mark the boundary between the two phases ( fig. S1). Atomic configurations from the next 30,000 to 50,000 simulation steps are used to calculate the average compositions of the liquid iron and silicate melt phases (tables S1 and S2).
Partitioning of plutonium and iodine between silicate melt and metal liquid takes place according to the following reactions (see Materials and Methods) PuO silicate The corresponding equilibrium constants, K Pu and K I , are of the form, ln(K) = a + b/T + cP/T, where T is temperature (in Kelvin) and P is pressure (in gigapascal), and were derived by fitting the exchange coefficients (K D ) of our data and previous laser-heated diamond anvil cell (LH-DAC) experimental data of iodine (11) using least-squares regression, yielding lnðK I Þ ¼ ð1:34 + 1:51Þ À ð13:4 + 4: lnðK Pu Þ ¼ ð11:1 + 1:1Þ À ð56:8 + 3:5Þ 1000 T ð4Þ The first-principles calculation data that we obtained for iodine are fully consistent with previous experimental data (Fig. 1A), further justifying the adequacy of our FPMD method for estimating elemental partition coefficients (18,20). Our results show that, similar to iodine, Pu becomes less lithophile at a higher temperature (Fig. 1B). This temperature dependency of K Pu is in line with the similar behavior of another actinide, uranium, observed in LH-DAC experiments (24,25).

DISCUSSION
Both I and Pu tend to partition into the metal phase under high pressures and temperatures relevant to core formation (Fig. 1). However, given that the I and Pu partition coefficients (D = X metal /X silicate ) have different P-T-X metal O (oxygen molal fraction in metal phase) dependency (figs. S2 and S3), the silicate mantle would lose I and Pu to the core in variable proportions during metal-silicate equilibration events. To assess whether core formation could sufficiently fractionate I and Pu to explain the 129 Xe*/ 136 Xe* Pu offsets between the OIB and the MORB mantle sources, we run a series of multistage core formation simulations (22) (see Materials and Methods). To explain the lower 129 Xe*/ 136-Xe* Pu in the plume mantle, a factor of 2.8 depletion of I/Pu would be required, assuming a synchronous closure time of Xe loss for the whole mantle (9,11). This is a conservative assumption because (i) prolonged Xe loss from the mantle to the atmosphere would tend to decrease the upper mantle 129 Xe*/ 136 Xe* Pu and require an even higher initial I/Pu ratio in the upper mantle to account for the observed contrast between MORB and OIB 129 Xe*/ 136 Xe* Pu ratios (13) and (ii) a later closure time of the plume mantle source compared to the upper mantle would be physically difficult to explain.

Homogeneous volatile accretion
First, we assume a homogeneous volatile accretion scenario, meaning that in each simulation, all building blocks are assumed to have the same I/Pu ratio (i.e., volatile content). Equilibrium pressure, temperature, and oxygen fugacity ( fO 2 ) during accretion were varied to test the sensitivity of the simulation outcomes (see Materials and Methods and figs. S4 to S6). Under the assumption of homogeneous volatile accretion (i.e., constant I/Pu in Earth's building blocks), we find that core formation processes do not result in substantial fractionation of I from Pu in Earth's mantle. I and Pu partition into the metal phase very similarly during core formation processes (e.g., Fig. 2A and fig. S7), and mantle I/Pu ratios barely vary during accretion, with the mantle I/Pu ratio of the fully accreted Earth being only 8% lower than that of the initial protomantle (Fig. 2B). This is the case even when the equilibrium pressure of last accretion is assumed to be the upper bound (80 GPa) needed to generate a distribution of other moderately siderophile elements (e.g., Ni, V, Co, and Cr) consistent with observations (22) (fig. S4). Even assuming an extreme accretion scenario with higher equilibration pressure in the middle stage of accretion compared to the late stages (11), the maximum difference between the accreting and the final mantle I/Pu ratios is only~25% ( fig. S4, B and C). Hence, the isotopic Xe signature in the plume mantle source is unlikely to be readily explained by episodes of a high equilibrium pressure during any stage of a homogeneous accretion.
The above calculations were run assuming accretion starting from reduced building materials, but changes in fO 2 during accretion [simulated here using prescribed evolution curves for the FeO content of the magma ocean (22)] could have a notable effect on the I/Pu ratio of the mantle of proto-Earth ( fig. S5). This effect of changes in fO 2 depends on the assumed T of metal-silicate equilibration. When equilibration occurs along the mantle liquidus geotherm based on the melting of a chondritic mantle (cold liquidus geotherm hereafter) (21), the I/Pu ratio for Earth's mantle monotonically decreases throughout Earth's accretion ( fig. S5A). Regardless of the fO 2 conditions, under the cold liquidus geotherm, isolation of an early formed mantle reservoir would yield higher I/Pu ratio in OIBs than MORBs, the opposite of what is observed. In contrast, at the very high temperatures following the liquidus geotherm constructed through the melting profiles of peridotite (hot liquidus geotherm hereafter) (26) and for the most oxidized evolution scenarios [log( fO 2 )~IW-0.6, where IW is the ironwüstite buffer], metal-silicate equilibration would lead to a 4.5 times lower I/Pu ratio in the middle stages of accretion compared to the final mantle ( fig. S5B). If some part of the mantle in this intermediate stage was not involved in the subsequent accretion process and preserved throughout Earth's history, then these relics of proto-Earth's mantle could explain the lower 129 Xe*/ 136 Xe* Pu in the plume mantle. However, elemental and isotopic evidence suggests that Earth mainly accreted under less oxidizing conditions (i.e., between IW-4.5 and IW-1.5) (1,2,22,23,27,28), not to mention that accretion under oxidizing conditions would yield a much smaller terrestrial core than observed (29). Considering only the plausible range of fO 2 relevant to Earth's accretion (starting from a value between IW-4.5 and IW-1.5), even in the hot liquidus geotherm case, our simulations reveal that iodine is unlikely to be sufficiently fractionated from Pu during a homogeneous volatile accretion to explain the 2.8 times lower I/Pu recorded in the plume mantle ( fig. S5B).
Incorporation of intrusions from Earth's liquid outer core and/ or suspended liquid metallic droplets into part of the silicate mantle has been proposed as the possible source of the plume mantle (30). If metal-silicate equilibration happened at temperatures following the cold liquidus geotherm, then both Earth's core and the liquid metal in the magma ocean would have consistently higher I/Pu ratios than the values in the bulk mantle (Fig. 2B). In this scenario, any core/metal contribution to the plume mantle would lead to a higher 129 Xe*/ 136 Xe* Pu after decay of 129 I and 244 Pu, again the opposite of what is observed between MORBs and OIBs. Accretion under the hot liquidus geotherm, treated as the upper limit for the equilibrium temperatures, can produce I/Pu ratio two times lower in the magma ocean liquid metal (i.e., the metal after silicate-metal equilibration but before segregation into Earth's core) or in the protocore compared to the mantle during the middle and late stages of accretion ( fig. S6). This value is only 30% lower than the factor of 2.8 difference observed between OIBs and MORBs but assumes that the plume mantle Xe derives entirely from the core/metal. In the much more likely scenario where the plume mantle initially contained some Xe, incorporation of Xe from the metal/core could not produce the 2.8 times offset observed between the MORB and plume mantle reservoirs. Overall, these results suggest that the Xe anomalies in the plume mantle are unlikely to be due to incorporation of Xe from the core and/or the presence of suspended droplets of magma ocean metal liquids that never settled to Earth's core (30,31).

Heterogeneous volatile accretion
We now consider a scenario of heterogeneous volatile accretion with core formation processes, whereby the volatile content (i.e., I/Pu ratio) of Earth's building blocks is no longer kept constant over time. Moderately volatile siderophile elements in Earth's mantle suggest that volatile elements may be delivered to Earth in the later stages of the accretion process (2,(32)(33)(34), implying that the material accreted to Earth during early stages of accretion (phase 1) likely had a lower I/Pu ratio than material accreted in the later stages of accretion (phase 2) (see Fig. 2C). Specifically, our phase 2 includes not only the later stage of core formation process but also the late veneer (i.e., the last~0.5% mass of accretion that dominates the highly siderophile element budget of bulk silicate Earth), as they are thought to have delivered the same volatile-rich (high I/Pu ratio) building block materials to Earth (1). Note that the late veneer does not participate in the core formation process, and the additional volatiles it delivers are directly added to the mantle. In these heterogeneous accretion scenarios, the I/Pu ratio in Earth's mantle at the end of its accretion history (end of phase 2) could be much higher than the I/Pu ratio in the proto-Earth's mantle at the end of phase 1. At the CMB, a small portion of the proto-Earth's mantle before addition of volatile-rich material, if isolated from the subsequent accretion process and Earth's dynamical evolution (30,35,36) thereafter, could be the source of the depleted I/Pu mantle reservoirs sampled by OIBs (Fig. 2C), while the rest of the mantle enriched by the subsequent phase 2 accretion processes would represent the reservoirs with higher I/Pu ratio sampled by MORBs. Figure 2C shows an illustrative calculation of a heterogeneous volatile accretion history where addition of volatile-rich materials (with high I/Pu) starts after 60% of Earth's accretion. In this example, the material accreted in phase 2 only needs to have~6 times higher I/Pu than that accreted in phase 1 to explain the 129 Xe*/ 136 Xe* Pu contrast between OIBs and MORBs. Evidently, onset of phase 2 at a later time would require a higher I/Pu in the late accreted material to still explain the 2.8 times difference in Xe isotope ratios between MORBs and OIBs. The relationship between the timing of onset of phase 2 (no later than 99.5 wt % of Earth's final mass accreted as it includes late veneer) and the mantle's I/Pu enrichment (i) after core formation processes (which could be treated as the endmember scenario in which the late veneer happened after 129 I extinction and cannot deliver any 129 Xe to Earth) and (ii) after the entire accretion history are respectively shown in fig. S8 and Fig. 3, along with independent constraints on these variables. Although some earlier works suggested that a large fraction / X silicate ) of iodine and plutonium between silicate melt and metal liquid during core formation processes. The simulation used the fiducial multistage model (i.e., cold liquidus geotherm, 60-GPa maximum pressure, and reduced conditions; see Materials and Methods for details). On the basis of these partition coefficients, change in the I/Pu ratio in terrestrial reservoirs was calculated assuming (B) a homogeneous and (C) a heterogeneous volatile accretion scenario. In (B), all materials accreted to Earth have identical I/Pu ratio, and core formation results in nearly invariant mantellic (mantle and silicate melt in magma oceans) I/Pu ratios during accretion. The high I/Pu in metal fractions (core and liquid metal in magma oceans) rules out a core/metal contribution to the plume mantle as the source of the low 129 Xe*/ 136 Xe* Pu . The I/Pu ratios (y axis) are normalized to the bulk I/Pu ratio of the building block materials. In (C), materials in the later stages of accretion (phase 2) are enriched in iodine compared to those in the earlier stages of accretion (phase 1; see enrichment factor on each colored curve). The y axis shows the mantle I/Pu ratio normalized to the bulk I/Pu ratio of the material accreted during phase 1. Accretion of volatile-enriched materials in phase 2 can result in an increase in I/Pu ratio in the final mantle relative to the mantle at the end of phase 1 by a factor of 2.8 (horizontal dashed line) or more, which is the lowermost value needed to explain the 129 Xe*/ 136 Xe* Pu offsets between the MORBs and OIBs.
(30 to 40 wt %) of Earth's building blocks were volatile-rich, recent estimates (3) based on a multielement comparison between the isotopic composition of bulk Earth and meteorites indicate that volatile-rich materials account for <~15 wt % of the bulk Earth and most likely only~4 wt %. On Fig. 3, the full range of allowable mass fraction of accretion of volatile-rich materials is thus shown to the right of the gray vertical line and reveals that I/Pu ratios in phase 2 materials must be more than an order of magnitude higher than in phase 1 materials to explain the minimum I/Pu differences between MORBs and OIBs (horizontal black line).
The high iodine enrichment, relative to Pu, in phase 2 materials indicates that accretion histories for Earth involving only chondrites are very unlikely. Deconvolving I-decay and Pu-fission Xe (see Materials and Methods) from thoroughly compiled available Xe isotopic data of different meteorites (tables S3) only suggests a small range (within a factor of~10) of 129 Xe*/ 136 Xe* Pu among chondrites (tables S4 and S5). Independently, the latest data in carbonaceous, enstatite, and ordinary chondrites show that iodine abundances do not vary substantially across the board (37), and the vast majority has iodine content within a factor of 10 of each other (38). Last, while some studies also suggested that comets could accrete to Earth during the late veneer (39,40), the fact that the iodine content of comets (67P/Churyumov-Gerasimenko) appears to be of same order of magnitude as that of chondrites (37,38,40) further strengthens the conclusion that chondrites are very unlikely to be the primary building block materials during phase 1.

Accretion model of differentiated planetesimal
Compared to chondrites, achondrites have 40 to 7500 times lower 129 Xe*/ 136 Xe* Pu (tables S4 and S5) and, thus, much lower I/Pu ratios. This high extent of iodine deficiency compared to plutonium in achondrites is in line with the idea that most parent bodies of achondrites are highly volatile-depleted and refractory-enriched (5). In a scenario where the parent bodies of achondrites represent Earth's main building blocks in phase 1, the Xe isotopic offsets between MORBs and OIBs can be readily generated by late accretion of any type of chondrites. Most achondrites originated from volatile-depleted and refractory-enriched asteroids (5), which are themselves the remnants of differentiated planetesimals. Meanwhile, volatile-poor differentiated planetesimals have been suggested to be the archetypal first-generation planetary embryos formed in the innermost Solar System within <2 Ma of its formation (as defined by the age of calcium, aluminum-rich inclusions) (41). Our results thus support a heterogeneous accretion history of Earth, whereby early-formed, volatile-poor differentiated planetesimals represent Earth's main building blocks (6) (>85%), and chondrites, which are comparatively volatile-rich, represent late-accreted materials and only account for a small fraction (<15%) of Earth's mass (Fig. 4).
Isotopic constraints indicate that Earth accreted from both noncarbonaceous (volatile-poor inner Solar System material) and carbonaceous (volatile-rich outer Solar System material) materials (3,42). Our results further suggest that the noncarbonaceous material primarily consists of early differentiated planetesimals so as to create large enough volatile difference among Earth's building blocks. Because silicon preferentially partitions into Earth's core over magnesium at high pressure, the Mg/Si ratio in the mantle would increase during the accretion process ( fig. S9). As a result, at the CMB, the portion of proto-Earth's mantle isolated from the later stages of accretion would have a higher viscosity (due to lower Mg/Si ratio) compared to the final mantle (43): a feature that could help preserve the chemical records of early accretion of volatilepoor differentiated planetesimal at the base of the mantle until today. Our heterogeneous accretion model also readily explains the fact that present-day plume mantle reservoirs have a uniformly low ratio of potassium (a volatile, incompatible element) over uranium (a refractory, incompatible element) compared to the MORB mantle, without resorting to complex geochemical or geodynamical processes such as crust recycling or mantle convection (44,45). Last, our model sheds light on the origin of Earth's water, as it requires that chondrites represent the main material delivered to Earth in the last 1 to 15% of its accretion (3,6). Independent constraints from Mo nucleosynthetic anomalies require these late accreted materials to come from the carbonaceous supergroup (3). Together, these results indicate that carbonaceous chondrites must have represented a non-negligible fraction of the volatile-enriched materials in phase 2 and, thus, play a substantial role in the water delivery to Earth. This is in agreement with recent evidence from H and N isotopes that independently point to delivery of 4 to 15% of CI-like materials to Earth during the late stages of accretion Each curve represents a different extent of iodine enrichment (curve labels) in the phase 2 materials relative to phase 1 materials. As phase 2 includes late veneer (last 0.5 wt % of accretion process), its starting time should be earlier than 99.5 wt % of Earth's final mass accreted. Simulation conditions as in Fig. 2. The horizontal black line denotes a I/Pu ratio in the final mantle that is 2.8× higher than at the end of phase 1. The vertical line shows the upper limit (15 wt %) on the mass fraction of volatile-rich material accreted by Earth late in its history, as constrained by nucleosynthetic anomalies data (εMo and correlations with εCr, εNd, and εZr) on meteorites (3). Large I/Pu enrichments in the MORB mantle (compared to the mantle at the end of phase 1) are only achieved if volatile-poor, differentiated planetesimals (with I/Pu ratios of 40 to 7500× lower than chondrites; tables S4 and S5) constitute the main building blocks of Earth during phase 1 (blue area). In contrast, accretion histories involving undifferentiated chondrites (whose I/Pu ratios only differ by a factor of~10; tables S4 and S5) as the main material in phase 1 (red area) fail to yield a 2.8× difference in mantellic I/Pu ratio. To date, there are no known planetary materials with I/Pu enrichment between 12 and 40× (tables S4 and S5) although these materials could produce large I/Pu enrichments in the MORB mantle.
(4). A coherent picture is therefore emerging where Earth accreted mostly dry and gained its water and volatiles only late in its accretion history, highlighting the impact that small differences in formation history (e.g., changes in the nature of minor building blocks) can have on a planet's evolution and habitability.

FPMD simulations
Calculations were performed using the two-phase FPMD method following (18)(19)(20) and using the ab initio total energy and molecular dynamics program VASP (46). The projector-augmented wave potentials (47,48) were used together with the generalized gradient approximation of the exchange-correlation potential (49) in which 6 valence electrons are considered for O (2s 2 2p 4 ), 2 for Mg (3s 2 ), 14 for Fe (3p 6 3d 6 4s 2 ), 4 for Si (3s 2 3p 2 ), 7 for I (5s 2 5p 5 ), and 16 for Pu (6s 2 6p 6 6d 2 5f 4 7s 2 ). The plane-wave basis set cutoff was 400 eV. The accuracy for electronic self-consistent iteration was 10 −4 eV. The Brillouin zone sampling was performed only at the gamma point, and the Fermi-Dirac smearing was used to consider the temperature effect. Molecular dynamics simulations were performed in the NVT canonical ensemble (i.e., number of atoms, volume, and temperature remain constant during simulation) with a time step of 1 fs. Because of the complex nature of the 5f electron states for Pu, the splitting of the orbital energy level under symmetrical crystal structures might cause spin crossover (50,51). However, it would be hard to form symmetrical configuration around Pu in silicate melt or iron metal liquid phases considering the size of Pu. Furthermore, the high pressure and temperature would favor Pu in the nonmagnetic low spin state even if the 5f could split in liquid (51). Therefore, in consideration of the time consumptions of our long-time (260 atoms and more than 60,000 steps) FPMD running, we did not perform spin-polarized calculations in this work.
The simulation box contained 260 atoms. Oxygen (O), Mg, Fe, and Si represent 256 of these atoms, and the other 4 atoms are I or Pu. The relative proportions of O, Mg, Fe, and Si were set to match bulk Earth compositional models from (52). Most simulations used the O-bearing Earth composition model, which contains 127 O, 51 Mg, 44 Fe, and 34 Si atoms, while the simulation run named "I3" (see table S1) used the Si-bearing Earth composition model, which has same number of Mg and Fe atoms but contains 123 O and 38 Si atoms. The size effect (i.e., number of atoms) has been tested in (19).
The atoms were randomly distributed in a simulation cell, and 30,000 time steps were used to allow for segregation of liquid iron and silicate melt and for the system to reach equilibrium. Another 30,000 to 50,000 steps were used to calculate the average compositions of the liquid iron and silicate melt phases. For calculating composition of the Fe phase, a polyhedron alpha shape (a term used in computational geometry with smaller values of alpha describing more details of an object) was first constructed for the Fe cluster using the randomized incremental algorithm (53). The alpha shape is derived from the convex hull of iron cluster using Delaunay triangulation of all Fe atoms and describes the detail of the iron cluster surface by sifting the Delaunay tetrahedral sets according to the radius of their circumscribed sphere (fig. S1). The periodic boundary condition of the simulation box must be considered here because the iron cluster tends to spontaneously conform to a cylindrical shape, especially at high pressures and temperatures, probably due to a smaller surface area of the cylinder (considering periodic boundary) compared to sphere. Determination of whether an atom (Mg, Si, O, Pu, and I) resided inside the polyhedron was then done using the random ray crossing algorithm (54). More details could be found in (18,19).
As discussed in (18), because of the small size of the system size, one needs to consider the surface effect of iron cluster. The surface effect refers to the potential ambiguity in determining whether the iron atoms on the surface of an iron cluster belong to the silicate cluster or the iron cluster. The irregular shape of the Fe cluster makes it difficult to estimate the number of Fe irons in each phase. In keeping with (20), we used valence balance of the silicate phase to consider this effect. By distinguishing whether an atoms is included in the polyhedron, we already know the number of Si, O, Mg, Pu, and I that belong to the iron cluster, N Metal

. Schematic representation of the heterogeneous accretion history of Earth that is consistent with the more siderophile behavior of I and Pu at high P-T conditions (this work).
As core formation alone does not result in I/ Pu fractionations sufficient to explain the~3 times lower 129 Xe*/ 136 Xe* Pu ratio observed in OIBs compared to MORBs (Fig. 2B), a scenario of heterogeneous accretion has to be invoked in which volatile-depleted differentiated planetesimals constitute the main building blocks of Earth for most of its accretion history (phase 1), before addition of, comparatively, volatile-rich undifferentiated materials (chondrite and possibly comet) during the last stages of accretion (phase 2). Isolation and preservation, at the CMB, of a small portion of the proto-Earth's mantle before addition of volatile-rich material would explain the lower I/Pu ratio of plume mantle, while the mantle involved in the last stages of the accretion would have higher, MORB-like, I/Pu ratios. Because the low I/Pu mantle would also have an inherently lower Mg/Si ( fig. S9), its higher viscosity could help to be preserved at the CMB until today. silicate cluster, we obtain N Silicate where the valence of Pu in the silicate melt is assumed to be +3 (55,56) and that of iodine is 0 (11). Conservation of the number of atoms means that the number of Fe atoms in the iron cluster is then simply As a result, we know all the numbers of atoms in the two phases at a single snapshot. Here, we use 30,000 to 50,000 snapshots to calculate the average compositions of the liquid iron and silicate melt phases.

Iodine and plutonium partitioning
Iodine partitioning between silicate melt and iron metal liquid can be described by the reaction in Eq. 1 following (11). Although IO 3− , I 2 , and I − could coexist in the silicate glass (57), the exchange coefficient of iodine derived by using the reaction in Eq. 1 could well describe the high P-T experiment data (11). The iodine equilibrium exchange partition coefficient (K I D ) and equilibrium constant (K I ) are defined as where X denotes the mole fraction and γ is the activity coefficient. The natural logarithm of the equilibrium constant K I is the change in Gibbs free energy, which is where P is pressure (in gigapascal) and T is temperature (in Kelvin). The parameters a, b, and c are related to the changes in the standard-state entropy (ΔS 0 I ), enthalpy (ΔH 0 I ), and volume (ΔV i ) of the reaction in Eq. 1. It is noted that equilibrium constants only vary with pressure and temperature and do not depend on fO 2 conditions. We therefore directly regressed the equilibrium constant (instead of the partition coefficient) as a function of pressure and temperature. Once this functional form has been established, the exchange partition coefficient K I D can then be expressed, for any fO 2 condition, as In the above equation, the γ silicate term is incorporated into the expression of K I because activity coefficients in the silicate phase are not strong functions of composition (58) following (11). The epsilon formalism of (59) is used to describe the activity coefficient in the iron-rich metal liquid (the terms of γ Fe and γ 0 i are also incorporated into parameters a, b, and c) 29. Because the model with or without ɛ C I fits the data equally well based on F test (probability value: P > 0.05), ɛ C I was set to 0 in this work. Although the uncertainty of a I is higher than its value, the a I term cannot be neglected because it account for the entropy change of reaction (58).
Plutonium partitioning between silicate melt and iron metal liquid can be described as a dissociation reaction in Eq. 2. As trivalent Pu is verified to be enriched in silicate minerals and glass at high temperature (55,56), we assume Pu to behave similarly in metal-silicate partitioning situations (25). The plutonium equilibrium exchange partition coefficient (K Pu D ) and equilibrium constant (K Pu ) are defined as Because there are only three FPMD data points for Pu and no experimental data with detailed composition data to calculate K Pu D (17), too many parameters would result in overfitting, and we therefore assumed Parameters a Pu and b Pu were fitted to exchange coefficients by the least-squares method and yielded a Pu = 11.1 ± 1.1 and b Pu = − (56.8 ± 3.5) × 10 3 .
In theory, because equilibrium constants do not depend on fO 2 , experiments and FPMD simulations could be performed at any fO 2 conditions and still be used in core formation modeling. In practice and to avoid extrapolation, it is preferable to run experiments and FPMD simulations under conditions relevant to the core formation modeling. This is why our simulations were conducted at fO 2 values between IW-1.9 and IW-2.9. These values cover the range of values used in typical core formation models (1,22,23) and encompass our preferred accretion model whereby Earth mainly accreted from differentiated planetesimals thought to be very reduced. That range is similar to the conditions of the more oxidized piston cylinder experiment (IW-1.9) in (11) and is only~1 to 2 log units below (i.e., more reduced than) those of the DAC experiment (IW-0.5 to IW-2.0). Overall, all available data (FPMD and experiments) focused on the quite reducing conditions thought to have been relevant to core formation on Earth [IW-1.5 to IW-2.5; e.g., (23,28)]. Thus, when possible (i.e., for iodine), all data are considered to build the equilibrium constant fit as function of pressure and temperature.
A final quality control on our simulations was performed by calculating the major element equilibrium constants obtained in our FPMD simulations. As shown in fig. S10, the equilibrium constants for Si, O, and Mg from our new FPMD data are fully consistent with previous FPMD data and experimental data obtained at fO 2 ranging from IW-0.4 to IW-4.0. This excellent agreement further demonstrates (i) the adequacy of our FPMD method and (ii) the feasibility of regressing equilibrium constant to obtain partition coefficients, even with a limited number of simulations.

Core formation modeling
Core formation was modeled as a multistage process during Earth's accretion (22,60). In the fiducial model of this work, Earth was accreted to its present mass in 0.1% increments (i.e., 1000 steps). The metal fraction in the impactor was set to 0.325 in the fiducial model. Each impact was assumed to generate a magma ocean at the top of the proto-Earth, within which the metal in the impactor fully equilibrated with the entire magma ocean (including the silicate in the impactor) under the pressure and temperature of the base of the magma ocean (1,22). After equilibrium, the equilibrated liquid metal was assumed to descend rapidly to the core through the solid mantle without further equilibration (60) and the equilibrated magma ocean mixed with the solid mantle to form a new mantle (1). The above assumptions, while necessarily simplistic, were made as extreme examples to enable the expression of the largest core formation-induced I/Pu fractionation. The rationale being that if even this scenario cannot explain the I/Pu contrast between MORBs and OIBs, then more realistic scenarios (e.g., invoking partial equilibration) would not either.
During this incremental accretion process, the pressure at the base of the magma ocean increased with the fraction of mass accreted (f ), following (22) where P fin is the equilibrium pressure of the magma ocean at the last accretion. Using this equilibrium pressure, the temperature was determined on the basis of a cold liquidus geotherm constructed through the melting of chondritic mantle (21) in the fiducial model. The fraction of mantle mass involved in the magma ocean ( f mo ) was parameterized as a function of P fin following (1) where the radius of Earth (r earth ) is 6371 km, the radius of Earth's core (r core ) is 3485 km, and the pressure at CMB (P CMB ) is 135 GPa. In the fiducial model, P fin is assumed to be 60 GPa. During the multistage core formation process, the distribution of elements was calculated by mass balance using their partition coefficients. The concentrations of O, Si, Co, Cr, V, and Ni were first derived at each step. The partition coefficient of these elements were taken from (61,62) following (22). In each accreting building block, the bulk concentrations of these elements (except for O and Fe) relative to Al were those of the bulk Earth's composition assuming 5 wt % of Si in core (63), following (1). In principle, the fO 2 of Earth would vary as accretion proceeds because of loss of O to space (64), which would affect the size of Earth's core (29). Because it is hard to precisely evaluate the O escape to space during accretion, to test how fO 2 affects the accretion process, we assume that a constant core mass fraction during Earth's accretion and the concentration of O and Fe were instead controlled by predetermined evolution scenarios of the FeO content in the magma ocean. These evolution histories were constructed by linear interpolation between different meteorites as a starting point and present Earth's mantle as an ending point [ figure 1 in (22)]. The fiducial model used the FeO in high-iron enstatite (EH) chondrite as the starting point [path 2 of figure 1 in (22)]. This fiducial model illustrates a scenario in which Earth started accretion from a plausible reduced redox condition but does not mean that Earth formed through accretion of EH chondrites. Then, the distribution of Pu and I were determined on the basis of the pressure, temperature, and metal compositions for each step using Eqs. 3 and 4 and their activity coefficients.
Besides the fiducial model, a series of simulations were performed in which key input parameters were varied to test the sensitivity of these parameters on the simulation outcomes. These runs tested different (i) P fin values (60, 70, and 80 GPa) that determine the equilibrium pressure, (ii) mantle liquidus geotherm [hot (26) and cold (21)] that determines the equilibrium temperature, and (iii) evolution scenario of FeO content of the magma ocean [paths 2, 9, 10, 11, 12, 13, and 14 of figure 1 in (22)], ensuring that the fO 2 conditions are varied. A special accretion scenario of equilibrium pressure was also constructed to verify whether highpressure episodes in the middle stages of accretion could efficiently fractionate I and Pu. In this scenario, equilibrium pressures during the first 50% of the accretion was the pressure at the CMB of the proto-Earth (P fin = 135 GPa), and the pressure during the last 50% accretion linearly decreased from the highest pressure (85 GPa) to a final equilibrium pressure (40 GPa). We also tested whether the change of core size would affect I/Pu fractionation during accretion process. Heterogeneous volatile accretion simulations were performed by multiplying iodine concentration of building block materials in phase 2, while Pu concentration was assumed to remain constant throughout Earth's accretion.

Meteoric Xe analysis
Xenon has nine isotopes (i.e., 124 Xe, 126 Xe, 128 Xe, 129 Xe, 130 Xe, 131 Xe, 132 Xe, 134 Xe, and 136 Xe) whose relative abundances represent a mixture of six components (see table S3): (i) spallogenic xenon (noted here sp), (ii) initial or primordial xenon (init) that can be represented by solar wind or Phase-Q, (iii) atmospheric xenon (air), (iv) uranium-fission xenon (U), (v) plutonium-fission xenon (Pu), and (vi) iodine-decay xenon. Except for the iodinedecay component, 132 Xe is present in all other components. So, we can normalize to 132 Xe (avoiding issues from divisions by zero) and decompose the isotopic composition in a sample as the computations using Tianhe-2 supercomputers. Funding: Y.Z. is grateful for the support from Strategic Priority Research Program (B) of the Chinese Academy of Sciences (XDB18000000).