The influence of δ-(Al,Fe)OOH on seismic heterogeneities in Earth’s lower mantle

The high-pressure phases of oxyhydroxides (δ-AlOOH, ε-FeOOH, and their solid solution), candidate components of subducted slabs, have wide stability fields, thus potentially influencing volatile circulation and dynamics in the Earth’s lower mantle. Here, we report the elastic wave velocities of δ-(Al,Fe)OOH (Fe/(Al + Fe) = 0.13, δ-Fe13) to 79 GPa, determined by nuclear resonant inelastic X-ray scattering. At pressures below 20 GPa, a softening of the phonon spectra is observed. With increasing pressure up to the Fe3+ spin crossover (~ 45 GPa), the Debye sound velocity (vD) increases. At higher pressures, the low spin δ-Fe13 is characterized by a pressure-invariant vD. Using the equation of state for the same sample, the shear-, compressional-, and bulk-velocities (vS, vP, and vΦ) are calculated and extrapolated to deep mantle conditions. The obtained velocity data show that δ-(Al,Fe)OOH may cause low-vΦ and low-vP anomalies in the shallow lower mantle. At deeper depths, we find that this hydrous phase reproduces the anti-correlation between vS and vΦ reported for the large low seismic velocity provinces, thus serving as a potential seismic signature of hydrous circulation in the lower mantle.

www.nature.com/scientificreports/ depths, causes a softening of the isothermal bulk modulus, which consequently decreases v Φ

22
. Moreover, a theoretical study proposed a negative correlation of pressure (P)-shear modulus (μ) for low-spin ε-FeOOH, which causes a decrease of v S /v Φ with increasing pressure 26 . These studies point to the possibility that the solid solution may in part be responsible for regionally seismic heterogeneities observed in the lower mantle 27,28 .
As discussed above, δ-(Al,Fe)OOH is likely stable in subducted slab materials throughout the lower mantle 18 , yet the possible effects on the seismic wave velocities of a phase assemblage containing the δ-phase are not well constrained. In this study, we report the elastic wave velocities of the δ-(Al , Fe)OOH phases having an Fe/(Al + Fe) ratio of ~ 0.13 determined by nuclear resonant inelastic X-ray scattering (NRIXS). The Fe/(Al + Fe) ratio of our sample is within the rations of δ-phase formed in the MgO-Al 2 O 3 -Fe 2 O 3 -SiO 2 -H 2 O system that simplifies natural basaltic compositions (Fe/(Al + Fe) = 0.10-0.17) 18 . By combining our results from NRIXS with the equation of state 22 , we discuss the potential relationships between δ-(Al,Fe)OOH and seismic anomalies in the lower mantle.

Results
Determination of elastic wave velocities. NRIXS spectra of the two δ-(Al 0.873 57 Fe 0.127 )OOH (δ-Fe13-r1) and δ-(Al 0.867 57 Fe 0.133 )OOH (δ-Fe13-r2) samples were obtained up to 79 GPa (Fig. S1). The 57 Fe-partial projected phonon density of state (pDOS) was subsequently extracted from each NRIXS spectrum using the PHOENIX software 29 (Fig. 1). Generally, vibrational peak positions in the pDOS for minerals steadily shift to higher energies with compression 30,31 . However, the peak positions of pDOS of δ-Fe13 shift to lower energies with increasing pressure from ambient pressure to 10.8 GPa, above which they shift to higher energies (Fig. 1). This pressure at which the peak position of pDOS of δ-Fe13 reaches the lowest energy is very close to that of the structural transition from P2 1 nm with asymmetric hydrogen bonds to Pnnm with symmetric hydrogen bonds inferred from XRD experiments on δ-(Al,Fe)OOH (Fe/(Al + Fe) = 0.047(10) (hereafter, δ-Fe5) and 0.123(2) (δ-Fe12)) 22 . Static calculations for pure δ-AlOOH predicted that phonon-softening occurred as a result of this transition involving hydrogen bonding 32 . Through this transition, the lowest and highest A 1 and B 2 optic modes soften and the OH stretching frequency approaches nil 32 . The Raman spectra of polycrystalline δ-AlOOH show broadening or disappearance of B 2 mode of P2 1 nm and the appearance of A g mode of Pnnm near 5.6 GPa 20 . Therefore, the shifting of the pDOS to lower energies below about 10 GPa and the trend observed in the Lamb-Mössbauer factor (f LM ) of δ-Fe13 (Fig. 2) are likely related to this transition. Due to the likelihood that one or more optical branches cross into the accessible low-energy regime in this low-pressure range (~ 2 < E (meV) < ~ 15), the assumption that the low-energy region of the pDOS is occupied purely by acoustic modes breaks down. Thus, we restrict our determination of the Debye sound velocity (v D ) to ambient pressure and pressures higher than 20 GPa.
A representative Debye sound velocity determination is shown in Fig. 2a. In this figure, two different models for phonon dispersion are plotted; one model is a 'Debye-like' model that plots as a constant value, and the other is an empirical power law model for phonon dispersion defined as Eq. 1: where A 1 and A 2 are constants that are optimized by a standard least squares fitting process. In this power law model, v D is defined as v(0) in this profile. Note that 'Debye-like' dispersion is accommodated in Eq. (1).
However, low-count rate experiments such as these pose challenges to determining v D . For example, the values determined from these two dispersion models may be sensitive to the fitting energy range for pDOS and the phonon dispersion model, especially if there is scatter in the data. To place better constraints on v D without relying on the choice of an energy fit range, the analysis of a probability distribution function (PDF) of v D proposed by Morrison et al. (ref. 33 ) was applied in this study and an example at 67 GPa is shown in Fig. 2. The PDF method accounts for hundreds of physically reasonable fit ranges and phonon dispersion behavior in the determination of v D (Eq. (1)), and therefore places a more reasonable constraint on v D and its uncertainty than would be calculated only from a single energy range 33 . The details of the PDFs, including the energy ranges considered and bin size are provided in Table S1. The reported v D values and their uncertainty are the peak positions and full width at half maximum of fitting an asymmetric function to the PDFs, respectively (Figs. 2b, S2; Tables S1 and S2). The 57 Fe concentration of the sample was used in the determination of v D .

Behavior of elastic wave velocities under compression.
The v P , v S , and v Φ values of the two δ-Fe13 samples are summarized in Fig. 3 and Table S2. We make the following observations and compare with previous results. The concentration of iron in our sample is three times higher than the δ-(Al 0.956 Fe 0.044 )OOH polycrystalline sample (hereafter, δ-Fe4) measured using Brillouin light spectroscopy, thus when compared with the Al end-member values and the δ-Fe4 sample, our ambient pressure v P and v S values are consistent with this trend. The v S of δ-Fe13 increases up to the pressure conditions of the Fe 3+ spin transition 22 , although the slope is shallower than that reported for δ-Fe4 23 and δ-AlOOH 20 (Fig. 3). The lower v S values we report are reasonably explained by the difference in Fe concentration, as noted above. However, other factors such as preferred orientation could offer explanations. The NRIXS measurements used crushed single grains of a more iron-rich composition and the Brillouin scattering measurements used a potentially finer-grained polycrystalline iron-poor sample 23 . These differences might enhance the effect of orientation and/or intergranular stresses on measured data in the vicinity of phase transitions. Interestingly and unlike other materials undergoing a spin crossover, the polycrystalline Brillouin scattering results do not show any softening of v P near this transition 23 . At 32 and 37.1 GPa, both v P and v Φ decrease due to the reduction of K influenced by Fe 3+ spin transition 22 (Fig. 3a,c and Table S2). A decrease in K as a result of the spin transition is also observed in other low-Fe content minerals (e.g., ferropericlase 34,35 ). The degree of reduction of v P and v Φ of δ-Fe13 due to the spin transition is significant. Specifically, the v P and v Φ values at 37.1 GPa are decreased by 13% and 23% compared to those at 27 GPa. After  2). We note that the fit models plotted here are purely demonstrative and do not directly reflect our reported Debye velocities. Families of these models are considered, fit to the data, and plotted as probability distributions (ref. 33 The energy ranges applied in the PDF analysis at 67 GPa fall within, but are not limited to, the E min and E max values given in Table S1. See text for a detailed description of our Debye velocity determination method. www.nature.com/scientificreports/ the completion of the Fe 3+ spin transition (45 GPa and higher pressures), both v P and v Φ increase steadily with compression ( Fig. 3a,c). Above 45 GPa, the pressure-velocity profiles of δ-Fe13 do not show any abrupt changes (Fig. 3). Pressure dependencies of velocities, however, are significantly different among v P , v S , and v Φ . The v P and v Φ exhibit relatively large gradients with pressure, similar to the Preliminary Reference Earth Model (PREM) 36 . In contrast, the v S values are essentially invariant with respect to pressure, such that the values of v S determined in this study above 45 GPa are equal within the associated errors (6.38-6.55 ± 0.10-0.18 km/s). This leads to a weak pressure dependence of the shear modulus (μ). The μ increases from 190 to 206 GPa between 45 and 79 GPa, which is within the range of errors: ± 0.06 to ± 0.11 (Table S2).

Discussion
The spin transition of Fe 3+ in the δ-Fe13 samples results in a high v S /v Φ (0.86 at 37.1 GPa) (Table S2), which is ~ 19% higher than PREM at this corresponding depth (~ 970 km). The spin transition of Fe 3+ in the octahedral site is also observed in the Fe-bearing NAL phase 37 . However, pressure condition and width of spin crossover involving volume collapse are slightly lower and narrower in the δ-Fe5 and δ-Fe12 (32-40 GPa) 22 37 . A computational study by ref. 38 shows that the spin transition pressure of Fe 3+ in the octahedral site of the NAL phase is ~ 40 GPa, and it remains mostly invariant to temperature and the width moderately increases with temperature. At 300 K, the isothermal bulk modulus of the Fe 3+ -bearing NAL phase (NaMg 2 (Al 4.67 SiFe 0.33 )O 12 ) in the pressure range of spin crossover is up to ~ 25% lower than the bulk modulus-pressure trend expected without a spin transition 38 . The reduction of bulk modulus is limited to ~ 8% at 1200 K 38 , the estimated temperature of a slab penetrated into the top of the lower mantle 39,40 , resulting in the ~ 4% reduction of v Φ . At 300 K, the isothermal bulk modulus of δ-Fe13 throughout the spin crossover is up to ~ 55% lower than the bulk modulus-pressure trend expected without a spin transition, resulting in the ~ 33% reduction of v Φ . This would imply that for δ-Fe13 at about 1200 K the pressure of the spin-crossover is likely to be unchanged from that measured at 300 K and the resultant softening of the bulk modulus remains appreciable compared to the Fe 3+ -bearing NAL phase 37 . This behavior is different than the Fe 2+ spin crossover in ferropericlase, in that although ferropericlase exhibits elastic softening, the transition occurs at deeper depths 34,35 and exhibits a stronger temperature dependence 41 than that for the Fe 3+ -bearing phases discussed above. It is also different from that of Fe 3+ -bearing bridgmanite 42 and stishovite 43 .
Seismic observations in this region just below the transition zone (< 1000 km in depth) suggest a range of anomalies 44 , inferred to be correlated with the presence or stagnation of slab debris 45 . Therefore, the high v S /v Φ in the spin crossover region of δ-Fe13 (and NAL) could be related to these seismic anomalies, especially in regions of relatively lower temperatures, such as those calculated for subducted slabs. Further studies on the temperature dependence of these properties across the spin transition would help quantify this effect.
Above 45 GPa, the gradients of v P and v Φ (dv P /dP and dv Φ /dP) of δ-Fe13 are close to those of PREM (Fig. 3a,c). By contrast, that of v S (dv S /dP) is practically constant above 45 GPa (Fig. 3b). To discuss velocity relations between δ-Fe13 and PREM at deep-lower mantle pressure conditions, we calculated the elastic wave velocities for pressures higher than 80 GPa by using the following extrapolation procedure. First, we determined the elastic wave velocities at room temperature (Table S2). The v P values were extrapolated using a linear relationship (Birch's law 46 ) between the ρ calculated from the thermal EoS (see 'Determination of elastic wave velocities' section) and the measured v P of low-spin δ-Fe13 at 45 GPa, 59 GPa, 67 GPa, and 79 GPa (Table S2). The v Φ was calculated from those thermal EoS. The v S was then calculated from the extrapolated v P and the calculated v Φ by using the Eqs. (5)- (7).
The v S of δ-Fe13 at room temperature is 5% slower than PREM at 1871 km in depth, and then becomes 9% slower than PREM at 2771 km in depth, 120 km above the core-mantle boundary. The difference of v S between  www.nature.com/scientificreports/ δ-Fe13 and PREM becomes larger with depth because of the very week pressure dependence of shear modulus of δ-Fe13 (Fig. 4). On the other hand, v P and v Φ at room temperature are 1% and 6% faster than PREM at 1871-2771 km in depth, respectively (Fig. 4, Table S3). The difference of v S and v Φ between δ-Fe13 and PREM suggests that an incorporation of 17 vol.% δ-Fe13 into PREM (representing average mantle) can account for − 0.9 to − 1.5% anomaly for v S and + 1% anomaly for v Φ in this depth range if the temperature dependences of velocities of δ-Fe13 are not considered. Such an anti-correlation between the v S and v Φ has been reported for the large low-shear-velocity provinces (LLSVPs) in the depth range of 2000 km down to the core-mantle boundary, beneath the Pacific Ocean and the region spanning underneath the Atlantic Ocean to the western and southern part of the African continent 27,28 . It has been suggested that the LLSVPs are surrounded by downwellings, driven by slab subduction 47 . Therefore, if slabs contain a hydrous δ-phase and these slabs are transported into the deep-lower mantle, then slab debris containing the δ-phase could accumulate in the LLSVPs over geologic time, although it is not clear what the exact thermal state of this assemblage would be. To compare the elastic wave velocities of δ-Fe13 and PREM along an average mantle geotherm 48 , we apply the temperature dependence of v S , which was reported for MgSiO 3 -bridgmanite at 100 GPa (i.e., dv S /dT = − 1 × 10 -4 (km/s)/K) at 100 GPa 49 . Another study demonstrated that the dv S /dT of newly discovered FeO 2 H x P-phase is also approximately − 1 × 10 -4 (km/s)/K, at 133 GPa 50 . The v Φ along the mantle geotherm was calculated from the thermal EoS parameters used in the calculation of v Φ at room temperature, and v P was subsequently calculated from the extrapolated v S and v Φ along an average mantle geotherm by using the Eqs. (5)- (7). In this case, the v P and v S of δ-Fe13 along the mantle geotherm 48 are 1% and 8 to 11% slower than PREM in the depth range of 1871-2771 km, respectively, whereas v Φ of δ-Fe13 is 4 to 5% faster than PREM (Fig. 4, Table S3). The anti-correlation between v S and v Φ results in ~ 30% lower (v S − v Φ )/v Φ ratio for δ-Fe13, compared with the PREM model (Fig. 4).
In this discussion, we used γ 0 = 0.64 (with q = 1.8 and Θ 0 = 1485 K reported for δ-AlOOH 11 ). However, the effect of spin transition on the Grüneisen parameter of δ-Fe13 has not been investigated. To address the uncertainty of the Grüneisen parameter in our extrapolations, we tested a value of 1.28 as γ 0 which is twice as large as the γ 0 used in the present discussion (0.64) . The doubled value of the Grüneisen parameter decreases the v Φ , the most sensitive parameter of all three velocities (v P , v S , and v Φ ) to changes in the Grüneisen parameter, by ≤ 0.12 km/s (1.2% reduction) at the pressure and temperature conditions of 1871-2771 km depth, which are smaller than the errors of the extrapolated velocities ( Fig. 4 and Table S3). Such a significant change in the Grüneisen parameter is unlikely for δ-Fe13. For example, a recent calculation for (Mg 0.75 Fe 0.25 )O ferropericlase showed that at 120 GPa the Grüneisen parameters of high-and low-spin states were 1.94 and 2.05 at 1400 K, respectively 51 . As discussed above, δ-Fe13 is primarily in the low-spin state at 80 GPa under elevated temperatures. If the difference in the Grüneisen parameter between high-spin and low-spin δ-Fe13 is also ~ 0.1, the velocity reductions are smaller than the example we tested above. To place a better constraint on the Grüneisen parameter of Fe-bearing δ-phase through the spin transition, high-temperature experiments and computational studies are required.
When assuming an incorporation of 7 vol.% of δ-Fe13 into the average mantle (PREM) at 1871-2771 km depth, − 0.6 to − 0.8 ± 0.2% anomaly of v S and + 0.3 ± 0.2% anomaly of v Φ would be generated, which are likely regarded as typical anomalies where the edge regions of LLSVPs have been sufficiently sampled. Although this discussion involves simplifications for the petrology and significant extrapolations, as well as noting that the magnitudes of anomalies in the seismic wave velocity and density depend on the seismic models 52,53 , our new data suggests that the presence of δ-(Al,Fe)OOH can produce seismic heterogeneities in the deep lower mantle. The total mass of LLSVPs is estimated to be 1.91 wt.% of the entire mantle 54 . On the other hand, the amount of  www.nature.com/scientificreports/ recycled oceanic crust that entered and stayed in the lower mantle since 4 Gyr is estimated as ~ 3-5 wt.% of the present mantle 55,56 , which is at least 1.5 times larger than the mass of LLSVPs. It has been suggested that about one ocean mass of water (0.037 wt.% of the mantle) has been transported via slabs to the deep mantle over the age of Earth 57 . If this amount of H 2 O is completely stored in LLSVPs, a hypothetical and unlikely end-member scenario, via a dense oxyhydroxide like the delta phase, the H 2 O content of LLSVPs is calculated to be 1.9 wt.%. This scenario suggests that a hydrous mafic system, such as the system examined in ref. 17 where the Al-rich Phase H (i.e., the Mg, Si, Fe-bearing δ-phase) could accumulate in LLSVPs. If this 1.9 wt.% H 2 O in the LLSVPs is completely stored in δ-Fe13, then this hydrous phase accounts for 12 vol.% of the LLSVP assemblage. This amount is larger than the 7 vol.% concentration of δ-Fe13 that can reproduce the typical anomalies at the edge regions of LLSVPs, as discussed above. The incorporation of 12 vol.% δ-Fe13 into "PREM" at 1871-2771 km depth can cause − 1.0 to − 1.4 ± 0.3% anomaly for v S and + 0.5 to + 0.6 ± 0.3% anomaly for v Φ , which are comparable to the anomalies reported in the central regions of LLSVPs 27 . If the temperature is higher than the phase stability, then the δ-phase will dissociate. If these conditions are met, then it is possible that released hydroxyl migrates upward and hydrates regions at lower temperatures 58 , although much work is still needed to understand the phase relations of such processes. Our elastic wave velocity results shed a new light on lower mantle seismic observations. At shallow lower mantle depths (800-900 km), δ-(Al,Fe)OOH exhibits low-v Φ and v P due to the spin transition, which may explain some seismic anomalies in this depth range. At deeper depths, the elastic behavior of δ-(Al,Fe)OOH is unlike typical lower mantle phases, and we discuss its relevance to LLSVP observations. Compositional differences between the LLSVPs and the surrounding mantle have been proposed as an origin of LLSVPs' seismic structure 47,[59][60][61] . Our study proposes that if δ-(Al,Fe)OOH is formed by the reaction between oceanic crust and water 17,18 , and accumulates in deep local regions such as the edge regions of LLSVPs, those regions can reproduce the seismically observed negative anomaly of v S and the positive anomaly of v Φ .

Methods
Single crystals of Fe-bearing δ-phase were synthesized by the high-pressure hydrothermal method using a 1000ton multi anvil apparatus (Hymag installed at Bayerisches Geoinstitut, University of Bayreuth), as previously reported by ref. 62 . The starting material of synthesis was a mixture of reagent-grade Al(OH) 3 (3) )OOH with dimensions of ~ 40 × 50 × 20 µm (Run 2, δ-Fe13-r2) were selected for this study. The chemical compositions of the two samples were determined using an electron microprobe operating in wavelength-dispersive mode (JEOL, JXA-8800 installed at Tohoku University) operated at 15 kV and 10 nA. The numbers in parentheses in the chemical compositions indicate one standard deviation of ten measurement data for each sample. The H 2 O contents were calculated from oxide deficit in total mass in the microprobe analyses, as reported in ref. 62 . The weight deficit analyses indicated that δ-Fe13-r1 and δ-Fe13-r2 contain 2.3 wt.% and 2.9 wt.% excess H 2 O, respectively, which are included in calculation of the density of the two samples.
High-pressure NRIXS experiments were performed using panoramic diamond anvil cells (pDACs). Beryllium disks are used as gaskets. The disks were pre-indented to thicknesses of 40-50 μm, and then 210 μm or 165 μm (in diameter) holes were drilled in the center of gaskets. The former and latter gaskets were used for the pDACs with 400 μm-culet anvils (Run 1, δ-Fe13-r1) and 300 μm-culet anvils (Run 2, δ-Fe13-r2), respectively. A mixture of 10-20 μm thick boron epoxy (amorphous boron powder:epoxy = 4:1 by weight) is attached to the gasket holes, which can stabilize the sample room at high pressures. We put the sample in the center of gasket hole and two ruby spheres beside the sample as pressure markers, and then loaded a compressed neon gas as a pressuretransmitting medium into the sample chamber of pDACs. The gas-loading was conducted at the Seismological Laboratory, California Institute of Technology. The cell design for the pDAC using 300 μm-culet anvils is same as that used in the previous synchrotron Mössbauer spectroscopy experiments by Ref. 22 .
We conducted NRIXS measurements at room temperature and high pressures up to 79 GPa at Sector 3-ID-B at the Advanced Photon Source. During the measurements, the storage ring was operated in 24-bunch top-up mode. Each bunch was separated by 153 ns. Incoherent inelastic X-ray scattering was detected with three avalanche photodiode detectors (APDs) positioned radially around the pDACs. In addition to the three APDs, another APD located at a downstream parallel to X-ray path was used for collecting forward elastic scattering, which helps constrain the resolution function. The energy of incident X-ray was tuned around 14.4125 keV, the energy of nuclear resonance of 57 Fe, by using a high-resolution monochrometer 63 with a typical energy resolution (full width at half maximum at zero energy transfer) measured to be about 1.1 meV. The details of the setup of NRIXS experiments at Sector 3-ID-B are presented elsewhere 33,63,64 . Here, we note a brief summary of conditions of the NRIXS experiments. Energy scans were performed over the ranges of − 100 to + 150 meV (ambient conditions in air), − 80 to + 110 meV (at 5.6 and 10.8 GPa), − 80 to + 120 meV (18.1, 21.1, 27, 32 and 37.1 GPa) and − 80 to + 150 meV (45,59,67 and 79 GPa), with an energy step-size of 0.25 meV. The pressure at each compression point was determined using the ruby fluorescence method 65 . The ruby fluorescence spectra were obtained before and after collecting NRIXS spectra at each pressure point using an off-line Raman spectrometer at the Sector 3-ID-B. Pressure errors arise from the standard deviation between the four ruby measurements obtained before and after NRIXS measurements for two ruby spheres, and the error of ruby scale 66 . The raw NRIXS data were analyzed by using the version 3.0.0 PHOENIX software package (www. NRIXS. com) to obtain the partial phonon density of states (pDOS), v D , and the Lamb-Mössbauer factor 33,67,68 . From the low-energy region of pDOS, v D is calculated from the following equation: www.nature.com/scientificreports/ where ρ is the density of sample, m is the mass of the nuclear resonant isotope, and D(E) indicates the pDOS in the low energy region. v(E) is equal to v D in the limit when energy (E) approaches zero. The pDOS shown in Fig. 2a are scaled by this Eq. (2). For a single fit range, a v D value is determined from the fitting using an empirical power law model (Eq. (1)) of the phonon dispersion curve, noting that 'Debye-like' dispersion is accommodated in this form. To determine a PDF for v D , binning of hundreds of possible fit ranges was used (Table S1). The reported v D and its uncertainty is defined as the peak position and FWHM of the fitting by asymmetric function in the probability distributions, respectively (Table S1). The details of this PDF method used in determination of v D was discussed in ref. 68 . The elastic wave velocities (v P , v S , and v Φ ) are then calculated using the K and ρ determined with XRD measurements 22 and the v D value, via the following equations.
Here, v P , v S , and v Φ are defined as in the following equations: The isothermal bulk modulus (K T ) obtained from the EoSs was converted to the adiabatic bulk modulus (K S ) by using the Grüneisen parameter (γ 0 = 0.64 ± 0.05), the exponent of its volume-independent (q = 1.8 ± 0.3) and Debye temperature (Θ 0 = 1485 ± 300 K) of δ-AlOOH reported by ref. 11 .
To calculate the velocities and density of the δ-Fe13 samples having natural isotopic enrichment of iron, we applied the following equations to the velocities and density determined with the 57 Fe-enriched experimental mass: where M nat amd M enr are the molecular mass of the samples having natural isotopic enrichment and 57 Fe-enrichment, respectively.

Data availability
The data that support the findings of this study are available from the corresponding author(s) upon reasonable request. (2) v(E) = mE 2 2π 2 3 ρD(E) 1 3 , www.nature.com/scientificreports/