Hexagonal close-packed iron (hcp-Fe) is a main component of Earth’s inner core. The difference in density between hcp-Fe and the inner core in the Preliminary Reference Earth Model (PREM) shows a density deficit, which implies an existence of light elements in the core. Sound velocities then provide an important constraint on the amount and kind of light elements in the core. Although seismological observations provide density–sound velocity data of Earth’s core, there are few measurements in controlled laboratory conditions for comparison. We report the compressional sound velocity (V P ) of hcp-Fe up to 163 GPa and 3000 K using inelastic x-ray scattering from a laser-heated sample in a diamond anvil cell. We propose a new high-temperature Birch’s law for hcp-Fe, which gives us the V P of pure hcp-Fe up to core conditions. We find that Earth’s inner core has a 4 to 5% smaller density and a 4 to 10% smaller V P than hcp-Fe. Our results demonstrate that components other than Fe in Earth’s core are required to explain Earth’s core density and velocity deficits compared to hcp-Fe. Assuming that the temperature effects on iron alloys are the same as those on hcp-Fe, we narrow down light elements in the inner core in terms of the velocity deficit. Hydrogen is a good candidate; thus, Earth’s core may be a hidden hydrogen reservoir. Silicon and sulfur are also possible candidates and could show good agreement with PREM if we consider the presence of some melt in the inner core, anelasticity, and/or a premelting effect.

Keywords

The sound velocity of pure hcp-Fe at high temperatures provides an important constraint on the composition of Earth’s core: any model of the impurity composition of the core must both satisfy the density deficit from seismology and give an appropriate velocity. Here, we report the V P of hcp-Fe up to 163 GPa and 3000 K ( Fig. 1 ) based on a combination of the laser-heated diamond anvil cell (LHDAC) ( 17 ) and IXS measurements ( 18 ).

The sound velocity (V P ) profile of Earth’s core is important information obtained from seismological observations, and the density-velocity relation for Earth’s core is well known. To constrain the abundances of light elements in Earth’s core, Birch’s law, a linear relation between sound velocity and density ( 5 ), has been used to extrapolate measured sound velocities to inner core conditions. To do this, the sound velocities of Fe and Fe alloys have been measured at high pressure and moderate temperatures (<2000 K) by inelastic x-ray scattering (IXS) ( 6 – 13 ) and nuclear inelastic scattering (NIS) ( 14 – 16 ). It is then desirable to extend these measurements as close as possible to core conditions. However, this is difficult, and sound velocity measurements with static compression at temperatures above 2000 K have not been reported.

Seismological data reveal important physical properties of Earth’s core, such as density and elasticity. The core is mainly composed of iron (Fe), and the stable crystal structure of Fe under inner core conditions is hexagonal close-packed (hcp) ( 1 ). The density of the inner core in the Preliminary Reference Earth Model (PREM) ( 2 ), which is based on seismological data, is about 2 to 5% smaller than that of hcp-Fe (a core density deficit) ( 3 ), and it is accepted that the inner core consists of iron and light elements, such as silicon, sulfur, oxygen, hydrogen, and carbon ( 4 ).

Figure 2B compares the present results to previous IXS measurements ( 6 , 10 ) and ab initio calculations ( 20 , 21 ). The present isothermal ρ-V P fitting lines are consistent with previous IXS ( 6 , 10 ). In the case of ab initio calculation, the lines have a good agreement with the results ( 20 , 21 ) above 4000 K, although some discrepancies exist in the lower temperature range (<4000 K).

We choose T 0 to be 300 K; M and B are the coefficients of Birch’s law at room temperature, whereas A and ρ* embody the temperature dependence. This form is convenient for the consideration of errors: the large data set at room temperature fixes M and B relatively precisely, whereas the sparse data at high temperatures are reflected in the uncertainty in A and ρ*. The appearance of the crossover density ρ* also explicitly indicates a maximum upper bound for the validity of this relation. Fitting our data, combined with shock compression data ( 19 ) along the Hugoniot, gives M = 1.160 ± 0.025, B = −3.43 ± 0.29, A = 7.2 × 10 −5 ± 3.6 × 10 −5 , and ρ* = 14.2 ± 1.5. Details of the fitting are presented in the Supplementary Materials.

( A ) Results of the present study and of Ohtani et al. ( 12 ) based on IXS and shock compression experiments as reported by Brown and McQueen ( 19 ). Isothermal ρ-V P fitting lines are expressed as V P = [1.160ρ − 3.43] + [7.2 × 10 −5 × (T − 300) × (ρ – 14.2)]. ( B ) Comparison of obtained isothermal ρ-V P fitting lines with previous studies of NIS reported by Mao et al. ( 10 ), diamonds; IXS reported by Antonangeli et al. ( 6 ), inverted triangles; ab initio calculation reported by Vočadlo et al. ( 20 ), open crosses; and report of Sha and Cohen ( 21 ), crosses.

We fit our data by allowing the coefficients in the linear ρ-V P relation to vary with temperature, which preserves Birch’s law at any fixed temperature ( Fig. 2A ). A good fit to the data is possible by taking a first-order, linear, temperature dependence of the coefficients. The preferred orientation of the sample in the DAC is estimated using x-ray diffraction (XRD) pattern, and the velocity variation is less than 1.5% at 163 GPa and 3000 K (see also the discussion in the Supplementary Materials). We then parameterize the dependence as (1)

The values of hcp-Fe and each Fe–light element compound [hcp-Fe 92 Ni 8 ( 14 , 24 ), purple; dhcp-FeH ( 11 , 25 ), blue; hcp-Fe 85 Si 15 ( 10 , 26 ), aqua; Fe 3 S ( 13 , 27 ), green; FeO (B1/rhombohedral phase) ( 7 , 28 ), yellow-green; Fe 3 C ( 9 , 29 ), pink; Fe 7 C 3 ( 30 ), yellow] show distributions from the ICB to the COE. The obtained ρ-V P of hcp-Fe at 5500 K is indicated by orange, and the star is at the expected ICB conditions. The other stars denote ρ-V P of Fe–light element compounds under ICB conditions assuming that the temperature effects of these compounds are the same as those of hcp-Fe (the dashed lines are the temperature effect). Triangles connecting three stars (hcp-Fe–hcp-Fe 92 Ni 8 –Fe–light element compound) indicate the potential ρ-V P region given by mixing three components.

Several sound velocity and density measurements of pure iron and iron alloys demonstrate the effect of light elements on the properties of iron. The ρ-V P plots of hcp-Fe ( 3 , 12 ), hcp-Fe 92 Ni 8 ( 14 , 24 ), dhcp-Fe ( 11 , 25 ), hcp-Fe 85 Si 15 ( 10 , 26 ), Fe 3 S ( 13 , 27 ), FeO (B1/rhombohedral phase) ( 7 , 28 ), Fe 3 C ( 9 , 29 ), and Fe 7 C 3 ( 30 ) are summarized in Fig. 4 . The ρ-V P plots of hcp-Fe at 4200 and 5500 K are also plotted in this diagram, and a star symbol indicates the ICB condition. Each ρ-V P plot of iron alloys and iron–light element compounds is shown from the center of Earth (COE) to the ICB. The stars indicate the extrapolation to the ICB conditions (T ICB = 5500 K) assuming that the temperature effect on these compounds is the same as that on hcp-Fe. This assumption is made on the basis of the few studies that show the temperature effects on iron alloys based on Birch’s law. To account for the composition of Earth’s inner core, the triangles (hcp-Fe–hcp-Fe 92 Ni 8 –Fe–light element compounds) need to overlap with PREM. In this discussion, we assumed Fe alloy to be an ideal mixture among various Fe-rich compounds, such as hcp-Fe, hcp-Fe-Ni, and Fe–light element alloys to estimate the average ρ-V P of Fe alloy with light elements. Oxygen and carbon may not be suitable as major light elements in the core. On the other hand, silicon, sulfur, and hydrogen decrease both density and V P and are thus potential candidates. Hydrogen is a good candidate because the mixing of hcp-Fe, hcp-Fe 92 Ni 8 , and dhcp-FeH can be accountable for the ρ-V P plot of PREM. Thus, the inner core may be a hidden hydrogen reservoir in Earth. Silicon and sulfur can also meet an important requirement if their temperature effects on V P are larger than that of hcp-Fe. In addition, the presence of melt in the inner core ( 31 , 32 ), anelasticity ( 32 ), and/or a premelting effect ( 33 ), which causes the decrease in sound velocities, can support a silicon- or sulfur-rich inner core.

Stars indicate V P and ρ of hcp-Fe and PREM at 330 GPa. The difference in the ρ-V P plot between PREM and hcp-Fe at 5500 K (orange) shows a 4 to 5% core density deficit and a 4 to 10% core velocity deficit and that of 4200 K (blue) shows a 5 to 6% core density deficit and a 7 to 13% core velocity deficit.

To demonstrate the impact of this work for understanding core composition, we focus on ρ and V P at the inner core boundary (ICB), where the pressure has been determined to be 330 GPa by a seismological study ( 2 ). The recent discussion of the thermal state of the outer core based on the melting experiments of iron–light element systems suggests that the temperature at ICB is around 5500 K ( 22 ). The density of hcp-Fe can be calculated using the equation of state ( 3 ). The density-velocity relation from the present work extrapolated to this temperature is shown in Fig. 3 . Also shown are the results of the PREM ( 2 ). The phenomenological Eq. 1 may be used for extrapolation to the inner core conditions (see also the discussion in the Supplementary Materials). As is apparent from Fig. 3 , both the density and V P of hcp-Fe are higher than PREM values: Earth’s inner core has a 4 to 5% smaller density and a 4 to 10% smaller V P than hcp-Fe at 5500 K. Thus, on the basis of the present sound velocity results, we conclude that light elements, or a combination of light elements and nickel, in the inner core decreases not only density but also simultaneously the V P of hcp-Fe under inner core conditions. Assuming the lower core temperature model (T ICB = 4200 K) proposed on the basis of the melting temperature of pyrolite ( 23 ), Earth’s inner core has a 5 to 6% smaller density and a 7 to 13% smaller V P than hcp-Fe at 4200 K.

MATERIALS AND METHODS

Samples and cells The IXS experiments and XRD experiments were conducted up to 163 GPa and 3000 K at the same beamline, BL35XU. High-pressure and high-temperature conditions were generated using a double-sided LHDAC with 150-and 250-μm culets. The used laser heating system is COMPAT (17), which has been developed for an LHDAC in a limited experimental space, such as the beamline of a synchrotron radiation facility. Rhenium gaskets were preindented to 35 to 70 μm in thickness, and a hole 50 to 100 μm in diameter was drilled in the gasket as the sample chamber. The iron sample (99.99% purity) was compacted from a powder into a foil 45 to 90 μm in diameter and 10 to 30 μm in thickness. The sample was sandwiched between two NaCl layers that are 10 μm in thickness, which were used as a pressure medium and a thermal insulator.

Laser heating The average temperature of the heating spot obtained at 118 GPa and 3000 K is shown in fig. S1. The typical temperature uncertainty is 200 K due to fluctuation. The experimental duration at one pressure-temperature condition was about 4 to 6 hours. The temperature was monitored and controlled for 4 to 6 hours and recorded every 15 min during the heating experiments. Sample temperature under laser heating was determined by fitting Planck’s formula to a spectrum of thermal radiation from a sample between 600 and 800 nm. The optics used to collect thermal radiation was calibrated with an intensity standard lamp. The hot spot is greater than the x-ray beam size at full width at half maximum (FWHM). We calculated the sample temperature by averaging the variation in the heating area probed by x-rays during the IXS measurements that lasted 4 to 6 hours.

IXS and XRD The IXS and XRD experiments were performed at BL35XU of SPring-8 in Japan (18). We used the Si (9 9 9) backscattering optics, which provides an incident photon energy of 17.794 keV with an energy resolution of 2.8 meV at FWHM. The beam size was focused to 16 × 16 μm by a Kirkpatrick-Baez mirror pair for this work (34). The scattered x-rays were analyzed by 12 crystals, which are arranged in a two-dimensional (3 × 4) array. The momentum transfer Q = 2k 0 sin(2θ/2), where k 0 is the wave vector of the incident photons and 2θ is the scattering angle, was selected by rotating the spectrometer arm in the horizontal plane. Here, we collected IXS spectra in the range of Q = 6.0 to 10.0 nm−1 at each experimental pressure. The momentum resolution was set to about 0.4 nm−1 full width. IXS spectra of hcp-Fe at 163 GPa and 3000 K are shown in fig. S2. The XRD patterns were obtained under the same experimental conditions as IXS using a flat panel area detector installed at BL35XU. An example of XRD pattern is shown in fig. S3. Rediffraction peaks may come from tails of x-ray beams. Pressure was calculated on the basis of XRD patterns using the equation of state of hcp-Fe (3). The lattice preferred orientation of this compressed sample was estimated using this diffraction pattern and discussed below. The experimental conditions are summarized in table S1.

Fits to the dispersion Figure 1 shows the IXS spectra collected at 163 GPa and 3000 K. The spectra are characterized by an elastic contribution centered at zero energy and inelastic contributions from hcp-iron and diamond. The peak of the TA phonons of diamond appears at higher energies compared to hcp-iron because of the much higher sound velocity of diamond. The zero energy was well determined by the elastic peak. The energy positions of photons were extracted by fitting the data with a set of Lorentzian functions. The compressional sound velocity (V P ) was determined by fitting the phonon dispersion with a sine function (S1)where E and Q are the energy and the momentum of the acoustic mode, and Q MAX is the first Brillouin zone edge (35). The expression can be derived from the solution of the dynamical matrix limited to the nearest-neighbor interaction within the framework of the Born–von Karman lattice dynamics theory (36). Here, Q MAX was taken to be a free parameter in the fit, as in previous studies (9). V P can change when Q MAX is fixed to be Qa/2π = 0.6, but it remains within the error range for the fits to V P when Q MAX is free (fig. S4). Because our sample was polycrystalline, the determined velocity is an aggregate sound velocity, averaging over crystal orientations. We made critical assessments of texture effects and deviation from ideal average velocity as shown below. Figure S5 shows the dispersion curve at each experimental condition and the fits using eq. S1. The phonon energy of inelastic scattering increases with pressure, which indicates an increase of the compressional sound velocity with pressure. The results are summarized in table S2. Our 12 independent data points (albeit at mainly four different momentum transfers) provide a good constraint on the sound velocity. In addition, we conducted extra experiments to check the Q dependence of the obtained sound velocity. The difference in velocity at higher (IXS101) and lower (IXS102) Q in the same conditions is ~0.35%, as shown in fig. S6.

Birch’s law Birch’s law is a linear relation of the compressional sound velocity V P with density ρ as V P = aρ + b (4, 5), which arises in the physics of solids as a linearization of a power law (37). No regard is given to the temperature effect on the equation. To cover ρ-V P data under a wide range of temperature conditions, the temperature effect on Birch’s law should be considered. A recent ab initio calculation study (33) reports the premelting effect on sound velocities of hcp-Fe. Their results indicate that the temperature dependence of the velocity is linear below melting temperature (T/T m < 0.9). Assuming the linear temperature effect on the slope, a, and intercept, b, we modify Birch’s law as (S2) Equation S2 can be rewritten in the following form (S3) Here, we replaced a 0 with M, b 0 with B, and a 1 with A, and we refined –b 1 /a 1 as ρ* (S4) Equation S4 can be divided into a temperature-independent term, [Mρ + B 0 ], and a temperature-dependent term, [A(T − 300)(ρ − ρ*)]. The first term shows typical Birch’s law at 300 K. The second one exhibits the temperature effect on V P at each density (pressure). At the same density, V P decreases with temperature. In addition, the decrease in V P at the same temperature becomes smaller with an increase in density. However, the linear temperature dependence implies that these relations should only be applied below a density, given by ρ* in eq. S4: at ρ*, V P becomes independent of temperature, and above ρ*, the temperature dependence becomes unphysical, with velocity increasing as temperature increases. ρ* corresponds to a pressure of about 500 GPa, far beyond the center pressure of Earth. To approximate our results and extrapolate them to different conditions, we used the linear form discussed above and simultaneously fit the IXS results at high temperatures reported here and at room temperature (12) and the shockwave data (19). The error bars of the shockwave data are not given precisely; hence, we consider a velocity error of 1 or 2%. Table S3 gives the resulting parameters. To check the validation of fitting parameters, we compared fit results with experimental data. Figure S7 shows that the fit results are within error bars of each experimental data; thus, the fitting parameters are reasonable. Notably, as is often the case for linear fits, the normalized covariance matrix of the fit parameters tends to be highly correlated, with diagonal values that approach unity magnitude. To estimate the error in the value at the core-mantle boundary (CMB), we therefore used the error propagation equation (38) with the full covariance matrix of the fit parameters, which, in this case, is (S5)where σ AB is the covariance of A and B and the error in the velocity V(ρ,T) is given by . The results given in table S3 suggest that the velocity determination is fairly robust to small changes in the errors of the shockwave data. Fitting using our high-temperature data (2300 and 3000 K) and shock compression data provided the parameters of eq. S4 in table S3. However, in the case of the combination of ambient temperature data (300 K) and shock compression data, fitting was not completed successfully. The fitting curves in Fig. 2 were obtained from the data at 2300, 3000 (this study), and 300 K (12) and from shock data (19). We checked the fitting by a power law as suggested by Mao et al. (10). Our analyses revealed that the fitting by the power law shows almost no curvature and is very close to the linear extrapolation, suggesting that our linear extrapolation is valid. Equation 1 gives the isothermal ρ-V P relation by plugging the temperature. To compare with PREM, we assumed that the CMB temperature is 5500 K (22). The density at 330 GPa and 5500 K was calculated to be 13.42 ± 0.08 g/cm3 using the equation of state of hcp-Fe (2). Thus, the velocity of hcp-Fe at 330 GPa and 5500 K is 11.85 ± 0.40 km/s. In the case of 5000 K, ρ 330GPa, 5000K = 13.49 ± 0.08 g/cm3 and V P, 330GPa, 5000K = 11.98 ± 0.40 km/s, whereas at 6000 K, ρ 330GPa, 6000K = 13.34 ± 0.08 g/cm3, that is, V P, 330GPa, 6000K =11.69 ± 0.40 km/s. The ρ and V P of PREM are 12.764 cm3 and 11.028 km/s, respectively (2). In the temperature region between 5000 and 6000 K, the hcp-Fe indicates larger density and compressional sound velocity than PREM at 330 GPa. Therefore, the density and velocity deficits in Earth’s inner core compared to hcp-Fe are [Δρ (hcp-Fe – PREM) /ρ hcp-Fe ] 5500K = 4 to 5% and [ΔV P(hcp-Fe – PREM) /V Phcp-Fe ] 5500K = 4 to 10%, respectively. In an analysis for the conclusion, we assumed that the temperature effect on Birch’s law for all iron–light element compounds is the same as that for hcp-Fe. The present conclusion on the light elements in the inner core does not change without unusual reduction of the sound velocity of FeO and Fe 3 C with increasing temperature, which is not likely (39–41).

Ab initio calculation To evaluate the effect of the preferred orientation on the sound velocity measured in this work, we calculated the anisotropy of the elastic constants of hcp-Fe under the conditions of the present work. We calculated the elastic constants of hcp-Fe at 163 GPa and 3000 K based on an ab initio calculation: 1321.5 GPa for C 11 , 661.3 GPa for C 12 , 620.4 GPa for C 13 , 1419.7 GPa for C 33 , and 211.7 GPa for C 44 . The C ij calculated here at 163 GPa and 3000 K is consistent with previous work under similar conditions. Details of the ab initio calculation are given by Tsuchiya and Fujibuchi (42). The difference between C 11 and C 33 of hcp-Fe becomes smaller with temperature because C 11 is 1552.8 GPa and C 33 is 1767.0 GPa at 180 GPa and 0 K. The difference will become much smaller at higher temperatures and hcp-Fe will become isotropic, as reported by Sha and Cohen (21).