Formation of large low shear velocity provinces through the decomposition of oxidized mantle

Large Low Shear Velocity Provinces (LLSVPs) in the lowermost mantle are key to understanding the chemical composition and thermal structure of the deep Earth, but their origins have long been debated. Bridgmanite, the most abundant lower-mantle mineral, can incorporate extensive amounts of iron (Fe) with effects on various geophysical properties. Here our high-pressure experiments and ab initio calculations reveal that a ferric-iron-rich bridgmanite coexists with an Fe-poor bridgmanite in the 90 mol% MgSiO3–10 mol% Fe2O3 system, rather than forming a homogeneous single phase. The Fe3+-rich bridgmanite has substantially lower velocities and a higher VP/VS ratio than MgSiO3 bridgmanite under lowermost-mantle conditions. Our modeling shows that the enrichment of Fe3+-rich bridgmanite in a pyrolitic composition can explain the observed features of the LLSVPs. The presence of Fe3+-rich materials within LLSVPs may have profound effects on the deep reservoirs of redox-sensitive elements and their isotopes.

T he large low shear velocity provinces (LLSVPs) are two massive and mysterious regions sitting beneath Africa and the Pacific 1-3 and occupy~3-9% of the volume of the Earth 3,4 . They are characterized by their lower-than-average seismic wave velocities 4 and extend by thousands of kilometers laterally and up to >1000 km vertically above the core-mantle boundary (CMB) 3,4 . As the largest seismic heterogeneities in the lower mantle, they may hold the key to understanding the thermal, chemical, and dynamical evolution of the Earth 5,6,7 . Different shear-wave tomography models 4 have reached agreement that the shear-wave velocity anomaly (dlnV S ) ranges from −0.5 to −1.0% in the shallow part of LLSVPs, while it could be up to −3.0% in the bottom part 6 . The compressional-wave tomography models also reveal negative anomalies of compressional-wave velocities (V P ) (refs. 8,9 ), although the amplitude, shape, and geographical location of V P anomaly vary widely among different models 5 . The V P anomaly generally has a smaller amplitude than the V S anomaly, causing a high dlnV S /dlnV P ratio 9 . In particular, waveform and travel-time seismic studies 1,3,10-12 reveal that the LLSVPs have sharp edges along their margins, which is also supported by the large lateral dV S gradients at the boundary of LLSVPs 6 .
The systematic and discontinuous contrasts in seismic properties indicate that the LLVSPs are likely composed of distinct chemical materials from the surrounding mantle 11,13 . A chemically distinct origin of the LLSVPs may be implicated by their density anomalies as well; however, large discrepancies for the density anomaly associated with the LLSVPs exist in the literature 5 . Recent tidal tomography based on body tide displacements 14 found that the mean density of the lower twothirds of the two LLSVPs is~0.5% higher than that of the surrounding mantle. On the contrary, a study using Stoneley modes suggested an overall negative density anomaly within LLSVPs, without excluding the possibility of a high-density anomaly within the lowermost LLSVPs 15 . It is still unknown whether the regional differences in density anomaly are caused by the choice of observations used to constrain density models or reflect the nature of LLSVPs associated with their origins.
Hypotheses for the origin of chemically distinct LLSVPs include processes associated with the accumulation of subducted oceanic crust over Earth history 16 and the differentiation and solidification of an ancient basal magma ocean 6 . Sunken piles of subducted oceanic crust, which is compositionally different from and significantly denser than the pyrolitic lower mantle 17 , was proposed to explain LLSVPs because of the low velocity of calcium silicate perovskite (CaSiO 3 , CaPv) (ref. 18 ). However, there are significant discrepancies in the velocity of CaPv between two experimental studies 18,19 and between experiments and theoretical results 20 . Elastic properties from ab initio simulations for the entire MORB assemblage indicate that subducted oceanic crust has relatively higher velocities than the ambient mantle 21 . Moreover, geodynamic simulations 22 suggested that the presentday subducted oceanic crust is too thin to provide enough negative buoyancy to survive viscous stirring and it hence is difficult to amass coherent thermochemical structures and shapes at the CMB similar to LLSVPs 23 . Alternatively, LLSVPs may be composed of primordial residues from basal magma ocean crystallization or core-mantle differentiation that have not yet been fully homogenized by the mantle convection [24][25][26] . These primordial materials would need to be intrinsically more dense than the surrounding mantle to overcome mantle stirring 27 . Dense Fe-Ni-S liquid, for instance, was proposed to explain the LLSVPs 28 , but the amount of this liquid remaining in the deep mantle, which depends on the drainage of melt to the core 29 , is under debate.
Consistent with both efficient drainage of metallic melt and a primordial origin of the LLSVPs is chemical heterogeneity produced by redox reactions in the magma ocean. Ferrous iron (Fe 2+ ) in silicate melts has been observed to disproportionate to ferric iron (Fe 3+ ) plus metallic iron (Fe 0 ) at high pressures 30 . Segregation of precipitated Fe 0 from the magma ocean into the core would enrich Fe 3+ in the mantle. Bridgmanite (Bdg), the dominant Fe 3+ -bearing mantle mineral, hosts Fe 3+ through the Fe 3+ -Fe 3+ or Fe 3+ -Al 3+ charge-coupled substitution in the deep mantle [31][32][33][34][35] . In particular, a Fe 3+ -rich Bdg with the chemical composition of (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 was recently synthesized by Liu et al. 36 . Oxidized domains enriched with such Fe 3+ -rich Bdg would be distinct from a pyrolitic lower mantle in sound velocity and density 37,38 and may be responsible for the origins of lower-mantle seismic and geochemical heterogeneities, such as the LLSVPs 1-3,6 . The conditions of formation of such Fe 3+ -rich Bdg in a mantle phase assemblage and its elastic properties at lower-mantle-relevant temperatures are vital to test this hypothesis, but these questions remain unclear.
In this work, we combine high pressure-temperature (P-T) experiments, ab initio calculations, and geodynamic simulations to study the formation of Fe 3+ -rich Bdg, its thermoelastic properties, and its dynamics in the mantle, to evaluate whether enrichment in Fe 3+ -rich Bdg can explain the seismic signatures of LLSVPs.

Results
The coexistence of Fe 3+ -rich and Fe-poor bridgmanite phases. We conducted a series of multi-anvil experiments with the bulk composition of 90 mol% MgSiO 3 -10 mol% Fe 2 O 3 from 10-24 GPa along a mantle geotherm (Supplementary Table 1). At 10 GPa and 1573 K (Fig. 1d) (Fig. 1c). This indicates that the iron-rich Aki coexisting with iron-depleted oxides/silicates are energetically more stable than a single-phase MgSiO 3 -Fe 2 O 3 solid solution with intermediate iron content. At 24 GPa and 1873 K, which corresponds to the P-T conditions around the uppermost lower mantle, the products consist of Fe-poor Bdg and Fe-rich Aki with 44-49 mol% Fe in both Mg and Si sites (Fig. 1a and  Supplementary Table 1), instead of forming a single phase of (Mg 0.9 Fe 0.1 )(Si 0.9 Fe 0.1 )O 3 Bdg. The Fe-rich Aki is evenly distributed in the matrix of the Fe-poor Bdg (Fig. 1a). The run products of experiments running at the same P-T conditions for 8 and 24 h have the same compositions within analytical uncertainty (Supplementary Table 1), confirming that the experiments reached equilibrium. In situ XRD measurements coupled with a diamond anvil cell (DAC) show that this Fe-rich Aki phase transforms to Bdg phase at 23.5 ± 1.0 GPa and 300 K (Supplementary Fig. 1). Moreover, this Fe 3+ -rich Bdg phase completely transforms back to Aki with the same lattice parameters as the starting Aki phase after the decompression of the DAC (Supplementary Fig. 1). The reversible phase transition of this Fe-rich phase means that for our multi-anvil experiments at the P-T conditions of the uppermost lower mantle, Fe 3+ -rich Bdg coexists with Fe-pool Bdg.
A previous multi-anvil study 39 synthesized Fe 3+ -only bridgmanite with 2-4 mol% Fe 3+ in the cation sites but did not observe the Fe 3+ -rich Bdg phase, possibly because the bulk Fe content of their experiments is not high enough to enable the formation of such Fe 3+ -rich Bdg. Moreover, the presence of unreacted MgO and SiO 2 in their run products ( Fig. 1 in ref. 39 ) suggests that their starting materials may not be homogenous or their experiments did not reach chemical equilibrium. Another study 40 synthesized Fe 3+ -only Bdg with the starting material of 90 mol% MgSiO 3 -10 mol% Fe 2 O 3 using laser-heated diamond anvil cell (LH-DAC) 40 . However, the chemical composition of Fe 3+ -only Bdg was not reported, possibly because there was a significant loss of Fe and Mg during melting 40 , and some Fe 3+ was reduced through reaction with diamond during laser heating 41 . In addition, the proportion of the Fe-rich Bdg phase is much smaller than the Fe-poor Bdg (Fig. 1), and therefore it is difficult to detect without a detailed analysis of the run products in ref. 40 .
We also performed ab initio calculations (see methods and supplementary materials) to investigate the stability of (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 Bdg under lower-mantle conditions. Our results show that the assemblage of (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 and MgSiO 3 Bdg has a lower Gibbs free energy than a singlephase (Mg 0.875 Fe 0.125 )(Si 0.875 Fe 0.125 )O 3 Bdg under the P-T of the whole lower mantle regardless of the spin state ( Supplementary  Fig. 2), indicating that the mixed two phases are more stable than the single-phase Bdg with a homogeneous composition. The theoretical results support our experimental observations and reveal that this Fe 3+ -rich Bdg with the chemical composition of approximately (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 should form as a separate phase coexisting with Fe-poor Bdg in the bulk composition of 90 mol% MgSiO 3 -10 mol% Fe 2 O 3 due to the miscibility gap.
Elastic properties and sound velocities of (Mg 0.5 Fe 0.5 ) (Si 0.5 Fe 0.5 )O 3 bridgmanite. Determining the seismic signature of a separate Fe 3+ -rich phase in equilibrium with the mantle phase assemblage requires elastic properties of this phase as a function of pressure and temperature conditions in the mantle. The elastic properties of (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 Bdg up to 130 GPa and 3000 K were theoretically obtained from ab initio calculations.
Because Bdg may also accommodate Al in the octahedral site 31,38,42 , we also conducted calculations on an end-member composition (Mg 0.5 Fe 0.5 )(Si 0.5 Al 0.5 )O 3 , to quantify the effect of Al on the elasticity of Bdg. In these calculations, the octahedral site (B-site) Fe 3+ in (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 Bdg undergoes a highspin (HS) state to a low-spin (LS) state transition with increasing pressure (Fig. 2a), while the dodecahedral-site (A-site) Fe 3+ in both compositions remain in the HS state throughout the lowermantle conditions 43 .
Our calculated volumes of HS-and LS-(Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 ) O 3 Bdg agree well with experimental measurements at 300 K (ref. 36 ) (Fig. 2a). The calculated spin transition of the B-site Fe 3+ occurs between 49 and 55 GPa at 300 K 36 , which is slightly higher and narrower than experimental results 36 . The predicted volume collapse (ΔV HS-LS ) caused by the spin transition of B-site Fe 3+ is 4.3% at 300 K, higher than experimental measurements (2.7%) 36 but consistent with previous theoretical calculations 44 on (Mg 0.5 Fe 0.125 )(Si 0.5 Fe 0.125 )O 3 Bdg assuming that ΔV HS-LS is linearly dependent on Fe 3+ content. The ΔV HS-LS discrepancy between experimental and theoretical studies is probably caused by the difference in the pressure range for the mixed-spin (MS) state. The predicted volumes of (Mg 0.5 Fe 0.5 )(Si 0.5 Al 0.5 )O 3 Bdg also show excellent agreement with experimental results at 300 K 42 . These comparisons demonstrate the high reliability of our DFT + U calculations in predicting elastic properties, as suggested by previous studies [43][44][45] .
The spin transition of B-site Fe 3+ in (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 Bdg generates a strong effect on bulk modulus (K S ) and V P , which both show deep valleys that broaden and decrease in magnitude with increasing temperature (Supplementary Fig. 3). The magnitude and width of the K S and V P anomalies are controlled by the fraction of LS B-site Fe 3+ (n LS ) and the pressure and temperature dependences of n LS (see Supplementary Materials).

Discussion
Combining elastic data from previous studies 20,45-47 with our results, we modeled the density and velocity anomalies caused by the presence of Fe 3+ -rich Bdg relative to the pyrolitic composition, which can effectively reproduce the reference seismic velocities and density of PREM 37,38 . The modeled chemical assemblage has pyrolitic mineral fractions (15% ferropericlase (Fp) + 78% Bdg + 7% CaPv) in which a portion of Fe 2+ -bearing Bdg was substituted by (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 Bdg. We find that the enrichment of (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 Bdg in the assemblage can explain the seismic features of the LLSVPs 4,6,8,9 .
The V S anomalies of −1.5% to −3.0% and the large dlnV S /dlnV P ratio >2.0 observed in LLSVPs 4,6,8,9 can be reproduced by the enrichment of 10-15% (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 Bdg in pyrolite at 110 GPa (Fig. 3). If LLSVPs are hotter (ΔT LLSVPS > 0), the required proportion of (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 Bdg would accordingly decrease; for example, it will decrease by~2% if ΔT LLSVPS is +400 K. For a pyrolitic composition, Bdg also con-tains~5 mol% Al 2 O 3 , which does not significantly change its velocities and density 48 . When 5 mol% Al 2 O 3 is incorporated into Bdg, the V P and V S anomalies caused by the presence of 15% (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.4 Al 0.1 )O 3 Bdg at ΔT LLSVPS equal to +400 K will be −1.5% and −3.1% (Fig. 3), respectively, which can also reproduce the dlnV S /dlnV P ratio >2.0. In addition, we find that the required proportion of (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 Bdg increases with the decreasing of Fe 2+ content in the modeled assemblage (noted by the FeO content in Fp, Fe 2+ Fp ) to explain the same V P and V S anomalies (Fig. 4). If there is no Fe 2+ , the (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 fraction that can reproduce the V S anomaly of −3.0% would be raised tõ 17% at ΔT LLSVPS equal to +400 K ( Fig. 4 and Supplementary  Fig. 5), which corresponds to a V P anomaly of~−1.0% and a dlnV S /dlnV P ratio of~3.0. The dlnV S /dlnV P ratio >2.0 tends to be reproduced at relatively low Fe 2+ contents (Fig. 4). If Fe 2+ Fp /Fe 2 + Fp, NM > 0.9 (Fe 2+ Fp, NM is the FeO content of Fp in a normal pyrolitic lower mantle, 18 mol%), the dlnV S /dlnV P ratio is less than 2.0. In contrast, when Fe 2+ Fp /Fe 2+ Fp, NM is <0.25, the dlnV S / dlnV P ratio >2.0 can be reproduced for different V S anomalies (Fig. 4). Our modeling implies that the enrichment of (Mg 0.5 Fe 0.5 ) (Si 0.5 Fe 0.5 )O 3 Bdg in the pyrolite assemblage with relatively lower Fe 2+ content can explain the velocity anomalies and the dlnV S / dlnV P ratio observed in LLSVPs 4,6,8,9 . In contrast to velocity anomalies, the modeled density anomaly (dlnρ) could be positive, zero, or negative, depending on the (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 abundance, temperature anomaly, and the fraction of Fe 2+ in Fp. To produce a V S anomaly of −3.0% for the assemblage enriched in (Mg 0.5 Fe 0.5 )(Si 0. 5 Supplementary Fig. 5). In general, the density anomaly is correlated with the magnitude of V S anomalies. For dlnV S < −1.0%, the dlnρ could be positive if ΔT LLSVPS is 0 K; however, if ΔT LLSVPS equals to +400 K, the dlnρ could be positive only when dlnV S is <−2.5% (Fig. 4). Our modeling suggests that the lowermost parts of LLSVPs with large negative velocity anomalies (<−3.0%) 4,6,8 could be denser than the ambient mantle, while the relatively shallow part of LLSVPs with −0.5 to −1.0% V S anomalies on average likely have slightly lighter density than the ambient mantle. Recent work proposed that the bottom two-thirds of the two LLSVPs are~0.5% denser than the surrounding mantle 14 , while Koelemeijer et al. 15 argued that the overall density of the LLSVPs is lower than the surrounding mantle. Such different conclusions may be related to different depth sensitivities of the datasets considered 49 , regardless of the input observations. Also, geodynamic modeling studies 50,51 suggested that the density anomalies of the LLSVPs relative to the surrounding mantle could be positive near the top and bottom of the LLSVPs but neutral or slightly negative in the middle of the LLSVPs. The density of LLSVPs could also be laterally inhomogeneous due to their internal convection and the entrainment of multiple compositional components into the LLSVPs 52 .
The present model for chemical heterogeneities within the LLSVPs is consistent with Fe-rich remnants of a basal magma ocean created early in Earth's history 6,[24][25][26] . Ferrous Fe in dense silicate melts associated with the basal magma ocean would partially disproportionate to Fe 3+ plus Fe 0 at high pressures 30 and segregation of precipitated Fe 0 from the magma ocean into core would enrich silicate melt in in Fe 3+ -rich bridgmanite piles. Our geodynamic modeling demonstrates that such Fe 3+ -rich piles with~18% (Mg 0.5 Fe 0.5 ) (Si 0.5 Fe 0.5 )O 3 Bdg, which is~1.5% intrinsically denser than the ambient mantle (Fig. 4c), could form large-scale thermochemical structures in the lowermost mantle without being mixed into the background mantle throughout Earth's history (Fig. 5), which may accumulate to form LLSVPs.
Our findings imply that the segregation of Fe 3+ -rich domains may cause heterogeneity in the redox states of the Earth's mantle, that is, the oxygen fugacity of the lowermost mantle may not be as low as inferred in previous studies 54,55 . We also expect these Fe 3+ -rich domains to be enriched in heavy iron isotopes because Fe 3+ has a larger Fe force constant than Fe 2+ (ref. 56 ), which may affect the iron isotopic features of the deep Earth 57 . The presence of lower-mantle oxidizing heterogeneities would have profound effects on the cycles of volatiles 30 and the deep reservoirs of redox-sensitive elements. For instance, the dense reduced Fe-C/H/S melts formed at the mantle transition zone and shallow  37 . The temperature anomaly is with respect to the normal mantle temperature from Brown and Shankland 74 . Data for elasticity at high pressure and temperature are derived from previous theoretical studies: Fp, ref. 47 ; Fe 2+ -Bdg, ref. 46 ; CaPv, ref. 20 . NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22185-1 ARTICLE lower mantle by slab/mantle interaction 58,59 , if they reach the lowermost mantle, could be converted back to an oxidized state instead of sinking into the core. Dynamic cycling with respect to mantle redox heterogeneity could provide new insights into the thermochemical evolution of the bulk silicate Earth and possibly the oxidation of the atmosphere.

Methods
High pressure-temperature experiments. The experiments were conducted using the 1000-ton multi-anvil apparatus at the University of Michigan. The COMPRES 8/3 and 10/5 cell assemblies were employed in the experiments. The starting material was a mixture of high purity (>99.99%) MgO, SiO 2 , and Fe 2 O 3 at a molar ratio of 9:9:1, which corresponds to a bulk composition of (Mg 0.9 Fe 0.1 )(Si 0.9 Fe 0.1 )O 3 . The mixture was heated at 1073 K overnight to remove the moisture and structural water before loading into a platinum capsule. The sample was compressed to target pressure and equilibrated at high temperature for 6-10 h to allow sufficient equilibrium. It was then quenched to room temperature and decompressed to 1 bar.
The recovered sample was polished, coated with carbon, and examined for texture and composition using the JOEL-7800FLV Scanning Electron Microprobe (SEM) and SX-100 Electron Microprobe Analysis (EPMA) at the Electron Microbeam Analysis Laboratory (EMAL) of the University of Michigan. An accelerating voltage of 15 kV and a beam current of 10 nA were employed for imaging and analysis. Forsterite and magnetite were used as standards for Mg, Si, and Fe quantification with EPMA.
First-principles calculations. Isothermal elastic tensors (C T ijkl ) of crystals in a Cartesian coordinate system usually can be calculated from Eq. (1) (ref. 60 ): where e ij (i,j = 1,3) are infinitesimal strains, P is the isotropic pressure, and F is the Helmholtz free energy, which can be expressed in the quasi-harmonic approximation (QHA) as: where V is the equilibrium volume of the crystal and T is temperature. Subscripts q and m refer to the phonon wave vector and the normal mode index, respectively. _ and k B are Planck and Boltzmann constants, and ω q,m is the vibrational frequency of the ith mode along with the wave vector q. U is the static energy at the equilibrium volume V. The second and third terms are the zero-point and vibrational energy contributions, respectively. Adiabatic elastic constants (C S ijkl ) can be derived from: where S is the entropy and C V is the constant volume heat capacity. Therefore, calculations of elastic tensors at high pressure and temperature using this usual method require a vibrational density of states (VDoS) of many strained configurations, which demand a tremendous amount of computational power to calculate based on the DFT. Wu and Wentzcovitch 61 proposed an analytical approach to calculate the thermal contribution to the elastic tensor, only requiring VDoS for unstrained configurations at different equilibrium volumes. This approach greatly reduces the computation cost to the level of <10% of the usual method without loss of accuracy. To obtain elastic tensors at static conditions and VDoS for unstrained configurations at different equilibrium volumes, we performed first-principles calculations using Quantum Espresso package 62 based on the DFT, plane wave, and pseudopotential. Local density approximation (LDA) was adopted for the exchange-correlation function. The energy cutoff for electronic wave functions was set as 70 Ry. The Mg pseudopotential was generated using the von Barth and Car method for all channels using a 2.5 Bohr cutoff radius and five configurations, 3s 2 3p 0 , 3s 1 3p 1 , 3s 1 3p 0.5 3d 0.5 , 3s 1 3p 0.5 , 3s 1 3d 1 , with weights of 1.5, 0.6, 0.3, 0.3, 0.2, respectively. The pseudopotentials for Si and O were generated using the Troullier-Martins method 63 with the cutoff radius of 1.47 Bohr for Si and 1.45 Bohr for O. Valence configurations for Si and O are 3s 2 3p 4 3d 0 and 2s 2 2p 4 , respectively. The pseudopotentials for Al and Fe were generated using the Vanderbilt method 64 with a valence configuration of 3s 2 3p 1 and a cutoff radius of 1.77 Bohr for Al, and a valence configuration of 3s 2 3p 6 3d 6.5 4s 1 4p 0 and a cutoff radius of 1.8 Bohr for Fe. To address the large on-site Coulomb interactions among the localized electrons (Fe 3d electrons) 65 , we introduced a Hubbard U correction to the LDA (LDA+U) for all DFT calculations. U values for Fe 3+ on A and B sites in bridgmanite were non-empirically determined using linear response method 66 in previous work 43 and adopted in this study. The initial structure was constructed by replacing one nearest-neighbor Mg 2+ -Si 4+ pair with one [Fe 3+ ] Mg -[Fe 3+ ] Si pair [43][44][45] . Crystal structures at variable pressures were well optimized on a 6 × 6 × 4 k-point mesh, and VDoS were calculated using the finite displacement method as implemented in the code PHONOPY 67 . The elastic tensors at static conditions were calculated from the linear dependence of stress on the small strain. Owing to the enormous computational cost of calculating the vibrational density of state based on LDA+U, we used the 20-atom unit cell for (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 bridgmanite [43][44][45] . Similar to previous studies on the elastic properties of Fe-bearing bridgmanite 45,46 , we only report results for aggregate elastic moduli, not individual elastic coefficients. The latter are sensitive to atomic configurations and therefore to supercell size, which can accommodate different configurations for the same composition. The aggregate elastic moduli, K S and G, are quite insensitive to the atomic configuration [68][69][70] .
The Helmholtz free energy calculated from Eq. (2) within the QHA versus volume was fitted by the isothermal third-order finite strain equation of state, and then we can obtain all thermodynamic properties, such as pressures at different temperatures and volumes. Our results show that (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 bridgmanite has a larger volume than (Mg 0.5 Fe 0.5 )(Si 0.5 Al 0.5 )O 3 bridgmanite (Fig. 2). The substitution of Al 3+ for LS Fe 3+ in the octahedral site causes a slight decrease of~1.0% in volume at >80 GPa. Compared to pristine Bdg (MgSiO 3 ), these two Fe 3+ -and Al 3+ -rich species have larger volumes; for example, at 90 GPa and 2000 K, the volumes of (Mg 0.5 Fe 0.5 )(Si 0.5 Al 0.5 )O 3 and (Mg 0.5 Fe 0.5 )(Si 0.5 Fe 0.5 )O 3 are 3.5% and 4.4% larger than that of MgSiO 3 , respectively.
Using the equation of states, we transferred volume-and temperaturedependent elasticity into pressure-and temperature-dependent elasticity. The adiabatic bulk modulus K S and shear modulus G can be obtained by computing the Voigt-Reuss-Hill averages 71 from elastic tensors. Thus, compressional and shear velocities can be calculated from the equations V P ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðKs þ 4 3 GÞ=ρ q and V S ¼ ffiffiffiffiffiffiffiffi G=ρ p (ρ is density). Bulk moduli (K S ), shear moduli (G), compressional-wave velocity (V P ), and shear-wave velocity (V S ) are also derived from LDA+U calculations as shown in Supplementary Fig. 3 Table 3).  ). In cases 1, 2, 3, and 4, the Rayleigh number is Ra = 1 × 10 7 , 1 × 10 6 , 1 × 10 8 , and 1 × 10 7 , respectively, and the iron-rich materials (shown by golden colors in the right panels) are 1.2%, 1.2%, 1.2%, and 1.5% intrinsically denser than the background mantle materials (shown by black colors in the right panels), respectively. Large thermochemical piles form after 4.5 Gyr for all cases.
Eqs. (3)(4)(5)(6) in Supplementary Materials. Thus, the Gibbs formation free energy of the decomposition reaction for pure HS/LS state can be expressed as: where H stat is the internal energy or the Gibbs free energy without the vibrational contribution.
As shown in Supplementary Fig. 2 In order to check the effect of the exchange-correlation function on the results, we also calculated the enthalpy change of this reaction using the generalized gradient approximation with Hubbard U correction (GGA+U). U values for Fe 3+ on A and B sites in bridgmanite are 3.3 and 4.5 eV, respectively. The ΔH predicted by the GGA+U calculations is similar to that from the LDA+U calculations ( Supplementary Fig. 2), both of which favor the mixed phases over the single phase.
Thermoelastic models for the estimations of velocity and density heterogeneities. Previous studies 37  In other words, the ferropericlase and Ca-perovskite contents are fixed to 15% and 7%, respectively. We also explore the effect of Fe 2+ content in the assemblage (noted by the Fe 2+ content in Fp, Fe 2+ Fp ) on modeling results because the incorporation of Fe 2+ into Fp and Bdg decreases their velocities to some extent 46,47  The elastic moduli and densities of the aggregate are calculated using: where ρ i , M i , and f i are the density, bulk modulus (K S ) or shear modulus (G), and the fraction of the ith mineral, respectively. Similarly, the compressional and shear velocities (V P and V S ) were derived from V P ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðKs þ 4 3 GÞ=ρ q and V S ¼ ffiffiffiffiffiffiffiffi G=ρ p , and hence the velocity and density anomalies relative to the pyrolitic lower mantle are estimated with the consideration of temperature anomaly.
Geodynamic models. We performed thermochemical calculations to study the dynamics of iron-rich materials in the lowermost mantle. The numerical simulations were conducted by solving the nondimensional equations of conservation of mass, momentum, and energy under the Boussinesq approximation. The intrinsic density anomaly is represented by the buoyancy number B, which is defined as the ratio between the intrinsic (compositional) density anomaly and the density anomaly caused by thermal expansion, B = Δρ/(ραΔT), where Δρ is the intrinsic density anomaly for the iron-rich materials compared to the background mantle. We used α = 1 × 10 −5 K −1 and ΔT = 2500 K in this study.
The whole-mantle dynamics is simulated in a 2D Cartesian box with an aspect ratio of 3:1. The model domain contains 1536 and 512 elements in the horizontal and vertical directions, respectively. Models are basally heated, with the top and bottom having a fixed temperature at T = 0 and T = 1, respectively. The top and bottom boundaries are both free-flip and the side boundaries are reflective. The viscosity is determined by η = η 0 exp[0.5−T], where η 0 is a viscosity pre-factor, and A is the activation energy. Here, η 0 is 1.0 and 30.0 for the upper mantle and the lower mantle, respectively, and A = 9.21 which leads to the viscosity changing four orders of magnitude as temperature increases from 0 to 1. Initially, the whole mantle was assumed to be hot with a temperature of T = 0.72 (or 1800 K) everywhere, and we introduced a global layer of iron-rich materials in the lowermost 300 km of the mantle.
We performed four cases. Case 1 is the reference case in which the iron-rich materials are 1.2% intrinsically denser (with a buoyancy number of B = 0.48) than the background mantle materials and the Rayleigh number, which controls the vigor of mantle convection, is Ra = 1 × 10 7 . Cases 2 and 3 have a Rayleigh number of Ra = 1 × 10 6 and Ra = 1 × 10 8 , respectively, while other parameters are the same as case 1. In case 4, the iron-rich materials are 1.5% intrinsically denser (with a buoyancy number of B = 0.6) than the background mantle materials, and other parameters are the same as case 1.

Data availability
The data are available in the main text, the supplementary materials, and from the corresponding authors.

Code availability
The open-source Quantum Espresso package used in this study is available at https://www.quantum-espresso.org/.