Earth’s core could be the largest terrestrial carbon reservoir

Evaluating carbon’s candidacy as a light element in the Earth’s core is critical to constrain the budget and planet-scale distribution of this life-essential element. Here we use first principles molecular dynamics simulations to estimate the density and compressional wave velocity of liquid iron-carbon alloys with ~4-9 wt.% carbon at 0-360 gigapascals and 4000-7000 kelvin. We find that for an iron-carbon binary system, ~1-4 wt.% carbon can explain seismological compressional wave velocities. However, this is incompatible with the ~5-7 wt.% carbon that we find is required to explain the core’s density deficit. When we consider a ternary system including iron, carbon and another light element combined with additional constraints from iron meteorites and the density discontinuity at the inner-core boundary, we find that a carbon content of the outer core of 0.3-2.0 wt.%, is able to satisfy both properties. This could make the outer core the largest reservoir of terrestrial carbon. A carbon content in Earth’s outer core between 0.3 and 2.0 % by weight, along with at least two other light elements, is compatible with observational constraints, according to molecular dynamics simulations, and could make the core Earth’s largest carbon reservoir.

E arth's core is chiefly iron (Fe) but must contain some proportion of lighter elements to have the density based on the seismological observation [1][2][3] . However, the identity and proportion of the light elements are largely unknown [4][5][6][7] . Carbon is one of the important potential light-element candidates in the Earth's core because of its chemical affinity to the metallic phase and its cosmochemical abundance 4,[8][9][10] . Carbon is highly soluble in metallic Fe-rich melt, with solubility reaching ≥5-8 wt.% at the core and core-forming magma ocean conditions 8,11,12 . Experiments that show such a high concentration of carbon in the Ferich melt are generally carbon saturated, whereas the carbon abundance in natural systems is expected to be well below the saturation limit imposed by graphite/diamond. To understand the behavior of carbon during large-scale melting during accretion and constrain its abundance in the mantle and core, the partition coefficient of carbon between Fe-rich alloy melt and silicate melt, D metal=silicate C has been estimated over the past decade. The high-pressure-temperature laboratory experiments show that D metal=silicate C is always >1, i.e., it behaves as a siderophile element under core-forming conditions [12][13][14][15][16][17][18][19][20][21][22][23][24] .
Although all studies point toward the siderophile behavior of carbon, there remains a debate on the extent of its siderophile behavior. Most studies, especially for alloy compositions that are poor in other light elements such as S, Si, and N determined D metal=silicate C values in the order of hundreds to thousands 13,[15][16][17][18]20 . Whereas some studies argued for D metal=silicate C values for systems containing S-poor Fe-rich alloy as low as~5-50 [23][24][25] . In a recent study, it is not only argued in favor of much lower D metal=silicate C at very deep magma ocean conditions, combining lower pressure data, but a strong negative pressure dependence of D metal=silicate C has also been reported 23 . Yet, based on experimental data to 15 GPa no such negative pressure effect has been observed. A weaker preference of C for Fe-rich metallic melt in the recent studies is at odds with the first-principles molecular dynamics (FPMD) simulation on the speciation of carbon in silicate melt with pyrolytic composition, which suggests that the significant amount of carbon can be partitioned into iron-rich metal based on the strong clustering between iron and carbon 26 . The uncertainty in the carbon budget of the Earth's outer core does not only stem from the uncertainty in the extent of siderophile behavior of C during magma ocean processes. It is also related to the poorly constrained concentration of carbon in inner Solar System planetary bodies during accretion as well as the extent of alloy-silicate equilibration during core formation. Because of these various sources of uncertainties, the estimated carbon content of the Earth's core based on models of accretion and alloy-silicate differentiation varies by about a factor of 30, i.e., from 0.1 to~3 wt.% 12,20,23,24,27 .Therefore, an independent estimate of the carbon content of the Earth's present-day outer core, based on other approaches, is highly desirable.
The carbon budget of the Earth's outer core has also been estimated using the partitioning behavior of carbon between liquid iron and condensed phase of iron, i.e., hcp Fe, which indicates that carbon preferentially partitions into liquid outer core 28 . Based on this strong partitioning of carbon in the liquid outer core, it is predicted that 1.5-2 wt.% carbon in the outer core can explain the density discontinuity of 4.5-6.5% at the inner and outer core boundary (ICB) if carbon was the only light element [28][29][30] . However, the density deficit in the Earth's inner core cannot be explained if no other light elements are present in the inner core. High-pressure and temperature-phase diagram of Fe-C alloy has also been used to constrain the carbon content of the Earth's core. Previous phase relation experiments on Fe-C system have shown that either carbon dissolved in hcp Fe as an interstitial defect or Fe 7 C 3 phase could crystallize at inner-core conditions 11,[31][32][33] . Based on the eutectic temperature in Fe-C phase diagram, Earth's outer core could contain <3 wt.% C and the inner core could contain <1 wt.% C 33 . The carbon content estimated using the phase diagram can partially account for the shear wave velocity anomaly, anisotropy, and Poisson's ratio of the inner core [34][35][36] . Yet, 1 wt.% C cannot adequately compensate for the density deficit in the inner core. In addition, the presence of other light elements could limit the amount of carbon in the Earth's core 37 . Based on the melting temperature of pure Fe, ICB temperature varies from 5500 to 6300 K 38,39 . Light-element impurities could reduce the melting temperature of Fe at core pressures 40 . Thus, constraints on the concentration of light elements can also be deduced from the reduction of the melting temperature. Previous experiments suggest that reduction of melting temperature due to sulfur (S) 41 , hydrogen (H) 42 could be as high as~1000 K at ICB. In contrast, the effect of oxygen (O) 43 and silicon (Si) 44 of reduction of melting temperature are limited to 100 K. Carbon could reduce the melting temperature of iron by 600-800 K 40,45 . However, in a multi-component system with several light elements in the outer core, it is challenging to ascertain the reduction of melting temperature and provide constraints of the light-element concentration solely based on the core temperature.
Independent constraints on the carbon content of the Earth's core can be obtained by examining the effects of carbon on the density (ρ) and compressional wave velocity (V P ) of the iron alloy at pressure (P) and temperature (T) conditions relevant for the Earth's outer. Owing to the challenges associated with experiments at the conditions relevant for the Earth's outer core, the experiments constraining the ρ and V P of Fe-C liquids are limited in pressure and temperature conditions (i.e., P < 10 GPa, T < 2500 K) [46][47][48][49][50][51] . Recently, in situ X-ray diffraction measurement have been performed at significantly higher P-T conditions of up to 58 GPa 52 . Similarly, the sound wave velocity of Fe 84 C 16 liquid was measured using inelastic X-ray scattering up to 70 GPa 53 . Yet, large extrapolations are still required to estimate the ρ and V P even at core-mantle boundary conditions, i.e., 136 GPa 52,53 . Since most of the experiments are performed at lower P-T conditions, an empirically established linear relationship between density and sound wave velocity, also known as Birch's law 54 , is usually utilized to extrapolate density and elastic properties at core conditions 55 . However, the validity of Birch's law for solid iron at very high P-T is debated 56 . Thus, the first-principles molecular dynamics (FPMD) simulations can provide valuable insight into the density and elastic properties of metallic liquids at extreme P-T conditions relevant to Earth's core [57][58][59][60][61][62] . Recent FPMD simulation of liquid Fe-Ni-C at temperatures of~1673 K and pressure ranging between 0 and 67 GPa demonstrated that the short-range liquid structure in FPMD simulation agrees well with the experiments at similar P-T conditions providing an essential benchmark for simulation 63,64 .
Previous FPMD simulations were used to provide constraints on possible light-element candidates for the Earth's outer core considering Si, O, S, and C elements 6 . More recent FPMD studies examined the core composition and also included hydrogen as a potential light-element candidate at outer core compositions 65,66 . However, these studies were performed only at core P-T conditions with limited data points. Because of sparse data, reliable thermal equations of state for a particular alloy system are not available. In addition, the empirical correction used for adjusting pressure also varies among studies, and nitrogen as a lightelement candidate was not considered in prior studies. Here, we explored Fe-C liquid binary systems with~4,~6, and~9 wt.% of dissolved carbon at a wide range of pressure and temperature conditions i.e., 0-360 GPa and 4000-7000 K, covering those that are relevant for Earth's core. We also provide the thermal EOS and explore the validity of Birch's law. In order to validate our simulation with high-pressure experiments 67,68 , we use thermodynamically calculated pressure corrections. By combining the effects of C on Fe-liquid alloy density and compressional velocity with similar effects of the other five possible light elements (Si, N, O, S, and H) in the core, we directly constrain the possible carbon budget of the Earth's present-day outer core. We combine our outer core C content estimates with the bulk silicate Earth (BSE) C estimates from literature and provide new estimates for the bulk Earth (BE) carbon abundance.

Results
Our results show that the pressure (P), volume (V), and temperature (T) relationship of Fe-C binary melt can be represented using the Mie-Grüneisen thermal EOS where P V; T ref À Á is the pressure at reference isotherm and dP dT À Á V is the thermal pressure. We used 7000 K as our reference isotherm and the pressure-volume (density) relation along 7000 K can be described using the third-order Birch Murnaghan equation of state ( Table 1). The thermal pressures dP dT À Á V are derived at constant volume or density using the linear relationship of pressure and temperatures along four isotherms, i.e., 4000-7000 K. Our result shows that the density of Fe-C liquid decreases with increasing carbon concentration along all isotherms ( Fig. 1 and Supplementary Fig. S1).
To test the predicted pressures for a given density of a melt, we also simulated Fe 75 C 21 (~5.7 wt.% carbon) at the constant temperature of 2273 K and compared with experimental melt density from X-ray absorption experiment 47 . We found that the predicted pressure from FPMD simulations is smaller than the experiments bỹ 8 GPa. So, we make an empirical positive pressure correction of 8 GPa to the calculated pressures. Corrected pressure-density data can be well described by a third-order Birch Murnaghan equation of state with ρ 0 ¼ 6:63 ± 0:04 gcm À3 , K 0 ¼ 65:3 ± 12:8 GPa, K 0 0 ¼ 8:4 ± 2:4 ( Table 1). Our zero pressure density (ρ 0 ) and the bulk modulus (K 0 ) are comparable with experimentally observed values from previous ultrasonic studies 48,50,51 and FPMD study of liquid Fe-Ni-C 64 ( Table 1). In recent FPMD studies, the slight discrepancy between the pressure predicted from FPMD simulation and that of the experiments was addressed by adding constant empirical pressure correction to the FPMD results 6,57,60,64,69 . However, the recent computational study showed that the pressure difference between FPMD and experiment could vary both as a function of pressure and temperature 66 . Thus, a constant pressure correction for the entire core P-T conditions may not be sufficient. In our study, we have adopted a pressure correction constrained by thermodynamics (P th corr ), which was implemented in previous FPMD studies 67,68 (Supplementary Note S1 and Supplementary Figs. S2-S4).
We estimate the Grüneisen parameter (γ) for carbon-iron alloys using the relation: Our results indicate that the thermal Grüneisen parameter (γ) for iron-carbon alloy are~1.0 at ambient pressure and increases upon further compression up to pressures of~50 GPa. At higher pressures, γ uniformly decreases with increasing pressure (Supplementary Fig. S1). Thus, the value of γ continuously decreases from core-mantle boundary (CMB) to inner-core boundary (ICB). Comparison of γ for carbon-bearing metallic melt with that of pure iron form previous study shows that the carbon reduces the value of γ ( Supplementary Fig. S1). The reduction of γ due to carbon is consistent with the effect of other light elements in previous studies of metallic melts, i.e., Fe-N and Fe-H liquids 59,61 . The pressure evolution of Fe-Fe coordination for all melt compositions shows an increase with increasing pressure up to~50 GPa and beyond 50 GPa, Fe-Fe coordination remains mostly unaffected upon further compression (Supplementary Fig. S1). At high pressures, Fe-Fe coordination for pure iron is 13, indicating the optimal packing Fe atoms similar to cubic closest packing (ccp) or hexagonal closest packing (hcp) of solid iron 59 . With increasing carbon concentration, Fe-Fe coordination decreases, from 12 to 10 for carbon concentration from 4 to 9 wt.% ( Supplementary Fig. S1). The short-range atomistic scale structure of iron with 4 wt.% C also exhibits a nearly close-packed structure. However, the local network of iron atoms becomes less close-packed with increasing carbon concentration. The structural change manifests itself in the reduction of the density of Table 1 Density (ρ 0 ), bulk modulus (K 0 ), and the pressure derivatives of bulk modulus (K 0 0 ) for iron-carbon alloys at room pressure and various temperature (T) conditions. iron-carbon melts compared to that of the volatile free molten iron ( Supplementary Fig. S1). P-wave velocity (V P ) of liquid iron-carbon can be calculated from where K S is the adiabatic bulk modulus, ρ is the density, G is the shear modulus. For the liquids where the frequency of the seismic probe is much smaller than the characteristic frequency of shear relaxation, shear modulus G ! 0 70 . The adiabatic bulk modulus (K s ) is related to the isothermal bulk modulus ðK T Þ by and, Using both the γ and dP dT À Á V from FPMD simulations, we estimate the V P for iron-carbon liquids at core conditions. Our result shows that the light-element carbon increases the V P of liquid bỹ 94 ms −1 for each wt.% of dissolved carbon at core conditions ( Figs. 2 and 3).
We estimate the ρ and V P of Fe-C alloys at core conditions along the outer core geotherm. We assume that the temperature gradient in the outer core is adiabatic, the outer core temperature wt.% of dissolved carbon as a function of pressure. Symbols are pressures estimated using first-principles molecular dynamics (FPMD) simulations at constant density, at 7000 K (red), 6000 K (blue), 5000 K (brown), and 4000 K (green). The pressure-density data from FPMD can be adequately described by a thermal equation of state (Table 1). Panel (b) shows the density of pure liquid iron, B19 59 and liquid Fe-4 wt.% C (this study) calculated along the geotherm with inner-core temperature (T ICB = 6000 K) ( Supplementary Fig. S2). Filled squares represent the density of liquid iron alloy with a stoichiometry of F 90 C 18 (~4 wt.% C) from FPMD simulations at core conditions 6 . Experimental results of the density of liquid iron alloy with a stoichiometry of F 84 C 16 (~4 wt.% C) liquid up to 200 GPa along geotherm (with T CMB ¼ 4300 K) is also shown (cyan line) for comparison 53   Integrating Eq. (4) and assuming a fixed temperature ðT ICB Þ at the inner-core boundary (ICB) and adopting the density at ICB ðρ ICB Þ from the preliminary reference Earth model (PREM) model, we derive an adiabatic temperature profile for the Earth's outer core 71 . Temperatures at ICB conditions are estimated to vary between 5500 and 6300 K based on the melting temperature of iron 38,39,72 . If the ICB temperature is held constant at 6000 K, based on the above approximation, the temperature at the core-mantle boundary (CMB) along the core adiabat is found to be~4400 K. Similarly, if the ICB temperature is held at lower temperatures of 5400 and 5000 K, temperature at the core-mantle boundary (CMB) along the core adiabat are found to be~4000 and~3700 K, respectively ( Supplementary Fig. S5). Our estimated outer core adiabat is in good agreement with the previous estimation based on thermodynamic parameters 57,59 . The temperature at the ICB will likely be lower if light elements are present. Previous work on Fe-C alloy estimated the melting temperature of Fe-C to be 5500 ± 500 K by considering the melting point depression of~600-800 K 40,45 .
For condensed matter phases such as solid minerals and alloys, the sound velocity often increases linearly with the density and is often referred to as Birch's law 54 . Birch's law is often used to extrapolate the elastic properties of the material from the lowpressure conditions to the core pressures. We tested the validity of Birch's law by comparing the bulk sound velocities of liquid iron and iron-carbon liquid along isotherms. Our results show that the isobaric temperature dependence of V P for metallic melt, i.e., dV P dT P is negligible. However, V P varies as a function of temperature at a constant density, i.e., V P increases with increasing temperature at a fixed density (Fig. 2 and Supplementary Fig. S5). Our observations are similar to previous FPMD simulations on Fe-bearing metallic liquids, which found similar dV P dT P for liquid Fe, Fe-S, Fe-N, and Fe-H liquid 57,59,61,73 . This observation indicates that Birch's law may not hold for pure liquid iron as well as other metallic liquids such as Fe-C, Fe-S, Fe-N, and Fe-H at core pressure and temperature conditions. Our result also agrees with the previous experiment on solid iron that indicated a nonlinear relationship between density and sound wave velocity for hexagonal closed packed iron up to 73 GPa and 1700 K 56 . However, other studies on solid and liquid metals do not agree with ref. 56 . For example, density and sound wave velocity measurement for solid Fe up to 93 GPa and 1100 K shows that the density and sound wave velocity data do follow Birch's law 74 . Also, the shock-wave experiment of liquid Fe and liquid Fe-O-S do support Birch's law at much higher pressure up to 800 GPa 75,76 .

Discussion
Our calculations allow us to estimate the ρ and V P of liquid Fe-C at core pressure conditions along the adiabatic geothermal gradient. The adiabatic geothermal gradient was estimated by using a reference ICB temperature at 6000 K. Because the exact ICB temperature is unknown and there is a range of estimates, to account for any differences caused by the ICB reference temperature, we also evaluated the melt properties with reference ICB temperature at 5400 and 5000 K. To predict the carbon content of the Earth's outer core, we compare the ρ and V P of iron-carbon liquid with the preliminary reference Earth model (PREM) values at CMB and ICB conditions 71 . We find a linear relationship for both ρ and V P of Fe-C liquid alloy and the carbon content (X = wt.% of carbon), with dρ dX C < 0 and dV P dX C > 0 (Fig. 3).
For the core adiabat with T ICB = 6000 K, we find that PREM density at core-mantle boundary (i.e., T CMB = 4400 K) could be explained by 6.9 ± 0.4 wt.% of carbon dissolved in liquid iron alloy (Fig. 3). However, the estimated V P for this amount of carbon exceeds seismological V P at CMB. Our estimates on the carbon content of liquid iron alloy at CMB with T CMB = 4400 K by matching V P is~1.3 ± 0.3 wt.% (Fig. 3). Similarly, at the innercore outer core (ICB) boundary with T ICB = 6000 K, we estimate the carbon content of 5.4 ± 0.2 wt.% and 1.3 ± 0.2 wt.%, by, respectively, matching ρ and V P of Fe-C liquid alloy and PREM (Fig. 3). If we consider lower core temperatures, i.e., along adiabatic geothermal gradient with T ICB = 5400 K,~7.5 and~5.8 wt. % carbon is necessary to explain ρ at CMB (T CMB = 4000 K) and ICB (T ICB = 5400 K) conditions, respectively (Supplementary  Table S1). Similarly, for the temperature profile with T ICB = 5000 K, a higher carbon amount, i.e.,~7.9 and~6.1 wt.% is necessary to explain ρ at CMB (T CMB = 3700 K) and ICB (T ICB = 5400 K). We find that the wt.% of carbon is less sensitive to the core Fig. 3 Effect of carbon on the density (ρ) and compressional wave velocity ðV P Þ of liquid iron. Change in (a) ρ and (b) V P of liquid Fe as a function of carbon content at core-mantle boundary (CMB: 136 GPa, 4400 K) and inner-core boundary (ICB: 330 GPa, 6000 K) conditions. Results for the metallic melts are from previous experimental studies, AA94, BM86, K20 [106][107][108] and FPMD studies, B14, I14 6,58 . Lines show the linear relationship of ρ and V P with the increasing amount of carbon (in wt.%). Filled gray squares represent the results after applying the thermodynamic pressure correction (P th corr ) (Supplementary Note S1). When the pressure is adjusted using P th corr , the wt.% of carbon required to satisfy the seismic ρ and V P becomes lower and higher, respectively. Error bars represent ±1σ uncertainties.
temperatures. For instance, the estimated wt.% of carbon required for explaining V P at lower temperatures, i.e., T ICB = 5400-5000 K and T CMB = 4000-3700 K are same as that estimated for higher temperatures of T ICB = 6000 K and T CMB = 4400 K. This is due to the fact that at constant pressures and in the explored temperature range relevant for the Earth's outer core, i.e., between 3700 and 7000 K, the dV P dT is negligible (Supplementary Fig. S5). Effect of pressure correction is such that the estimated carbon content to match the seismic density decreases and the carbon content needed to match the velocity increases (Supplementary Note S1 and Fig. 3). For the geotherm with T ICB = 6000 K and pressure corrected using thermodynamic constraints (P th corr ), Fe with~4.9 wt.% carbon will satisfy PREM density at ICB condition, whereas liquid iron with~5.6 wt.% carbon can reproduce PREM density at CMB pressure and 4400 K (Supplementary  Table S1). However, Fe with~2.1 wt.% and~3.9 wt.% carbon is required to match the V P at ICB and CMB conditions, respectively.
Our results at conditions relevant to the Earth's outer core indicate that Fe-C binary liquid alloy cannot simultaneously satisfy the ρ and V P of the outer core. High-pressure studies based on in situ inelastic X-ray scattering suggested a 2.9-3.8 wt.% of carbon could explain the density deficit of the core 53 . However, 2.9-3.8 wt.% of carbon would translate to a V P that is greater than the observed core V P . In fact, a lower concentration of carbon, i.e.,~0.9-1.2 wt.% was required based on the comparison between V P of the Fe-C alloy and the Earth's outer core 53 . More recently, in situ X-ray diffraction studies using a Paris-Edinburgh press and laser-heated diamond anvil cell of Fe-C liquid suggested that Fe with 1.8-3.7 wt.% carbon was sufficient to satisfy density profile at the core conditions 52 . Our results are in good agreement with an earlier FPMD simulations on Fe-C binary, which also found that Fe-C binary liquid alloy cannot simultaneously satisfy the ρ and V P of the outer core, and around~5 wt. % and~2 wt.% of carbon is required satisfy core ρ and V P , respectively 6 (Fig. 3). However, the actual carbon concentration in the core is expected to be lower than these estimates based on the study of Fe-C binary alloy due to the presence of other light elements.
Using the ratio of major and trace elements and assuming lowpressure core accretion under reducing conditions, it is estimated that the Earth's core may contain 7.35 wt.% Si, 2.30 wt.% S, and 4.10 wt.% O 5 . More recent core-formation models proposed that the Earth's core may contain 2 wt.% Si for a relatively oxidizing conditions 77 and 8-9 wt.% Si for reducing conditions 78 . Previous FPMD work indicated that 1.9 wt.% Si together with 3.7 wt.% O is the best match for both ρ and V P of outer core 6 . In follow-up work, it was argued that an outer core composition that contains 2.7-5 wt.% oxygen and 2-3.6 wt.% silicon can also satisfy the geochemical constraints of core formation, i.e., such coreformation scenario can also result in a mantle with the observed abundance of key siderophile elements, such as Ni, Co, V, and Cr. Based on partition coefficients, sulfur-rich core with 8-10 wt.% S is suggested 79,80 . However, geophysical constraints indicate sulfur poor core with 2-3 wt.% S 81 . It was also concluded that the maximum sulfur in the Earth's outer core is <2.4 wt.% and it might be necessary to have oxygen in all the proposed compositional models for the outer core 6 . Shock-wave experiments on Fe-Si-O liquid indicated that the Earth's outer core may be oxygen poor with less than 2.5 wt.% O 76 . More recent FPMD work of Fe-S binary alloy set the upper limit of sulfur to 3.5 wt.% 62 .
Since, the Earth's core likely contains other light elements such as oxygen, sulfur, silicon, hydrogen, and nitrogen in addition to carbon, we refine our estimates of the carbon content in the Earth's outer core based on the possibility of simultaneously explaining the ρ and V P of core by considering a ternary mixture of Fe-C-LE where C, and LE refers to carbon and the other light elements, respectively- ρ is plotted on the left axes (purple) and V P on the right axes (blue) in all panels. Pressure is corrected using thermodynamic constraints, P th corr (Supplementary Note S1). Thick gray line in each panel shows the density and compressional wave velocity at core conditions from preliminary reference earth model 71 . We note that the dρ dX LE and dV P dX LE for Si and N are very similar to that of carbon and hence it difficult to obtain a unique solution with Fe-Si-C and Fe-N-C ternary systems. Green bar in panels (a-f) indicates the unique solution in Fe-C-LE ternary that simultaneously explain ρ and V P at ICB or CMB conditions. No unique solution is found for (g, h) Fe-C-Si and (i, j) Fe-C-N ternary systems. Error bars represent ±1σ uncertainties.
where, X C and X LE are weight fraction of carbon (C) and other light elements (LE); Δρ refers to ρ Core À ρ Fe À Á and ΔV P refers to V Core P À V Fe P À Á . We assume cross composition terms i.e., d 2 ρ dCdX LE and d 2 V P dCdX LE are essentially negligible or "zero". This assumption has been successfully applied in prior studies 6,66,76 . We estimate the effect of light elements-O, S, Si, H, N on ρ and V P of liquid iron i.e., dρ dX LE and dV P dX LE calculated from prior studies 6,59,62,65,66 . We estimate the effect of carbon on ρ and V P of liquid iron i.e., dρ dX C and dV P dX C from this study (Fig. 3). In our analysis, we have not considered the effect of nickel but we do not anticipate a different result even if~5-10 wt.% Ni is present in the outer core because previous FPMD studies showed that the effect of Ni on the ρ and V P of liquid iron are small 6,65,66 .
Our analysis shows that we are unable to simultaneously explain ρ and V P for Earth's entire outer core by considering any Fe-C-LE ternary. However, we are able to explain ρ and V P , simultaneously and provide a unique solution for C either at the CMB or the ICB conditions, with LE as O, S, and H but such solution is also not obtained if carbon is considered with Si, and N as light elements (Fig. 4). For the core adiabat with T CMB = 6000 K, we are able to simultaneously explain ρ and V P for Earth's outer core at CMB conditions with (a)~2.7 wt.% C with 3.1 wt.% O for a Fe-C-O ternary; (b)~2.5 wt.% C with~5.6 wt. % S for a Fe-C-S ternary; and (c)~2.5 wt.% C with~0.5 wt.% H for a Fe-C-H ternary. Similarly, we can simultaneously explain ρ and V P for Earth's outer core along core adiabat of T ICB = 6000 K at ICB conditions with (a)~0.3 wt.% C with~6.2 wt.% O for a Fe-C-O ternary and (b)~0 wt.% C with~11.1 wt.% S for a Fe-C-S ternary. When the lower core adiabat with T ICB = 5400 K and 5000 K is considered, velocity and density can be explained with a lower concentration of carbon, which is compensated by a higher concentration of S and H. For example, with T CMB = 3700 K,~1.6 wt.% carbon with~9 wt.% S or~0.8 wt.% H can simultaneously explain ρ and V P at CMB. These analyses indicate that if carbon is a component in any model ternary outer core system, the closest solution is its concentration varying between 0 and 2.7 wt.% (Fig. 4). Our finding that the higher concentration of sulfur and hydrogen is favored at lower core temperatures is consistent with the large melting point depression of liquid Fe by S and H compared to that by C [40][41][42]45 . We find that carbon concentration between 0.3 and 2.7 wt.% can simultaneously explain ρ and V P in the Earth's outer core when oxygen is considered as another light element. However, when we consider several light elements, LE ¼ H, S, Si, and N in the Fe-C-LE ternary systems, our results show that the lower limit of carbon in the outer core could be as low as~0 wt.%. However, parental alloy melt for magmatic iron meteorites can contain up to~1100 ppm C [82][83][84] . Given iron meteorite parent bodies are the earliest formed protoplanets in the Solar System and have been argued to contribute volatile budget to later formed rocky embryos and planets 85 , it is unlikely that the core of large terrestrial planets such as Earth would not sequester some carbon during the accretionary process. Combining the evidence from iron meteorites and our analysis from Fe-C-O ternary, we place the lower limit of carbon in the Earth's outer core to 0.3 wt. %. Additional constraints on the upper limit of carbon content by comparing the density contrast between the outer and inner core with the density difference between solid iron 86,87 and liquid Fe-C alloys from our study suggest that the outer core should not have any more than~2 wt.% carbon (Supplementary Fig. S6) 29,30 . This is based on the assumption that carbon is the only a light element in the core and most of the carbon partitions into the liquid outer core. Our result on the ICB density discontinuity is consistent with recent constraints 28 . Thus, considering all of the above, we predict the expected range of carbon in the Earth's outer core to be as low as 0.3 wt.% and as high as 2.0 wt.%. Our predicted carbon concentration is consistent with the partitioning of carbon between silicate and metallic melts, which together with different accretional models and various assumptions of the extent of metal-silicate equilibration, provide a similar range of carbon content for the outer core~0.1 wt.% to~2 wt.% 12,20,23,27 . Our estimation is also similar to~1.5-2.0 wt.% carbon proposed for the Earth's outer core by matching the density discontinuity at ICB 28 . Estimation based on the eutectic temperature in Fe-C phase diagram is less than 3 wt.% 33 . The outer core with >3 wt.% C could crystallize Fe 7 C 3 phase at ICB 33 , and their densities are too small to compensate for the inner-core density deficit 88 . If the higher end of the outer core C content as constrained by our study and that of earlier studies is applicable, then the effective partitioning of C to the core during metal-silicate separation is expected to be those predicted by most C partitioning studies that produced partition coefficients in the range of hundreds or higher 14,15,17 .
Our study suggests that the Earth's outer core contains 0.3-2.0 wt.% C, i.e., 5:5 À 36:8 10 24 g of carbon, which is similar to some of the previous estimates 15,20,89 . There remains considerable uncertainty in the bulk Earth carbon budget, with recent estimates varying between~220 and 1100 ppm 7,23,90 . Assuming the inner-core carbon concentration to be negligible based on the previous estimation of strong preference of C toward liquid outer core relative to solid inner core 28 , and taking into account the recent estimates of the BSE carbon budget, i.e., 100-480 ppm 91,92 , we can re-assess the bulk Earth carbon (BE) budget. Summing up the BSE and the core carbon budget, with the constraint for the latter from the present study, the BE carbon concentration could be as high as~990-6460 ppm. Hence, the core could be the largest reservoir for carbon on Earth, accounting for as much as 93-95% of the total terrestrial carbon. Therefore, if the upper limit of core C suggested by our study is relevant, Earth's carbon budget maybe distinctly larger than estimated in various recent studies that suggest carbon acted only as a mild siderophile element during core formation. Although our study places a nonzero lower bound of 0.3 wt.% for the outer core C, constraining a strict lower bound will require future investigations likely combining mineral physics, cosmochemistry, geochemistry, and dynamics of the formation of Earth-like rocky planets.

Methods
We used FPMD in canonical (NVT) ensemble as implemented in Vienna ab initio simulation package (VASP) [93][94][95] . The calculation of forces and energies is based on density functional theory (DFT) with projector-augmented wave (PAW) method 96 . We used the generalized gradient-approximation formulation (GGA-PBE) 97 which has been successfully used in the other metallic systems 25,58,59,98 . We used Γ-point sampling to integrate the Brillouin zone with the time step of 1 femtosecond, fs (where 1 fs = 10 −15 s). We used a Nosé thermostat algorithm for constant temperature in our MD simulations 99 . We used an energy cutoff of 400 eV for the plane-wave basis set. To account the effect of limited cutoff energy, we added a volume-dependent Pulay correction (P Pulay ) to the pressure 100 . The P Pulay increases with increasing pressure and for the explored pressure range, i.e., 0-360 GPa, P Pulay varies from 0.5 to 1.5 GPa. To constrain the carbon content in the conditions relevant for the Earth's core, our calculations were performed in nonmagnetic liquid alloys. Previous experiments indicated that liquid Fe could exhibit magnetic moments at lower pressures of~1 bar 101 . Recent FPMD simulation on liquid iron and iron alloys (Fe-Si-S) indicated that magnetic moments in liquid were sustained at higher pressures 102 . It is well known that in crystalline iron alloys, magnetism influences elastic properties. At low pressures, crystalline iron alloys are often magnetic, and they undergo pressure-induced loss of magnetism i.e., magnetic to nonmagnetic phase, and concomitant with that there is stiffening of elasticity including bulk modulus 103,104 . To test the effect of magnetic moment on energetics of iron (Fe) and iron-rich alloys (Fe-C), we performed FPMD calculations with magnetic moments at core conditions. We set a magnetic moment of μ m ¼ 4μ B for each Fe atom in the liquid Fe and liquid Fe-C alloy. We found that the magnetic moments decrease to~0 within 1 ps of simulation time. At core pressures and temperatures, the magnetic moments do not have any effect on the predicted pressure ( Supplementary Fig. S7).
To explore how dissolved carbon affects the properties of liquid iron, we explored three different iron-carbon liquids, with a binary stoichiometry of Fe 66 C 30 (~9 wt.% C), Fe 74 C 22 (~6 wt.% C), and Fe 81 C 15 (~4 wt.% C). A total of 96 atoms are simulated in a cubic box with periodic boundary conditions. We consider several constant volumes (densities) along four different isotherms (4000-7000 K). In addition to the high temperatures explored for assessing the properties at core conditions, we also simulated Fe 75 C 21 (~5.7 wt.% carbon) at a temperature of~2273 K. This is to assess the difference between the FPMD predicted pressures and the experimental pressures for the same melt density 47 . The densities explored in this study, corresponding to the pressure range of 0-360 GPa, i.e., covering the entire range of core pressures. For each fixed volume, we conducted the FPMD simulation at 7000 K for 10-20 ps (1 ps = 10 −12 s). Temperature is then lowered to 6000 K within a shorter time of~4 ps. We use these equilibrated cells and performed MD simulation for longer time scales ranging between 10 and 45 ps, and this procedure is repeated for lower isotherms of 5000 K and 4000 K. Similar procedure is implemented for all iron-carbon liquids. Examination of the radial distribution function (RDF) and mean square displacement (MSD) indicate the simulated structure is in a liquid state. RDF for each atomic pair exhibits a short to mid-range order and long-range disorder with a strong first peak followed by a minima and fluctuation of RDF close to unity at a larger distance (Supplementary Fig. S8). Log-log plot of MSD vs time shows a slope of unity and the MSD exceeded 10 Å 2 within our simulation time indicating that the melt exhibits ergodic behavior ( Supplementary Fig. S9). We use the blocking method to determine the average energies and pressures for each simulation volume 105 . While taking an average, we discard the initial 10% of the data from the time series for the equilibration of the system. The time evolution of energy and pressure differences from the values averaged over the full simulation time shows that the results converge in a few picoseconds (ps) (Supplementary Fig. S9). Within~4 ps of simulation time, the difference in pressure and energy decreases to~0.5 GPa and~0.01 eV per atom, respectively ( Supplementary Fig. S9). We also explore the effect of pressure correction in calculated properties by considering a thermodynamically calculated empirical pressure correction (P th corr ) (Supplementary Note S1 and Supplementary Figs. S2-S4).

Data availability
All the data used in this study are presented in the main paper and in the supplementary document. Parameters used in the simulation are described in the "Methods" section. Raw data used to produce Figs. 1-4 and Supplementary Figs. S1-S9 can be accessed from the public repository Zenodo: https://doi.org/10.5281/zenodo.4903530. The repository also contains the files with simulation input parameters as well as the positions files to compute the coordination statistics ( Supplementary Fig. S1).