The density of the outer core is ~10% lower than that of iron-nickel alloy, and thus, the outer core must contain lighter elements in order to account for the density deficit (Birch, 1952(Birch, , 1964)). It is also possible that the light element abundance of the outer core gradually increases with inner core crystallization, assuming that most of the light elements in question are incompatible with solid Fe-Ni alloy (Buffett et al., 2000;Jacobs, 1992). However, the exact partitioning of the light elements between the inner and outer core is uncertain given that the inner core also is thought to contain ~2-4% light elements (Alfè et al., 2007;Anderson & Isaak, 2002;Hemley & Mao, 2001;Poirier, 1994;Vočadlo, 2007). The convection in the liquid outer core is one of the prime drivers for the core geodynamo, and the presence of light elements can assist in the chemical convection and sustenance of the geodynamo (Badro et al., 2016;Buffett et al., 2000;Lister & Buffett, 1995;Pozzo et al., 2013). Because it is evident that the light elements play a crucial role in the dynamics of the core, it is critical to put constraint on the concentration and the chemistry of the light elements and how it affects the structure, phase relations, elasticity, and transport properties of liquid and solid iron alloys at high pressures and temperatures (Badro et al., 2014(Badro et al., , 2016;;Buffett & Seagle, 2010;Helffrich, 2014;Huang et al., 2011). Owing to the large mass fraction of the Earth's core, even with a small concentration of a given light element, core may become a significant reservoir of that element in the bulk Earth (Dasgupta, 2013;Dasgupta & Hirschmann, 2010;Kaminsky & Wirth, 2017;Marty, 2012). Several candidates including H, C, O, Si, and S have received attention as possible light elements in the liquid outer and solid inner core (Birch, 1964;Chen et al., 2014;Hirose et al., 2013;Li & Fei, 2003;Poirier, 1994;Ringwood, 1977;Stixrude et al., 1997). Nitrogen is severely depleted in the bulk silicate Earth (BSE), and hence, it is also a possible candidate for the Earth's core (Marty, 2012;Marty & Dauphas, 2003). The consideration for nitrogen in the core is supported by the existing experimental data that suggest that nitrogen is a siderophile element in core-forming magma ocean conditions (e.g., Dalou et al., 2017;Kadik et al., 2011;Li, Dasgupta, et al., 2016;Roskosz et al., 2013). Recent estimates on nitrogen solubility in metallic melt predicts that several wt.% of nitrogen could be dissolved at high pressurestemperatures relevant for terrestrial magma oceans (Speelmanns et al., 2018). In fact, the depletion of nitrogen in BSE, relative to carbon, was used by Marty (2012) to argue that such superchondritic C/N ratio in BSE is a result of preferential partitioning of nitrogen to the core. More recent studies (e.g., Chi et al., 2014;Dalou et al., 2017), however, challenged such hypothesis given that the experimental alloy-silicate partition coefficient of carbon is Supporting Information:
• Supporting Information S1
• Table S1 • Table S2 • Data S1 significantly greater than that for nitrogen. Furthermore, nitrogen budget of the core is not only influenced by the alloy-silicate partition coefficient of N, D alloy=silicate N
, but also the available bulk nitrogen that participated in core-mantle fractionation relevant for terrestrial accretion. Owing to volatility during high-temperature accretional processes (e.g., Bergin et al., 2015;Trigo-Rodríguez, 2013), the budget of nitrogen available for sequestration in the core may have been limited. Similarly, the delivery of nitrogen to the Earth might have postdated core formation as has been suggested for many other volatile elements (Albarede, 2009;Dasgupta et al., 2013;Wang & Becker, 2013). Such late addition of nitrogen is not expected to influence the N budget of the core. Hence, arguments both against and in favor of nitrogen sequestration in the core exist from various geochemical considerations. Despite the geochemical arguments against significant concentration of nitrogen in the core, nitrogen cannot be ruled out as an important light element in the Earth's core on the basis of cosmochemical and early accretion process alone. For example, nitrogen may be delivered during core formation in the form of relatively refractory nitrides rather than less stable organics (Rubin & Choi, 2009), minimizing the volatility-induced loss. Moreover, a fraction of terrestrial core growth might have happened through near-disequilibrium merger of a differentiated planetary embryo (e.g., Li, Marty, et al., 2016;Rudge et al., 2010), in which case merger of a nitride-rich metallic core of more exotic composition formed from a different bulk composition could take place. The presence of nitrogen during core-mantle differentiation of Earth is also potentially supported by recent findings of nitride and carbonitride inclusions in lower mantle diamond (Kaminsky & Wirth, 2017) if lower mantle domains preserve memory of inefficient core formation. Furthermore, although N is shown to be moderate to mildly siderophile or even lithophile,
is known only at shallow magma ocean conditions and it remains an open question whether N becomes much more siderophile at deep magma ocean conditions. Hence, it would be important to place an independent constraint on nitrogen content of the Earth's core.
Estimates on the nitrogen in the liquid outer core could be obtained by examining its effect on the structure, density, elasticity, and bulk sound velocity and comparing such effects with the effects of other light elements such as Si, O, S, C, and H from literature (silicon: Badro et al., 2015;Sanloup et al., 2004;Yu & Secco, 2008;oxygen: Badro et al., 2014;Huang et al., 2011;sulfur: Jing et al., 2014;Sanloup et al., 2000;Umemoto et al., 2014;carbon: Kuwabara et al., 2016;Liu et al., 2016;Morard et al., 2017;Nakajima et al., 2015;Shimoyama et al., 2016;hydrogen: Umemoto & Hirose, 2015). Previous experimental studies showed that the solid iron-nitrides could possibly explain the geophysically observed density and sound velocities at the Earth's inner core (Adler & Williams, 2005;Litasov et al., 2017;Minobe et al., 2015;Niewa et al., 2009;Popov et al., 2015;Schwarz et al., 2009). Yet the effect of nitrogen on the physical properties of liquid iron in the Earth's core condition remains unknown. Hence, in this study we explore the effects of dissolved nitrogen on physical properties of liquid Fe-alloy at high pressures and temperatures to evaluate the possibility of N being an important light element in the Earth's outer core. In particular, using first-principles molecular dynamic (FPMD) simulation, we examine the effect of nitrogen on density, equation of state, and bulk sound velocity of iron alloy melt at the Earth's core conditions. Using such data, we evaluate the relative merit of nitrogen versus other light elements such as C, H, Si, O, and S as the dominant alloying component in Fe-rich core.
We used FPMD simulation with projector augmented-wave (Kresse & Joubert, 1999) method as implemented in Vienna ab initio simulation package (Kresse & Furthmüller, 1996a, 1996b;Kresse & Hafner, 1993). Our calculations are at Γ-point and with a plane wave energy cutoff of 400 eV. All FPMD simulations are based on the NVT canonical ensemble with a fixed number of atoms (N), constant volume (V), and a constant temperature (T). We used Nosé thermostat to maintain the constant temperature (Nosé, 1984). We used a time step of 0.5 fs (where, 1 fs = 10 À15 s) in all the FPMD simulations. We estimate the electronic exchange-correlation energy using the generalized gradient approximation (GGA) formulated by Perdew-Burke-Ernzerhof (Perdew et al., 1996). Previous density functional theory-based simulations showed that GGA method is more accurate compared to local density approximation for the description of ground-state electronic properties of iron (Alfè et al., 2000;Stixrude et al., 1994;Vočadlo et al., 2000). GGA-Perdew-Burke-Ernzerhof pseudopotential was successfully implemented in the FPMD simulations of crystalline iron up to 15 Mbar pressure and at 10.1029/2018GL080555 Geophysical Research Letters 7000 K (Godwal et al., 2015). In addition, density functional theory-based studies with GGA methods gave excellent agreement with experimental data for iron-alloy liquids (Ichikawa et al., 2014;Umemoto & Hirose, 2015;Zhang & Yin, 2012). In this study, we have used both local density approximation and GGA at a reference isotherm of 6000 K. We found that pressure determined by GGA method (P GGA ) for the pure iron liquid is in better agreement with the previous shock wave experiments (Anderson & Aherns, 1994;Brown & McQueen, 1986). In order to explore the effect of nitrogen on molten iron, we consider four distinct compositions that include pure iron (Fe), iron alloy with ~7, ~10, and ~15 wt.% of nitrogen. In our FPMD simulations, we use a cubic cell with 96 iron atoms for the pure iron composition, 74 iron atoms + 22 nitrogen atoms (Fe 74 N 22 ) iron alloy with ~7 wt.% nitrogen, 66 iron atoms + 30 nitrogen atoms (Fe 66 N 30 ) for the iron alloy with ~10 wt.% nitrogen, and 56 iron atoms + 40 nitrogen atoms (Fe 56 N 40 ) for the iron alloy with ~15 wt.% nitrogen composition. We used periodic boundary conditions for all the FPMD simulations. We explored pressures up to 400 GPa covering the Earth's core pressure range. We explored several densities, for example, for liquid iron we explored 6.01 to 13.50 g/cm 3 to explore Earth's core pressure (Figure S1 to S4). A more detailed discussion of the methodology can be found in the supporting information (Alfè et al., 2000;Alfè, 2009;Allen & Tildesley, 1987;Balog et al., 2003;Belonoshko et al., 2000;Boehler, 1993;Brown & McQueen, 1986;Ichikawa et al., 2014;Komabayashi, 2014;Laio et al., 2000;Shen et al., 2004;Shimoyama et al., 2013).
We explored pressure volume (density) relation for liquid iron and liquid iron-nitrogen alloy at isotherms ranging in temperature from 4000 to 7000 K. We have used a finite strain Birch-Murnaghan (BM) equation of state (Birch, 1978;Murnaghan, 1944) to describe the pressure-density results for the reference isotherm (T ref = 6000 K) of liquid iron and liquid iron nitrogen alloy (Figure 1, Table S1, and Figure S5). The temperature dependence of the pressure density results could be explained by a Mie-Grüneisen thermal equation of state:
ρ is the temperature derivative of pressure at constant density. For a constant density, we find that pressures (P) and energies (E) increase linearly with increasing temperature. Using the isotherms ranging between 7000 and 4000 K, we estimated the temperature dependence of pressures and energies, that is, dP dT À Á ρ , and dE dT À Á ρ respectively. The heat capacity at constant density (C ρ ) is expressed as dE dT À Á ρ . We also estimate the Grüneisen parameter (γ) and the coefficient of thermal expansivity (α) for liquid iron and nitrogen-iron alloys using the relation:
the thermal expansion is also related to the dP dT À Á ρ by the thermodynamic relation
where K T is the isothermal bulk modulus and is expressed as
. The thermodynamic parameters,
, C ρ and γ, varies as a function of density, pressure, and also chemistry (Table S2). Our results indicate that the thermal Grüneisen parameter (γ) for metallic melts, that is, pure iron liquid, is of the order of ~1.0 and increases upon further compression up to pressures of ~45 GPa. At higher pressures, γ uniformly decreases with increasing pressure (Figure 1). Nitrogen bearing iron liquids also show similar behavior. The reduction of γ in liquid iron alloys as in Fe-N is consistent with recent results in Fe-H liquid (Umemoto & Hirose, 2015). The behavior of γ in pure liquid iron can be directly correlated with atomistic scale changes, such as the pressure dependence of the Fe coordination environment. The average Fe coordination number (CN Fe À Fe ) of pure iron is ~10 at 0 GPa and increases to ~13 at 45 GPa. At pressures >45 GPa, the CN Fe À Fe remains nearly constant (Figure 1). In silicate minerals the γ decreases upon compression; however, 10.1029/2018GL080555
Geophysical Research Letters at pressures where SiO 4 transitions to SiO 6 coordination, the γ increases across the discontinuity. In insulating silicate liquids, upon compression, there is a continuous change in the fraction of SiO 4 and SiO 6 coordination number and as a result the γ increases (Mookherjee et al., 2008;Stixrude & Karki, 2005). In contrast with the insulating silicate liquids, the metallic liquids are nearly close-packed at the outer core condition. Without continuous changes in the CN Fe À Fe , the γ decreases with increasing pressure, similar to crystalline solid iron and mantle silicate minerals (Alfè et al., 2001;Stixrude & Lithgow-Bertelloni, 2005).
At 6000 K, zero pressure bulk modulus (K 0 ) and its first derivative (K 0 0 Þ for liquid iron are 54.0 ± 1.0 GPa and 4.91 ± 0.05, respectively. We find that the bulk modulus, K 0 , decreases, whereas the K 0 0 increases with increasing nitrogen content in the iron liquid (Table S1 and Figure 1). This is similar to the liquid Fe-S where K 0 Figure 1. Structure and property of liquid iron alloy at core pressures. (a) Pressure dependence of the density of pure liquid iron is compared with the liquid nitrogen alloys. Generalized gradient approximation (GGA) results at 6000 K, for pure liquid iron, are shown in red lines. The GGA results for liquid iron with 7, 10, and 15 wt.% N are shown in blue, blue-dash-dotted, and blue dotted lines, respectively. The local density approximation (LDA) results at 6000 K for pure liquid iron (thin red dashed) and liquid iron with 10 wt.% N (thin blue dashed) are also shown for comparison. Symbols represent previous theoretical results for liquid Fe: grey filled triangles (Alfè et al., 2000), grey filled circles (Ichikawa et al., 2014), grey filled squares (Umemoto et al., 2014), and grey filled hexagon (Belonoshko et al., 2000). Shock wave experimental results for liquid iron at ~6000 K are also shown for comparison: grey open diamonds (Brown & McQueen, 1986) and grey open triangles (Anderson & Ahrens, 1994). Inset shows zero pressure bulk modulus, K 0 (open squares) and its pressure derivative, K 0 0 (open circles) for liquid Fe and Fe-N alloys as a function of nitrogen content at 6000 K. The bulk modulus, K 0 ; pressure derivative, K 0 0 ; and uncertainties are reported in Table S1. (b) Thermal Grüneisen parameter, γ, of liquid Fe (red line), Fe-7 wt.% N (yellow line), Fe-10 wt.% N (blue line), and Fe-15 wt.% N (green line). Uncertainties in γ are ~10-15% for pure liquid iron and ~5-10% for nitrogen bearing liquids (not shown in figure for the sake of clarity). The open circles represent Grüneisen parameter for pure iron from shock wave experiments (BM86, Brown &McQueen, 1986, andAA94, Anderson &Ahrens, 1994). The solid black line (V03) is the previous FPMD results for liquid Fe from Vočadlo et al. (2003). The γ for solid iron (A00, Alfè et al., 2001, andV09, Vočadlo et al., 2009). (c) Fe-Fe coordination number for liquid Fe (filled red circles) and liquid Fe-10 wt.% N (filled blue circles). Fe-N and N-Fe coordination numbers are for liquid Fe-10 wt.% N. Fe-Fe coordination for pure liquid iron (S00) is from Sanloup et al. (2000).
Geophysical Research Letters reduces with the increase in sulfur content (Sanloup et al., 2000). The pressure-density behavior of liquid iron from this study is in good agreement, that is, differs less than 1.5% with previous shock wave results (Anderson & Ahrens, 1994;Brown & McQueen, 1986; Figure 1).
We estimate the P wave velocities (V P ) for the liquid iron and the iron nitrogen alloys using the relation
, where K S is the adiabatic bulk modulus and is related to the isothermal bulk modulus, K T by the relation,
. We find that the V P increases as a function of pressure and nitrogen concentration in the iron liquid (Figures 2 andS6). The boundary between the outer and inner core referred to as inner core boundary (ICB) is hotter than the boundary between the silicate mantle and the outer core referred to as core mantle boundary (CMB) by ~2000 K. Thus, the estimates of physical properties such as density along an isotherm are likely to cause a density change Δρ of ~0.36 g/cm 3 (~4%) at CMB and ~0.26 g/cm 3 (~2%) at ICB (Figure S7). Instead, we estimate physical properties along the adiabatic geotherm dT dP À Á s ¼ γT KS (Figure S8). We use pressure and temperature conditions for the ICB to determine the adiabatic geotherm. The melting temperature of pure liquid iron at ICB pressures of ~329 GPa is 6300 ± 500 K (e.g., Alfè, 2009;Anzellini et al., 2013). However, it is likely that temperatures at ICB are lower than the melting temperature of pure iron due to the presence of light elements and associated depression in melting temperatures (Fischer, 2016;Morard et al., 2017). Since the exact amount and nature of light element and depression of melting temperature remain uncertain, we use a range of ICB temperatures, that is, T ICB = 6000 and 5400 K (Hirose et al., 2013), to estimate the adiabatic geotherm (Figure S8). We use the thermal Grüneisen parameter (γ) and adiabatic bulk modulus (K S ) of iron nitrogen alloy to estimate the adiabatic temperature profiles (Figure S8). Using these parameters and with T ICB = 6000 K, the estimated temperature at CMB, that is, T CMB ~4400 K, is in good agreement with previous reports of CMB temperature (Alfè et al., 2007;Fischer, 2016;Hirose et al., 2013).
We estimate the density (ρ) and P wave velocity (V P ) of liquid Fe-N using FPMD and compare with the preliminary reference Earth model (PREM; Dziewonski & Anderson, 1981) to predict the upper limit of nitrogen content of the Earth's outer core. The ρ and V P of nitrogen linearly increases and decreases respectively with increasing nitrogen concentration (Figure 2), that is, dρ dX < 0 and dVP dX À Á > 0, where X = N. Our results indicate that the liquid Fe-N alloy cannot simultaneously reproduce the density and the bulk sound velocity of outer core. Based on our results on Fe-N liquid alloy, the density of liquid iron, and PREM or seismological estimate
, where X = N (mol %), we determine that ~22.5 ± 0.8 mol %, that is, 6.8 ± 0.2 wt.%, nitrogen can satisfy the PREM density at core-mantle boundary for the temperature profile with T ICB = 6000 K (Figure 2). However, based on our results on Fe-N liquid alloy, the V P of liquid iron, and PREM or seismological estimate of the V P at CMB, that is,
, where X = N, we determine that only ~7.4 ± 1.4 mol %, that is, ~2.0 ± 0.4 wt.%, nitrogen can satisfy the PREM. At the ICB and outer core boundary conditions, Fe-N liquid alloy with ~5.7 ± 0.1 wt.% and ~2.6 ± 0.6 wt.% nitrogen can independently reproduce the density and V P , respectively (Figure 2). The concentration of nitrogen derived by matching the V P estimated from PREM and from the iron-nitrogen alloy is always less than that obtained by matching the density. This sets a strict upper limit of permissible nitrogen content of <2.0 ± 0.4 wt.% at CMB and <2.6 ± 0.6 wt.% at ICB conditions. If the nitrogen content was to exceed this, although the density will be better matched, the V P will be overestimated. Given other light elements need to be incorporated to Geophysical Research Letters match V P and density simultaneously, the allowable N content of the outer core is actually much less than 2.0-2.6 wt.%.
Estimated weight percentage of nitrogen to match the seismic density profile would decrease or increase based on the higher or lower temperature conditions of the outer core. If the outer core temperature is lower than 6000 K, then, based on our results, we anticipate even greater concentration of nitrogen to satisfy PREM density profile. For instance, along the temperature profile with T ICB = 5400 K, the estimated temperature at CMB, T CMB ~4000 K. The estimated nitrogen concentration at CMB and ICB are ~7.4 ± 0.2 and ~6.2 ± 0.1 wt.%, respectively. These concentrations are distinctly higher than the estimated N content of the inner core of 2.0-4.8 wt.% that are required to match the density of the inner core for plausible range of inner core temperatures (Adler & Williams, 2005;Litasov et al., 2017). Thus, if density was the only constraint, our finding on nitrogen would be in agreement with the expectation that Earth's outer core contains more light elements than the Earth's inner core (Anderson & Isaak, 2002;Badro et al., 2014). The P wave velocity (V P ) of pure liquid iron and liquid iron nitrogen alloy are insensitive to temperatures at CMB and ICB pressures (Figure S6). Thus, we find that the nitrogen concentration required to match V P remains similar irrespective of the possible range of 6000-5400 K, that is, outer core temperature profile (Fischer, 2016;Hirose et al., 2013).
We have also estimated the amount of other light elements including hydrogen, carbon, oxygen, silicon, and sulfur required to explain the ρ and V P observations of the Earth's outer core (Badro et al., 2014;Umemoto et al., 2014;Umemoto & Hirose, 2015; Figure 3). We determine the light element content in the outer core with the formalism: Δρ ð Þ FeÀCMB ¼ dρ dX ÂX ρ (mol %) where X ρ is the mol % of light element required to match the observed density and
, where X VP is the mol % of light element required to match the observed P wave velocity (V P ). For all the light elements considered, we note that dρ dX < 0 and with an adiabat pegged at T ICB = 6000 K, ~1.1 ± 0.2 wt.% hydrogen, ~5.2 ± 0.2 wt.% carbon, ~5.8 ± 0.03 wt.% oxygen, ~6.9 ± 0.5 wt.% silicon, and ~9.6 ± 0.9 wt.% sulfur could match seismologically observed ρ at CMB (Figure 3). However, only ~0.3 ± 0.2 wt.% hydrogen, ~2.1 ± 1.2 wt.% carbon, ~3.8 ± 2.8 wt.% oxygen, ~2.6 ± 0.8 wt.% silicon, and ~5.8 ± 3.5 wt.% sulfur would be necessary to explain the observed V P at CMB (Figure 3).
Owing to the dependence of density dρ dX and P wave dVP dX À Á velocity with light element composition, it is evident that no single element can simultaneously explain both (ρ) and P wave velocity (V P ) of the Earth's outer core. This is consistent with previous findings (e.g., Badro et al., 2015). It is therefore not surprising that the determined upper limit of N in the outer core in our study is higher than the estimated nitrogen concentration in the Earth's core based on combined cosmochemical and geochemical fractionation-based arguments using nitrogen in the Earth's building blocks (Kerridge, 1985;Pearson et al., 2006). Combining the liquid metal alloy-silicate partitioning coefficients, D alloy=silicate N e2:0 ð Þobtained from FPMD simulations (Zhang & Yin, 2012) and nitrogen abundance of the BSE (McDonough, 2003), the nitrogen content of the core could be estimated to be ~0.0001-0.01 wt.%. However, more recent estimate of alloy-silicate partition coefficient of nitrogen, D alloy=silicate N of ~15-20 (e.g., Dalou et al., 2017;Kadik et al., 2011;Li, Dasgupta, et al., 2016;Roskosz et al., 2013), carbonaceous chondrite nitrogen abundance of ~1,500 ppm as the bulk Earth nitrogen content (Kerridge, 1985;Pearson et al., 2006), and homogeneous accretion of Earth leads to bulk core N content of 0.41-0.42 wt.%. If volatility-induced nitrogen loss is considered from chondritic building blocks, then the bulk N content of the core would be even lower. However, our estimated upper bound is much less than the carrying capacity of nitrogen in high P-T Fe-alloy liquid as suggested by recent solubility experiments (Speelmanns et al., 2018). The tradeoffs and strict bound provided in this study on the nitrogen content of the Earth's core based on the effect of N on physical properties of liquid Fe-alloy have implications for the available budget of nitrogen during core formation, the nature of terrestrial building blocks, and the extent of siderophile behavior N may display at more extreme conditions than what the existing metal-silicate experiments covered thus far.
BAJGAIN ET AL.