# Nitrogen isotope evidence for Earth’s heterogeneous accretion of volatiles | Panda Anku

### High-pressure experiments

All experiments at 1–1.5 GPa and 1700 °C were conducted in an end-loaded, solid media piston cylinder apparatus, using 0.5-inch diameter talc–Pyrex assemblies with stepped graphite heaters (Supplementary Data 2). For these experiments, starting materials of 30–50 wt.% metallic mixture and 50–70 wt.% silicate powder were loaded into zirconia-lined or graphite-lined Pt95Rh05 capsules. The use of zirconia-lined Pt95Rh05 capsules was to avoid the presence of C in the sample. The hot piston-in method was used to pressurize the assembly. Pressure was calibrated against the quartz–coesite and kyanite–sillimanite transitions, and a friction correction of 18% was applied to the nominal pressure. The pressure uncertainty is estimated to be better than 5% relative. We used W97Re03–W75Re25 thermocouples and temperature was controlled to ±2 °C and is accurate to ±10 °C. The experimental durations ranged from 1 to 3 h. All experiments were quenched to below 100 °C within 10–20 s by turning off electric power to the graphite heaters.

The experiments at 2000–2200 °C and 5–8 GPa were conducted in a 1000-ton DIA-type multi-anvil apparatus. The starting material was loaded into graphite capsules. Tungsten carbide (WC) anvils with 8-mm truncated edge lengths, together with a Cr2O3-doped MgO octahedron (14 mm edge length) as the pressure medium, were used to generate high pressures. The furnace is composed of a stepped cylindrical graphite or LaCrO3 heater and a ZrO2 thermal insulator. Sample pressures were estimated from the hydraulic pressure using calibrations based on the phase transitions of Bi, ZnS, and Mg2SiO4 polymorphs. Temperatures were measured using a W95Re05–W74Re26 thermocouple whose junction was located near the end of the graphite capsule. The cell assemblies were first pressurized to target pressures at room temperatures in 4 h and then heated to target temperatures at a rate of 100 °C/min. The temperature fluctuation during the experiments was approximately ±10 °C, and the experimental durations ranged from 20 to 40 mins. All experiments were quenched by switching off the electric power, and then decompressed to ambient pressure for more than 16 h.

All recovered sample capsules were sectioned longitudinally into two halves. One half was mounted in epoxy and polished for electron microprobe and Raman spectroscopy analysis, and the other half prepared for noble gas mass spectrometry analysis.

### Electron microprobe analysis

Major element compositions of the quenched metallic and silicate melts were measured with a JEOL JXA–8200 microprobe. The analyses were performed in wavelength-dispersive mode, and a PAP matrix correction was applied to the raw data. The metallic melts were analyzed with 20 kV acceleration voltage and 20 nA beam current, whereas the silicate melts were analyzed with 15 kV/10 nA. Natural and synthetic standards were used to calibrate the instrument. For the metallic melt analysis, Fe, Si, Ni, and Cu were calibrated on pure metals, S was calibrated on a synthetic pyrrhotite with well-known Fe: S ratio, and O was calibrated on magnetite. For the analysis of C in metallic melt, the samples and standards were uncoated but surrounded with silver-bearing conductive varnish to avoid charging, and Fe3C was the standard. For the silicate melt analysis, Na was calibrated on albite, Ca on wollastonite, K on orthoclase, Ti and Mn on ilmenite, Si on enstatite, Mg on forsterite, Al on spinel, P on GaP, and Fe on metallic Fe. Sulfur in the quenched silicate melts was analyzed with 50 nA beam current and 60 s peak counting time using a pyrrhotite standard. A defocused beam of 10 or 20 µm diameter was used for all standardizations and sample measurements.

### Noble gas mass spectrometry analysis

The coexisting metallic melt and silicate melt in one half of all samples were physically separated and measured for N-content and -isotopic composition, using a modified noble gas mass spectrometry at the Atmosphere and Ocean Research Institute, the University of Tokyo, which has a capacity to precisely measure N-content and -isotopes at the sub-Nanomole levels62. Big metal blobs and clean silicates were separated readily under binocular microscope, except for run PYH01 synthesized at 5 GPa and 2200 °C (Supplementary Data 2), which only contained small metal grains (10–30 µm) being heterogeneously distributed in the quenched silicate melt. For this sample, we did not measure its N-content and -isotopic composition.

Silicate chips with a mass of ~1 mg and metal chips with a mass of ~0.1 mg were prepared for each sample. Each silicate or metal was loaded into a small quartz glass tube, which was then loaded into a large quartz glass container. The large quartz glass sample container was then placed in a furnace equipped with a tungsten wire. The sample container was baked out at 150 °C overnight under vacuum to remove any atmospheric N. A high-vacuum line was used for N gas extraction and purification. The vacuum line is basically composed of four parts: a gas-extraction part connected to the sample container, a gas-purification line (AQ-line) including cold trap-1, cold trap-2, a copper oxide finger (CuO) and a Pt finger, a vacuum line (AC-line), and a vacuum line directly connected to the mass spectrometer VG3600 (CE-line) for N-measurements. The AQ vacuum line is connected to a quadrupole mass spectrometer for checking whether N-intensity is appropriate for measurements. The vacuum lines are evacuated by a turbo molecular pump and an ionic pump. About 0.2 torr and 1 torr O2 was produced for silicate sample and metal sample, respectively, by heating a copper oxide finger (CuO) to 850–900 °C. The produced O2 was sealed in the quartz glass sample container and used for oxidizing the sample. All samples were firstly heated to 200 °C to remove any potential surficial contamination and to check the leakage of the system. During heating, the excess O2 in the AQ line was absorbed by the Cu finger through decreasing temperature of the finger to 600 °C. After the excess O2 was absorbed, the AQ line was evacuated, and the copper oxide finger (CuO) was reheated to 850 °C to produce O2. After oxidizing the samples at 200 °C for 30 min, the gases in the system produced at 200 °C were then directed to the purification vacuum line (AQ) by opening the valve of the sample container for 2 min. Condensable gases such as carbon dioxide and water were trapped using liquid N at cold trap-1. The non-condensable gases such as carbon monoxide, hydrocarbons, and hydrogen in the AQ line were then reacted for 5 min with pure O2, catalyzed by a platinum foil at 1000 °C. Carbon monoxide, hydrocarbons, and hydrogen were oxidized to CO2 and H2O by the pure O2. After reaction for 5 min, excess O2 was again resorbed by copper by decreasing the furnace temperature to 600 °C for 20  min and eventually to 450 °C for another 20 min. During this time, the produced gases such as CO2 and H2O were trapped using liquid N at cold trap-2. The pressure in the purification line, which is measured by a capacitance manometer, is no more than 0.0001 torr higher than the baseline pressure. A quadrupole mass spectrometer (QMS; HAL201, Hiden Analytical) was used to check whether the sample volume is proper for analysis through determining the sample size. If the N2 intensity in QMS is in the appropriate range, then the purified gases would be measured using the high-sensitivity static vacuum mass spectrometer (a modified VG3600, VG Micromass Ltd.). The determined concentration of N released at 200 °C was negligible (less than 1 ppm). After the measurements of gases released at 200 °C, following the same procedure, the samples were then heated to 1200 °C and kept for 30 min for glasses and 60 min for metals to ensure that the samples were completely oxidized for releasing N. The purification and measurement procedure for gases released at 1200 °C was the same as the procedure at 200 °C. After each measurement finished, the VG3600 mass spectrometer was evacuated for at least 20 min.

The running standard gas for N-isotopes was reserved in a large metal container attached to CE-line. The N-isotopic ratio of the standard gas was measured before and after the sample analysis, and the N-isotopic ratio of the standard gas was checked periodically by comparison with that of the local air. After the standard/sample gas was introduced into the VG3600 mass spectrometer, the intensities of 28 and 29 were measured for 15 times, and the 28/29 ratio was extrapolated by the mass spectrometer, with a system error of ~0.3‰. Repeated analysis of the standard in a day showed that the overall reproducibility was better than 0.5‰ for N-isotopes, and repeated analysis of a same sample yielded consistent results with uncertainties of <5% relative. The N-isotopic ratio of the air was analyzed for a number of times during the course of analysis, which gave an average δ15N value of −0.03‰. The concentration and isotopic composition of N in each sample, released at 1200 °C, was corrected by subtracting the hot blank (i.e., background), and the hot blank was analyzed with glass tube that contained no sample but using the same heating duration, same amounts of oxygen, and same purification procedure as in the case for measuring the samples. An average N-isotopic composition of +2.99‰ ± 0.70‰ is obtained for the hot blank.

The error (σ) of the measured N-concentration was calculated based on the error of the N-concentration measured by the VG3600 mass spectrometer, the error of the hot blank, and the error of weighing the sample mass, which can be defined as:

$${sigma }_{{{{{rm{N}}}}}}({{{{{rm{ppm}}}}}})={{{{{rm{N}}}}}}({{{{{rm{ppm}}}}}})times sqrt{{left(frac{{sigma }_{{{{{rm{N}}}}}2}({{{{{rm{mole}}}}}})}{{{{{{rm{N}}}}}}_{2,{{{{{rm{original}}}}}}}left({{{{{rm{mol}}}}}}right)-{{{{{rm{N}}}}}}_{2,{{{{{rm{blank}}}}}}}({{{{{rm{mol}}}}}})}right)}^{2}+{sigma }_{{{{{{rm{weighing}}}}}}}^{2}}$$

(3)

In Eq. (3), N (ppm) is the blank-corrected N-concentration in the sample; and the error of sample weighing (σweighing) is usually 1% relative; N2original (mol) is the uncorrected N-concentration measured by the VG3600 mass spectrometer; σN2 (mole) is the error of background-corrected N-concentration in the sample, which can be defined as:

$${sigma }_{{{{{{rm{N}}}}}}_{2}}({{{{{rm{mol}}}}}})=sqrt{{sigma }_{{{{{{rm{N}}}}}}_{2}{{{{{rm{original}}}}}}}^{2}+{sigma }_{{{{{{rm{N}}}}}}_{2}{{{{{rm{blank}}}}}}}^{2}}$$

(4)

$${sigma }_{{{{{{rm{N}}}}}}_{2}{{{{{rm{blank}}}}}}}({{{{{rm{mole}}}}}})={{{{{rm{N}}}}}}_{2,{{{{{rm{blank}}}}}}}left({{{{{rm{mol}}}}}}right)times {sigma }_{{{{{{rm{N}}}}}; {{{{{rm{sensitivity}}}}}}; {{{{{rm{in}}}}}}; {{{{{rm{air}}}}}}}}$$

(5)

where the ({sigma }_{{{{{{rm{N}}}}}; {{{{{rm{sensitivity}}}}}}; {{{{{rm{in}}}}}}; {{{{{rm{air}}}}}}}}) is about 9.9%. The error of N-isotopic composition of the sample was calculated based on the error of sample gas δ15N (‰) measured by the VG3600 mass spectrometer, the error of reproducibility of standard N2 (0.5‰), and error of δ15N of the air (0.43‰), which can be defined as:

$${sigma }_{{delta }^{15}N({{textperthousand }})}=sqrt{{sigma }_{{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}{left({{textperthousand }}right)}_{{{{{{rm{MS}}}}}}}}^{2}+{sigma }_{{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}{left({{textperthousand }}right)}_{{{{{{rm{Standard}}}}}}}}^{2}+{sigma }_{{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}{left({{textperthousand }}right)}_{{{{{{rm{Air}}}}}}}}^{2}}$$

(6)

$${sigma }_{{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}{left({{textperthousand }}right)}_{{{{{{rm{MS}}}}}}}}=sqrt{{sigma }_{{28/29{{{{{rm{N}}}}}}}_{{{{{{rm{Standard}}}}}}-1}}^{2}/2+{sigma }_{{28/29{{{{{rm{N}}}}}}}_{{{{{{rm{MS}}}}}}}}^{2}+{sigma }_{{28/29{{{{{rm{N}}}}}}}_{{{{{{rm{Standard}}}}}}-2}}^{2}}/2$$

(7)

where ({sigma }_{{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}{left({{textperthousand }}right)}_{{{{{{rm{MS}}}}}}}}) is the error of δ15N (‰) of sample gas measured by the VG3600 mass spectrometer; ({sigma }_{{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}{left({{textperthousand }}right)}_{{{{{{rm{Standard}}}}}}}}) is the error of the reproducibility of the standard gas; repeated analysis of the standard in a day showed that the overall reproducibility was better than 0.5‰; ({sigma }_{{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}{left({{textperthousand }}right)}_{{{{{{rm{Air}}}}}}}}) is the error of the δ15N (‰) value of the air, which is 0.43‰; ({sigma }_{{28/29{{{{{rm{N}}}}}}}_{{{{{{rm{MS}}}}}}}}) is the standard deviation of 15 times measurements of 28/29 ratio of sample gas extrapolated by the mass spectrometer, which is approximately 0.3‰; ({sigma }_{{28/29{{{{{rm{N}}}}}}}_{{{{{{rm{Standard}}}}}}-1}}) and ({sigma }_{{28/29{{{{{rm{N}}}}}}}_{{{{{{rm{Standard}}}}}}-2}}) are the standard deviation of 15 times measurements of 28/29 ratio of the standard gas before and after each sample analysis.

### Raman spectroscopy analysis

To obtain information on the speciation of N–C–H–O in the silicate melt, Raman spectra of some silicate glasses were collected. Micro-Raman spectra were recorded in back-scattering geometry using a Horiba LabRAM HR UV spectrometer with CCD detector, 1800 mm−1 grating, 50× objective, and confocal mode. A confocal pinhole of 100 μm was used, which limits spectral resolution to 3.5 cm−1. The 514.5 nm line of an Ar+ ion laser with an output power 0.1 W was used for excitation. Spectra were collected from 200 to 4500 cm−1 with an acquisition time of 2 × 300 s for each range, in order to detect N–C–H–O species with very low concentrations. Typical Raman spectra of some silicate glasses at fO2 of ~IW−0.5 are shown in Supplementary Fig. 2. We encountered significant fluorescence when measuring the reduced silicate melts. Such fluorescence could significantly mask the weak peaks of N–C–H–O species; we therefore only reported the Raman spectra that were not affected by fluorescence (Supplementary Fig. 2). Actually, the variation of N–C–H–O species in silicate melt as a function of fO2 has been extensively measured and discussed in the literature36,38,39,53,63,64, and their results are consistent with one another. Here we only cited these previous studies to discuss the effect of the variation of N-species in silicate melt on the N-isotopic fractionation (see the main text).

### Major element compositions of experimental products and calculation of sample fO2

In all experiments, the basaltic melts were quenched into glasses, while the mantle pyrolitic melts had a dendritic texture (Supplementary Fig. 1). Major element compositions of the quenched metallic and silicates are given in Supplementary Data 3, 4. The ratio of non-bridging oxygens to tetrahedral cations (NBO/T) of the quenched silicate melt was 0.2–3.1, and the Fe-rich metallic melts contained 87.0–98.8 wt.% Fe, 0–6.8 wt.% Si, 0–1.8 wt.% S, and 0–11.5 wt.% C.

The fO2 prevailing in our experiments was calculated from the coexistence of Fe-rich metallic melt and silicate melt with finite FeO content using the following equilibrium:

$${{{{{rm{FeO}}}}}}({{{{{rm{silicate; melt}}}}}})={{{{{rm{Fe}}}}}}({{{{{rm{metallic; melt}}}}}})+{1/2{{{{{rm{O}}}}}}}_{2}$$

(8)

from which the fO2 relative to the Fe-FeO buffer at any given PT can be defined as:

$$Delta {{{{{rm{IW}}}}}}=2{{log }}({a}_{{{{{{rm{FeO}}}}}}}/{a}_{{{{{{rm{Fe}}}}}}})=2{{log }}({X}_{{{{{{rm{FeO}}}}}}}{{{{{{rm{gamma }}}}}}}_{{{{{{rm{FeO}}}}}}}/{X}_{{{{{{rm{Fe}}}}}}}{{{{{{rm{gamma }}}}}}}_{{{{{{rm{Fe}}}}}}})$$

(9)

aFeO represents the activity of FeO in the silicate melt; aFe represents the activity of Fe in the metallic melt; XFeO and XFe are molar fractions of FeO in the silicate melt and Fe in the metallic melt, respectively; γFeO and γFe are activity coefficients of FeO in the silicate melt and Fe in the metallic melt, respectively. Calculations of fO2 using both ideal (γFeO  = 1 and γFe = 1; ideal fO2) and non-ideal solution models (non-ideal fO2) were performed. The fO2 calculation using the non-ideal solution model was performed assuming γFeO = 1.5 (refs. 17,65). γFe was calculated using the ε-approach and the online “Metal Activity Calculator” (http://www.earth.ox.ac.uk/~expet/metalact/) provided by the University of Oxford, which take into account the non-ideal interaction between all components in the Fe-rich metallic melt66. The calculated non-ideal fO2 values are between IW and IW–5. The calculated ideal fO2 values are between IW–0.7 and IW–6, which are 0.7–1 log units lower than the non-ideal fO2 values. The calculated fO2 values decrease with decreasing the FeO content of the silicate melt or increasing the Si content of the metallic melt, consistent with previous results17.

### Sample N-content and δ15N

The δ15N of the starting Fe7N3 was −7.9 ± 1‰. The N-contents of the metallic and silicate melts are 43–14293 ppm and 44–4620 ppm (by wt.), respectively (Supplementary Data 3, 4). In two experiments (LY24 and LY25; Supplementary Data 2), where ~50 wt.% Fe7N3 was added in the starting material, the metallic melts contained ~1.1 and 1.4 wt.% N, respectively. These values are comparable to the N-solubility in Fe-rich metallic melts determined under similar conditions37, implying that these two experiments must have been N2-saturated. The δ15N of metallic melts ranged from −7.0‰ to +7.6‰, and the δ15N of silicate melts ranged from −7‰ to +0.42‰ (Supplementary Data 3, 4). Except for the two N2-saturated experiments, mass-balance calculations of bulk sample N content and δ15N were performed for all other experiments (Supplementary Data 2). The results show that the bulk sample N-contents vary from 650 to 3000 ppm, with most values (780–2200 ppm) generally consistent with the N-mass added in the starting materials. Most of the bulk sample δ15N values vary between −3‰ and −7‰, with the δ15N values of three multi-anvil experiments above +0.2‰ (N-18, N-19, and N-20; Supplementary Data 2). The deviation of the bulk sample δ15N from that of the starting Fe7N3 could be explained by air N2 contamination when loading silicate powder into the graphite capsule, the presence of N in the starting metallic Fe, and/or δ15N inhomogeneity of the starting Fe7N3. Nevertheless, the positive δ15N values of the three multi-anvil experiments are difficult to explain. However, as shown in the main text, the calculated ({D}_{N}^{{{{{{rm{metal}}}}}}/{{{{{rm{silicate}}}}}}}) and ({triangle }^{15}{{{{{rm{N}}}}}}^{{{{{{rm{metal}}}}}}-{{{{{rm{silicate}}}}}}}) of these three experiments are fully consistent with those of the other experiments.

### N-loss during experiments?

Previous studies showed that N-loss occurred during the run in some of the high PT experiments36,37,40,41. The N-loss could be either related to the storage of some N in the porosity of the graphite capsules or to N diffusive loss through the graphite (-Pt) capsule walls41. We note that in previous experiments, ~0.5–2 wt.% N was put in the starting materials, so that the samples were usually saturated with a N2 vapor phase. However, in our experiments, we added only ~1000–2000 ppm N in the starting materials, which are far below the solubility limit37. Our experiments were therefore not saturated with a N2 vapor phase, except for two experiments (LY24 and LY25; Supplementary Data 2). As mentioned above, the reconstructed bulk N content in our samples ranged from 650 to 3000 ppm, with most values between 780 and 2200 ppm and consistent with the N-mass added in the starting materials (Supplementary Data 2). We therefore believe that N-loss in our experiments was limited. Some experiments, as shown in Supplementary Data 2, even gained N, which could be caused by air N2 contamination during loading sample materials into the capsules, and the presence of N in the starting metallic iron. And our experiments LY24 and LY25 were still saturated with a N2 vapor after running 90 mins. The factors that control the N-loss or the N-loss rates are elusive, but at least, the absence of a N2 vapor phase should effectively suppress the N-loss rates or the storage of N2 gas in the porosity of the graphite capsule walls.

### Attainment of equilibrium partitioning

The experimental durations ranged from 60 to 180 mins at 1700 °C, and from 10 to 20 mins at 2000–2200 °C. The observation that ∆15 Nmetal-silicate are time-independent over a small fO2 range (Supplementary Fig. 6) demonstrates 60 mins being sufficient for approaching metal/silicate N-isotopic equilibrium at 1700 °C. A time series of experiments are usually not performed at temperatures of 2000–2200 °C, because the diffusion of elements at such high temperatures are thought to be very fast. If metal/silicate N-isotopic equilibrium at 1700 °C was approached within 60–180 min, there is no doubt that N-isotopic equilibrium at 2000–2200 °C must have also been approached within 10–20 min, because the diffusion coefficient of an element in a given system increases by orders of magnitude with temperature increasing from 1700 to 2000–2200 °C. In addition, as shown in Fig. 2 and summarized in Eq. (1), our ({D}_{{{{{rm{N}}}}}}^{{{{{rm{metal}}}}}/{{{{{rm{silicate}}}}}}}) are consistent with previous data, also indicating equilibrium partitioning of N in our experiments.

### Model N-partitioning and -isotopic fractionation during Earth’s accretion and core-formation

Using the multistage core-formation model with inputs of the Grand Tack N-body simulation results15,49, and using our Eqs. (1) and (2), we modeled the N-partitioning and -isotopic fractionation between Earth’s core and mantle. We first considered Earth accreted its first 60% mass through the collisions of reduced, EC-like impactors, and then its last 40% mass through the collisions of increasingly oxidized impactors. Reduced, EC-like impactors formed at heliocentric distances of <0.9–1.2 AU with δ15N = −30‰, while increasingly oxidized impactors originated from great heliocentric distances (1.2–3 AU). Since the solar system δ15N increases with the heliocentric distance4,20, the increasingly oxidized impactor δ15N should increase from −30‰ (EC’s value) to a slightly positive value at 3 AU. We used a δ15N value of +5‰ for the last impactor, which added Earth’s last 10% mass through the Moon-forming giant impact, and a δ15N value of +5‰ represents the mixing of EC (δ15N = −30‰) and CC (δ15N = +40‰) with a same N content and mass ratio of ~1 : 1 as constrained by Mo isotopes for the Moon-forming giant impactor19. The N content in the impactors depends a number of conditions. For example, large and oxidized impactors may hold more N than small and reduced ones, and N is mainly stored in impactor cores42. For an internally differentiated planetesimal67, i.e., no surface magma ocean, more volatiles could be retained in the planetesimal, no matter the planetesimal is oxidized or reduced. However, even oxidized planetesimals may have undergone severe evaporative volatile loss, if the surface magma ocean lasted sufficiently long. Iron meteorites contain a few ppm to 200 ppm N21, whereas basalts in asteroids may contain up to 10 ppm N68. However, both iron meteorites and asteroid basalts may have experienced extensive volatile loss, considering the pressure-dependence of N solubility in both metallic and silicate melts37,69. In our model, we varied the bulk N-content from 50 ppm for small and reduced impactors to 150 ppm for large and oxidized impactors.

We also considered the other factors that potentially affect the N-content and -isotopic composition of Earth’s different reservoirs. Since N is a strong volatile element, we have to consider the degassing of the silicate magma ocean. Degassing of the silicate magma ocean may have inevitably formed a proto-atmosphere7. However, for Earth-sized planets, it is unlikely that the entire silicate magma ocean is in equilibrium with the proto-atmosphere in volatile partitioning70. We defined a factor Φ, which represents the mass percentage of silicate magma ocean that is in equilibrium with the proto-atmosphere. We fix Φ at 100% at Earth’s first ~50% accretion, but varied it from 50% to 5% at Earth’s ~50–100% accretion. In addition, the N-isotopic fractionations during degassing of the hot magma ocean are thought to be inconsiderable, as even the degassing of mid-ocean ridge basalts may have not caused N-isotopic fractionations71. A kinetic process during degassing of a hot silicate may cause significant N-isotope fractionation;72 however, it remains unfeasible to quantify such effect during the degassing of a hot magma ocean. We therefore did not consider kinetic process-induced N-isotope fractionation in our model. Hydrodynamic escape of light gases in the early atmosphere to space may have resulted in preferential loss of 14N in space; however, such effect on Earth’s δ15N is difficult to quantify and the hydrodynamic escape model is inconsistent with the abundance and isotopic composition of Earth’s atmospheric xenon22. We therefore did not consider hydrodynamic escape of Earth’s proto-atmosphere in our model. However, it is important to consider catastrophic loss of the Earth’s proto-atmosphere during impacting73, because it surely affects Earth’s total N budget. Depending on a number of parameters, such as impact velocity and angle and the impactor to target mass ratio, the loss fraction of N in the proto-atmosphere during each impacting is difficult to quantify but could vary from <5% to 100% (refs. 73,74). In the Moon-forming giant impact, ~10–50% of the growing Earth’s atmosphere could have been lost from the immediate effects of the collision73. We fix a N-loss fraction of 100% during Earth’s first ~50% mass accretion, and from 80% to 10% during the last ~50% accretion, because small planets more readily lose their atmosphere during impacting74. Finally, we considered the delivery of completely oxidized or CI chondrite-like materials from beyond 6–7 AU, which contain 1000 ppm N with δ15N = +40‰ (refs. 50,75), to Earth’s magma ocean after Earth accreted its 60% mass. No metal–silicate segregation occurred when CI chondrite-like material added in the magma ocean, but the N in the silicate magma ocean participated in the subsequent core-formation events when metal-bearing impactors were accreted15. Inefficient emulsification of impactor cores may have occurred during Earth’s accretion of large and oxidized impactors. We used the degree of core–mantle disequilibrium during accretion as that in ref. 49. We also considered the light element content of the metallic melt that segregated into Earth’s core; we used a C-content of 0.2 wt.% and S-content of 1.5 wt.% (see below), and Si- and O-contents following the model in refs. 15,49.

The PTfO2 conditions of metal–silicate equilibration are critical for the resulting mantle and core compositions because of the dependence of ({D}_{N}^{{{{{{rm{metal}}}}}/{silicate}}}) and ∆15 Nmetal-silicate on PTfO2. We used the approach used by refs. 15,49 for modeling multistage core formation. For each impact-induced core formation event, the metal–silicate equilibration pressure Pe is a constant fraction of the target’s evolving core–mantle boundary pressure:

$${P}_{e}={f}_{P}times {P}_{{{{{{rm{CMB}}}}}}}$$

(10)

where fP is a constant proportionality factor; PCMB, the core–mantle boundary (CMB) pressure at the time of impact. The metal–silicate equilibration temperature, Te, lies between the peridotite liquidus and solidus at the equilibration pressure Pe. The equilibration fO2 varied from ~IW-5 to IW-2 with the accreted impactors changing from reduced to oxidized composition.

In an ith stage of collisional accretion, the distribution of N among the reservoirs of Earth’s core, silicate magma ocean, and atmosphere must follow a mass balance:

$${M}_{i}^{{{{{rm{N}}}}}-{{{{{rm{atm}}}}}}}+{M}_{i}^{{{{{rm{N}}}}}-{{{{rm{silicate}}}}}}+{M}_{i}^{{{{{rm{N}}}}}-{{{{{rm{metal}}}}}}}={M}_{i}^{{{{{rm{N}}}}}-{{{{{rm{bulk}}}}}}}$$

(11)

where ({M}_{i}^{{{{{rm{N}}}}}-{{{{{rm{atm}}}}}}}), ({M}_{i}^{{{{{rm{N}}}}}-{{{{{rm{silicate}}}}}}}), and ({M}_{i}^{{{{{rm{N}}}}}-{{{{{rm{metal}}}}}}}) are the N mass in the atmosphere, silicate melt, and metallic melt, respectively; ({M}_{i}^{{{{{rm{N}}}}}-{{{{{rm{bulk}}}}}}}) is the total N mass that participates in the equilibrium partitioning. Equation (11) can be further written as:

$${{M}_{i-1}^{{{{{rm{N}}}}}-{atm}}-{M}_{i}^{{{{{rm{N}}}}}-{{{{{rm{los}}}}}s}}+C}_{i-1}^{{{{{rm{N}}}}}-{{{{{rm{silicate}}}}}}}{times M}_{i-1}^{{{{{{rm{silicate}}}}}}}+{C}_{i}^{{{{{rm{N}}}}}-{{{{{rm{imp}}}}}}}times {M}_{i}^{{{{{{rm{imp}}}}}}}\={C}_{i}^{{{{{rm{N-{silicate}}}}}}}{times M}_{i}^{{{{{{rm{silicate}}}}}}}+{C}_{i}^{{{{{rm{N}}}}}-{{{{{rm{metal}}}}}}}{times M}_{i}^{{{{{{rm{metal}}}}}}}+{M}_{i}^{{{{{{rm{atm}}}}}}}$$

(12)

$${C}_{i}^{{{{{rm{N}}}}}-{{{{{rm{metal}}}}}}}={D}_{{{{{rm{N}}}}}(i)}^{{{{{{rm{metal}}}}}}/{{{{{rm{silicate}}}}}}}times {C}_{i}^{{{{{rm{N}}}}}-{{{{{rm{silicate}}}}}}}$$

(13)

In Eqs. (12) and (13), ({M}_{i}^{{{{{rm{N}}}}}-{{{{{rm{loss}}}}}}}) denotes the mass of N lost during impacting; ({C}_{i-1}^{{{{{rm{N}}}}}-{{{{{rm{silicate}}}}}}}) denotes the N concentration in the silicate melt at the (i-1)th stage of accretion; ({M}_{i-1}^{{{{{{rm{silicate}}}}}}}) denotes the mass of silicate melt at the (i-1)th stage of accretion; ({C}_{i}^{{{{{rm{N}}}}}-{{{{{rm{imp}}}}}}}) denotes the N concentration in the impactor; ({M}_{i}^{{{{{{rm{imp}}}}}}}) is the mass of the impactor; ({C}_{i}^{{{{{rm{N}}}}}-{{{{{rm{metal}}}}}}}) and ({C}_{i}^{{{{{rm{N}}}}}-{{{{{rm{silica}}}}}{te}}}) denote the N concentration in the metallic and silicate melts, respectively. ({D}_{{{{{rm{N}}}}}(i)}^{{{{{{rm{metal}}}}}}/{{{{{rm{silicate}}}}}}}) are a function of pressure, temperature, fO2, and compositions of the metallic and silicate melts, as summarized in Eq. (1) in the main text. In order to model the N partitioning between the atmosphere and the silicate, the N solubility in silicate melt is needed. If the mantle is in full equilibrium with the atmosphere, its N concentration should be equal to the N solubility corresponding to the partial pressure of N2 in the atmosphere. Following the model of Libourel et al.53, the N solubility in silicate melt (({S}_{i}^{{{{{rm{N}}}}}-{{{{{rm{silicate}}}}}}},{{{{{rm{ppm}}}}}})) can be written as:

$${S}_{i}^{{{{{rm{N}}}}}-{{{{{rm{silicate}}}}}}}({{{{{rm{ppm}}}}}})=0.0611times {P}_{{{{{rm{N}}}}}}+5.97times {10}^{-10}times {{{f{{{{rm{O}}}}}}}_{2}}^{-frac{3}{4}}times {{P}_{{{{{rm{N}}}}}}}^{1/2}$$

(14)

where PN is the partial pressure of N in the atmosphere. Generally, the partial pressure of a volatile element E (PE) on the surface of magma ocean can be expressed as:

$${P}_{E}=rtimes {M}_{i}^{{{{{{rm{atm}}}}}}}times {g}_{i}/{A}_{i}^{{{{{{rm{surf}}}}}}}$$

(15)

where gi is the gravitational acceleration, ({A}_{i}^{{{{{{rm{surf}}}}}}}) is surface area of the planet, and r is the mass number ratio of the volatile species and the element of interest (for N2, r = 1). It is worth noting that in Eq. (14), the fO2 denotes to the fO2 of the surface magma ocean. Following the method of ref. 76, we calculated the surface magma ocean fO2, which varies from IW−3.6 to IW+3.0 during Earth’s early to late accretions in our model. The fO2 variation could significantly affect the composition of the proto-atmosphere77. The combination of Eqs. (12)–(15) can be used to constrain the N distribution in the atmosphere, silicate magma ocean, and core in the framework of a multistage core-formation model.

For N-isotopic fractionations among the atmosphere, silicate melt, and metallic melt at the ith stage of accretion, δ15N also follows the mass balance rule:

$${{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{silicate}}}}}}}times f{{{{{rm{N}}}}}}_{i}^{{{{{{rm{silicate}}}}}}}+{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{metal}}}}}}}times f{{{{{rm{N}}}}}}_{i}^{{{{{{rm{metal}}}}}}}+{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{atm}}}}}}}times f{{{{{rm{N}}}}}}_{i}^{{{{{{rm{atm}}}}}}}={{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{bulk}}}}}}}$$

(16)

where (f{{{{{rm{N}}}}}}_{i}^{{{{{{rm{silicate}}}}}}}), (f{{{{{rm{N}}}}}}_{i}^{{{{{{rm{metal}}}}}}}), and (f{{{{{rm{N}}}}}}_{i}^{{{{{{rm{atm}}}}}}}) are the N fractions in the silicate melt, metallic melt, and atmosphere, respectively; ({{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{bulk}}}}}}}) is the bulk N isotope composition of the silicate melt, metallic melt, and atmosphere. At the ith stage of accretion, ({{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{bulk}}}}}}}) can be expressed as:

$${{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{bulk}}}}}}}={{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i-1}^{{{{{{rm{silicate}}}}}}}{times} left({C}_{i-1}^{{{{{rm{N}}}}}-{{{{{rm{silicate}}}}}}}{times M}_{i-1}^{{{{{{rm{silicate}}}}}}}right)/\ left({C}_{i-1}^{{{{{{rm{N}}}}}}-{{{{{rm{silicate}}}}}}}{times M}_{i-1}^{{{{{{rm{silicate}}}}}}}+{C}_{0}^{{{{{{rm{N}}}}}}-{{{{{rm{imp}}}}}}}times {M}_{i}^{{{{{{rm{imp}}}}}}}+{M}_{i-1}^{{{{{{rm{N}}}}}}-{{{{{rm{atm}}}}}}}-{M}_{i}^{{{{{{rm{N}}}}}}-{{{{{rm{loss}}}}}}}right)\ +{{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{imp}}}}}}}times C}_{i}^{{{{{{rm{N}}}}}}-{{{{{rm{imp}}}}}}}times {M}_{i}^{{{{{{rm{imp}}}}}}}/Big({C}_{i-1}^{{{{{{rm{N}}}}}}-{{{{{rm{silicate}}}}}}}{times M}_{i-1}^{{{{{{rm{silicate}}}}}}}\ +{C}_{0}^{{{{{{rm{N}}}}}}-{{{{{rm{imp}}}}}}}times {M}_{i}^{{{{{{rm{imp}}}}}}}+{M}_{i-1}^{{{{{{rm{N}}}}}}-{{{{{rm{atm}}}}}}}-{M}_{i}^{{{{{{rm{N}}}}}}-{{{{{rm{loss}}}}}}}Big)\ +{{{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i-1}^{{{{{{rm{atm}}}}}}}times} left({M}_{i-1}^{{{{{{rm{N}}}}}}-{{{{{rm{atm}}}}}}}-{M}_{i}^{{{{{{rm{N}}}}}}-{{{{{rm{loss}}}}}}}right)/\ left({C}_{i-1}^{{{{{{rm{N}}}}}}-{{{{{rm{silicate}}}}}}}{times M}_{i-1}^{{{{{{rm{silicate}}}}}}}+{C}_{0}^{{{{{{rm{N}}}}}}-{{{{{rm{imp}}}}}}}times {M}_{i}^{{{{{{rm{imp}}}}}}}+{M}_{i-1}^{{{{{{rm{N}}}}}}-{atm}}-{M}_{i}^{{{{{{rm{N}}}}}}-{{{{{rm{loss}}}}}}}right)$$

(17)

where ({{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{imp}}}}}}}) is the N-isotopic composition of the accreting impactor. The N-isotopic fractionation between metallic melt and silicate melt is:

$${{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{metal}}}}}}}{-{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{silicate}}}}}}}={triangle }^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{metal}}}}}}-{{{{{rm{silicate}}}}}}}$$

(18)

({triangle }^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{metal}}}}}}-{{{{{rm{silicate}}}}}}}) is a multi-function of temperature and fO2, as summarized in Eq. (2) in the main text. The equilibrium N-isotopic fractionation between the atmosphere and silicate melt at high temperatures is negligible71, so the following equation was used:

$${{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{atm}}}}}}}{-{{{{{rm{delta }}}}}}}^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{silicate}}}}}}}={triangle }^{15}{{{{{rm{N}}}}}}_{i}^{{{{{{rm{atm}}}}}}-{{{{{rm{silicate}}}}}}}=0$$

(19)

The combination of Eqs. (16)–(19) can be used to constrain the N-isotopic composition in the atmosphere, silicate magma ocean, and core in the framework of a multistage core formation model.

### Test the sensitivity of varying the parameters used in the above N-accretion model

We recognized that during Earth’s multistage accretion, the processes such as collision-induced catastrophic atmosphere loss and the degree of silicate magma ocean–atmosphere equilibrium remain poorly constrained. We therefore performed additional modeling to test the effect of varying the parameters used in our model on the final N-content and δ15N of the proto-Earth’s atmosphere and mantle (Supplementary Figs. 8−14). Supplementary Fig. 8 shows that fixing Φ at 20–60% would result in N-content and δ15N of ~1.1 to 1.7 ppm and ~ −1.8‰ to −0.3‰ for the proto-atmosphere, and ~0.7 to 2.2 ppm and ~ −6.8‰ to −0.1‰ for the proto-mantle. These values cover the observed ones for Earth’s mantle and surface. Supplementary Fig. 9 shows that varying the degree of impact-induced atmospheric loss from 40% to 60% would vary the proto-atmospheric N-content from 0.7 to 1.5 ppm and δ15N from −3‰ to +4.2‰, but does not change the mantle values. These proto-atmospheric values also cover those observed values for Earth’s surface reservoir. In addition, Supplementary Fig. 10 shows that slightly increasing the degree of equilibrium between silicate magma ocean and overlying atmosphere or decreasing the added mass of CI chondrites would decrease the proto-atmosphere δ15N to ~0‰ and meanwhile increase the atmosphere N-content at ~1 ppm. Supplementary Fig. 11 shows that changing the relative timing of the delivery of CI chondrite-like materials would either not change the main conclusions, as long as the CI chondrite-like materials were delivered after Earth accreting its ~60% mass but before complete core–mantle segregation. Supplementary Fig. 12 shows that varying the C-content in the impactor core from 0.5 to 2 wt.% would either not change the conclusions. Supplementary Fig. 13 shows that varying the relative contributions of EC and CC to the Moon-forming giant impactor in δ15N would nearly not affect the δ15N of the proto-Earth’s atmosphere, mantle, and core. This is because most of the N in the last impactor was in the impactor core, and limited emulsification (5%) of the impactor core resulted in most of the N in the impactor delivered to Earth’s core. The N-content in the impactors must affect Earth’s bulk N-content. We compare the modeling results using impactors containing 80–200 ppm N (N-rich) and 20–100 ppm N (N-poor) in Supplementary Fig. 14. It shows that the final N-contents in the proto-atmosphere (1.8 vs. 1.4 ppm) and proto-mantle (2.9 vs. 2.1 ppm, respectively) are well comparable, while the N-content in the core is more different in these two cases (234 vs. 118 ppm), because N is siderophile during the whole accretion phase. Supplementary Fig. 14 also shows that the δ15N of the proto-atmosphere, proto-mantle, and core is nearly not changed in these two cases, and any slight change in δ15N can be counterbalanced by slightly increasing the mass of CI chondrites added.

### Constrain the distribution and origin of other major volatiles (C–H–S) during Earth’s heterogeneous accretion

The other major volatiles (C–H–S) would undergo similar accretion process as modeled above for N. Since these volatile elements have very different geochemical behaviors relative to N, we had to consider the C–H–S contents in the accreted materials, their solubilities in silicate melts at the atmosphere–mantle equilibrium conditions, and their metal/silicate partition coefficients at the core-formation conditions, before we applied our above N-accretion model to constrain the distribution and origin of C–H–S in Earth’s different reservoirs.

For the C–H–S contents in the accreted materials, we assumed that the completely oxidized materials have CI chondrite-like C–H–S contents, and the C–H–S contents in other impactors were inferred from the asteroid meteorites. Hydrogen is mainly concentrated in the silicate part of the impactors due to the lithophile nature of H under asteroid conditions where the interior pressure is <20 GPa. Since the H contents in the achondrites are in the range of 50–170 ppm78,79, in our model the bulk H contents were changed from ~50 ppm in the reduced impactors to ~150 ppm in the oxidized impactors. Carbon and S are more concentrated in the metal core of the impactors because of their siderophile nature under asteroids conditions80,81. The C contents in iron meteorites are no more than 1500 ppm82, and the C contents in achondrites are in the range of 10–140 ppm83,84. We hence set the C contents in the impactors at ~150–1000 ppm. The S contents vary greatly in iron meteorites (0.4–19 wt.%)82 and achondrites (110–5000 ppm)85,86. However, the log(C/S) of iron meteorites are well correlated with the log(C) (ref. 82). Hence, we set the S contents in the bulk impactors in the range of ~2000–11000 ppm.

The surface magma ocean fO2 changes from IW−3.6 to IW+3 during Earth’s accretion in our model, as stated above, which would significantly influence the species and hence the solubilities of C–H–S in silicate melts77. In silicate melt, the relative fraction of H existing as H2 or H2O follows a relationship as:87

$${{{{{{rm{H}}}}}}}_{2}({{{{{rm{wt}}}}}}.%)=0.0873times {{{{{{rm{H}}}}}}}_{2}{{{{{rm{O}}}}}}({{{{{rm{wt}}}}}}.%) – 0.0044times Delta {{{{{rm{IW}}}}}} – 0.0198$$

(20)

In our model, the H contents in the impactors are no more than 150 ppm (see above), Eq. (20) indicates that the fraction of H as H2 is less than 30% at fO2 of IW−3.6 and would be negligible if the fO2 is higher than IW−2.5. This means that the main H species in the surface magma ocean is H2O during the whole accretion process. We hence used the H2O solubility model from ref. 88 in our model:

$$2{ln}{X}^{{{{{{rm{H}}}}}}_{2}{{{{rm{O}}}}}-{{{{{rm{silicate}}}}}}}=frac{2565}{T}+sum {b}_{j}{X}_{j}frac{{P}_{{{{{rm{g{as}}}}}}}}{T}+1.171times {ln}{f{{{{{{rm{H}}}}}}}_{2}{{{{{rm{O}}}}}}}^{{{{{{rm{fluid}}}}}}}-14.21$$

(21)

where ({X}^{{{{{{rm{H}}}}}}_{2}{{{{{rm{O}}}}}}-{{{{{rm{silicate}}}}}}}) and ({f{{{{{{{rm{H}}}}}}}}_{2}{{{{{rm{O}}}}}}}^{{{{{{rm{fluid}}}}}}}) are the mole fraction of water in the silicate melt and the fugacity of water in the fluid, respectively. Pgas is the total pressure of the gas, bj and Xj are the constant factor and oxide mole fraction of component j in silicate melt (bAl2O3 = −1.997, bFeO(T) =  −0.9275, bNa2O = 2.376)56. Assuming that the fugacity coefficient of H2O is equal to 1, the ({f{{{{{{rm{H}}}}}}}_{2}{{{{{rm{O}}}}}}}^{{fluid}}) is equal to its partial pressure and could be calculated by using Eq. (15) with r(H2O) = 9. The Pgas could be obtained by summing the partial pressure of all volatiles in the atmosphere.

The species of C in H-bearing silicate melts are also controlled by fO2. At reduced conditions (fO2 < IW−1), the main C species is methane, while at oxidized conditions (fO2 >IW−1), the main C species is carbonate89. For C solubility in silicate melt (SC-silicate) at fO2 < IW−1, we used the model from ref. 89 as:

$${{log }},{S}^{C-{{{{{rm{silicate}}}}}}}({{{{{rm{ppm}}}}}})=0.96times {{{log }}X}^{{{{{{rm{H2O}}}}}}-{{{{{rm{silicate}}}}}}} – 0.25times Delta {{{{{rm{IW}}}}}}+2.83$$

(22)

At fO2 >IW−1, we used the model from ref. 90, which is more valid for peridotitic melt. For S solubility in silicate melt (SS-silicate), we used a recent model of ref. 77:

$${{{{{rm{ln}}}}}},{S}^{S-{{{{{rm{silicate}}}}}}}(ppm)=13.8426 – 26476/T({{{{{rm{K}}}}}})+0.124times {C}^{{{{{{rm{FeO}}}}}}}+0.5times {{{{{rm{ln}}}}}}({f{{{{{rm{S}}}}}}}_{2}/{f{{{{{rm{O}}}}}}}_{2})$$

(23)

where the CFeO is the FeO content in the silicate melt (in wt.%), fS2 is sulfur fugacity, which could be approximately regarded as partial pressure of S2 in the atmosphere, and then could be expressed in a form like Eq. (15).

As for the metal/silicate partition coefficients, we used the parameterized ({D}_{H(i)}^{{{{{{{rm{metal}}}}}}}/{{{{{{rm{silicate}}}}}}}}) from ref. 91, ({D}_{C(i)}^{{{{{{rm{metal}}}}}}/{{{{{rm{silicate}}}}}}}) from ref. 80, and ({D}_{S(i)}^{{{{{{rm{metal}}}}}}/{{{{{rm{silicate}}}}}}}) from ref. 43.

With the above parameters for C–H–S, we applied our N-accretion model to C–H–S to constrain the distribution of C–H–S in Earth’s different reservoirs during Earth’s heterogeneous accretion. Our results presented in Figs. 5, 6 in the main text show that our modeled C–H–S contents in Earth’s different reservoirs are consistent with the present-day observations. Our study suggests that Earth acquired its volatiles through very complicated processes, but these complicate processes happened naturally during Earth’s main accretion phase. Therefore, the establishment of Earth’s volatile inventory was a natural outcome of Earth’s heterogeneous accretion processes.