Iron-rich carbonates stabilized by magnetic entropy at lower mantle conditions

Constraining the ﬂux of carbon in and out of the interior of the Earth due to long-term geological processes is important, because of the inﬂuence that it has on climate change. On timescales of billions of years, host minerals such as carbonate phases could play a signiﬁcant role in the global carbon cycle, transporting carbon into the lower mantle as a component of subducting slabs. We use density functional theory based calculations to study the high-pressure, high-temperature phase stability of Mg 1- x Fe x CO 3 . Our results show that iron-rich phases, where carbon is in tetrahedral coordination, are only stable at lower mantle conditions due to their magnetic entropy, which is also responsible for the unusual shape of their phase boundary. Low-pressure carbonate phases are found to be highly anisotropic, but high-pressure carbonate phases are not, which has important implications for their seismic detectability. Our work conﬁrms that future discussions of the global carbon cycle should include the deep Earth.


Introduction
The global carbon cycle is of great importance, due to its influence on climate change (Dasgupta and Hirschmann, 2010). It describes the distribution and exchange of carbon between major reservoirs, such as the atmosphere, crust, mantle and core. Studies of the global carbon cycle have often focused primarily on the atmosphere, oceans, and shallow crustal environments, referred to as the "near-surface cycle" (Hazen and Schiffries, 2013). While this works well for short time scales, to predict long-term changes in the concentration of CO 2 in the atmosphere it is necessary to include exchange between the surface and interior of the Earth. Carbon enters the interior of the Earth via subduction of carbonatebearing slabs (Dasgupta and Hirschmann, 2010). Carbonates are believed to persist into the mantle, provided that a slab follows a very cold subduction geotherm (Syracuse et al., 2010), which precludes their melting (Kakizawa et al., 2015;Thomson et al., 2016) or reaction with SiO 2 (Drewitt et al., 2019;Kakizawa et al., 2015;Maeda et al., 2017). Low solubility of carbon in mantle minerals, suggests that, once in the mantle, carbon is stored as either carbonates or diamond (Panero and Kabbes, 2008;Shcheka et al., 2006), depending on the local oxidation state (e.g. Stagno and carbon is in tetrahedral coordination (Boulard et al., 2015(Boulard et al., , 2012(Boulard et al., , 2011Merlini et al., 2015). Three out of four studies report (C 3 O 9 ) 6− rings (Boulard et al., 2015(Boulard et al., , 2012(Boulard et al., , 2011. Experimental investigations of pure FeCO 3 report that the R3c structure is stable up to, at least, 130 GPa at 300 K, with a pressure-induced spin transition occurring at about 45 GPa (Cerantola et al., 2015;Farfan et al., 2012;Lavina et al., 2010bLavina et al., , 2009Mattila et al., 2007;Weis et al., 2017). In contrast, at 1500 K and above, Boulard et al. (2012) observed a high-pressure phase coexisting with other run products above about 40 GPa. This was assigned a Fe 4 C 3 O 12 composition, but the atomic positions were unresolved. Liu et al. (2015) observed a high-pressure phase at similar conditions, above about 50 GPa, at 1400 K, but reported it to have a FeCO 3 composition. Structure refinement showed that the Pmm2 space group best fit their X-ray diffraction pattern, but the atomic positions were unresolved. In a more recent study by Cerantola et al. (2017), FeCO 3 was found to break down above 70 GPa at 1400 K, forming a complex series of decomposition products as temperature and pressure are increased. The first of these was found to have the same chemical composition (Fe 4 C 3 O 12 ) as the high-pressure phase of Boulard et al. (2012) and a structure consistent with the X-ray diffraction pattern of the high-pressure phase reported by Liu et al. (2015), suggesting that they are the same. Furthermore, Cerantola et al. (2017) were able to resolve the crystal structure, finding it to be a tetrairon(III) orthocarbonate, containing CO 4 tetrahedral units. Despite apparent observation of the same high-pressure structure, the phase diagrams of the investigations differ somewhat. For example, Liu et al. (2015) observe the formation of their high-pressure phase at 50 GPa, while at this pressure Cerantola et al. (2017) still find FeCO 3 (R3c) to be stable. Between 70-120 GPa and 1500-2200 K, Liu et al. (2015) only observe FeCO 3 (R3c) and their high-pressure phase, whereas Boulard et al. (2012) observe Fe 4 C 3 O 12 co-existing with diamond and iron oxides, and Cerantola et al. (2017) observe various combinations of co-existing Fe 4 C 3 O 12 , Fe 4 C 4 O 13 and iron oxides. In the latter two studies the phase boundary was also poorly constrained.
In the present work, we perform density functional theory (DFT) calculations to examine the phase relations of Mg 1-x Fe x CO 3 at lower mantle conditions. Our results illustrate the importance of including temperature when investigating phase behaviour, with significant differences found in the results of 0 K and high temperature calculations. In particular, our results indicate that magnetic entropy plays a significant role in stabilizing Fe 4 C 3 O 12 + C (diamond) and is the reason for the unusual shape of its phase boundary. For lower iron concentrations, our results suggest a phase transition from the R3c to C 2/m structure at about 75 GPa, with a small binary phase loop of a few GPa. In addition, we find the calculated seismic anisotropy of the high-pressure phases to be much smaller than that of the low-pressure phases, which has implications for their potential seismic detectability.

First-principles calculations
First-principles calculations were performed using VASP Furthmüller, 1996a, 1996b), employing the projector augmented wave (PAW) method (Blöchl, 1994;Kresse and Joubert, 1999), within the framework of density functional theory. For most calculations the PBE exchange-correlation functional was used (Perdew et al., 1996), but we also performed some calculations using a modified version of the HSE06 exchange-correlation functional (Krukau et al., 2006), to study reactions containing both Fe 2+ and Fe 3+ (discussed below). The valence electron configurations for the potentials were 2p 6 3s 2 for Mg, 3d 7 4s 1 for Fe, 2s 2 2p 2 for C, and 2s 2 2p 4 for O. The kinetic-energy cut-off for the plane-wave basis set was set to 850 eV. For calculations of lattice parameters and internal energies unit cells were used and the Brillouin zone sampled using the following Monkhorst-Pack grids (Monkhorst and Pack, 1976): 6 ×6 × 6 for the R3c structures of MgCO 3 and FeCO 3 (10 atoms), 2 × 2 × 2 for the C 2/m, P 2 1 2 1 2 1 , P 2 1 /c and P -1 structures of MgCO 3 and FeCO 3 (30-60 atoms), 2 × 2 × 2 for Fe 4 C 3 O 12 (38 atoms (primitive cell)) and 6 × 6 ×6 for diamond (2 atoms). The break condition for the electronic self-consistent loop and ionic relaxation were 10 −6 and 10 −5 eV, respectively. These parameters ensured that energies were converged to within 1 meV/atom and elastic constants to within a few percent.
It is well-known that standard density functional theory can fail to predict the correct band structure of transition metal minerals and oxides, because of the strongly correlated d electrons involved (Anisimov et al., 1991). In order to accurately describe the properties of FeCO 3 and Fe 4 C 3 O 12 at high pressure, we employed the DFT+U method (Anisimov et al., 1997), in particular, the simplified scheme of Dudarev et al. (1998) in which only the difference between onsite Coulomb interaction parameter U and onsite exchange parameter J is meaningful. Throughout the present work, U is used to mean U-J. In the present study, U = 2 eV was found to give best agreement with experimental values for the spin transition pressure of iron and so was adopted for production calculations of FeCO 3 and Fe 4 C 3 O 12 .
The DFT+U method implemented in VASP uses a constant U, which is inappropriate for studying reactions involving Fe 2+ and Fe 3+ , as iron in different oxidation states requires a different value of U, e.g. in (Mg,Fe)SiO 3 post-perovskite the value of U for Fe 3+ is about 1 eV higher than that for (Fe 2+ ) (Yu et al., 2012). Hybrid exchange-correlation functionals offer an alternative, although at a higher computational cost. In view of this, for all calculations of phase boundaries that involved both Mg 1-x Fe x CO 3 and Fe 4 C 3 O 12 , we used a modified version of the HSE06 exchange-correlation functional (Krukau et al., 2006) to compute internal energies.
The modification involved decreasing the fraction of Hartree-Fock exchange from 25 to 10 percent, since a recent investigation showed that, for FeCO 3 , this leads to improved agreement with experimental observations (Sherman, 2009). The internal energies were combined with vibrational free energies calculated using the DFT+U method (see below) to determine their Gibbs free energies. Tests showed that the difference in the vibrational free energies calculated using the DFT+U method and modified HSE06 exchange-correlation functional was no more than 10 meV/atom.
For minerals containing iron, additional terms need to be included in the calculation of the Gibbs free energy (e.g. Tsuchiya et al., 2006), due to the magnetic and configurational entropy contributions. The fraction of iron in the low-spin state can be computed as where G LS-HS (P , T ) is the calculated difference in the Gibbs free energy of the low-spin and high-spin states, at a particular pressure and temperature. The total volume is where V HS is volume of the high-spin state and V LS the volume of the low-spin state.
In a spin-crossover region, high-spin and low-spin iron is treated as a solid solution with configurational entropy where k B is Boltzmann's constant and n the fraction of iron in the low-spin state. The total Gibbs free energy of a phase is thus where P is pressure, T is temperature, G HS (P , T ) the Gibbs free energy of the high-spin state, G LS (P , T ) the Gibbs free energy of low-spin state, n the fraction of iron in the low-spin state, S mag the magnetic entropy and S conf the configurational entropy.
For phase transitions in pure end-members, phase boundaries were calculated from the difference in their Gibbs free energies, at a given pressure and temperature. For solid solutions they were determined from the co-tangent of the Gibbs free energies, in a similar manner to the method reported by (Metsue and Tsuchiya, 2012), making the assumption that we have ideal mixing and that the spin transition pressure is independent of composition (Fu et al., 2017;Hsu and Huang, 2016;Lavina et al., 2010a;

Elastic and seismic properties
In order to determine the six independent elastic constants of MgCO 3 (R3c) and FeCO 3 (R3c), one rhombohedral (for c 13 and c 33 ) and two monoclinic strains (for c 11 , c 12 , c 14 and c 44 ) were applied to the optimized structures. Corresponding stresses were computed by relaxing atomic positions in the strained configurations, with the lattice constants kept fixed. Elastic constants were then calculated from simple stress-strain relations, by the fitting of a second-order polynomial.
It should be noted that MgCO 3 (R3c) and FeCO 3 (R3c) have two settings. One has 10 atoms in the unit cell (rhombohedral setting) and another has 30 atoms (conventional setting). For computational efficiency, the rhombohedral setting was used in all calculations. While the total energy is unchanged under coordinate transformations, elastic constants do change. For the convenience of comparing with previous results, the following cartesian components are used: where a R is the length of arbitrary axis (a R , b R or c R ) since they are equal, R means this setting is rhombohedral and θ is the angle between two axes. The relation between lattice parameters in rhombohedral and conventional setting is expressed as where a c and c c are the lengths of the a and c axes in the conventional setting, respectively. The choices of crystal orientation affects the sign of c 14 (Golesorkhtabar et al., 2013). In the present study, the calculated value of c 14 is negative, but in order to compare with previous results, it is reported positive. The elastic constants of MgCO 3 (C 2/m) were calculated from six triclinic strains. The thirteen independent elastic constants were divided into six groups (c 11 , c 12 , c 13 , c 15 ; c 22 , c 23 , c 25 ; c 33 , c 35 ; c 44 , c 46 ; c 55 ; and c 66 ), with each group being calculated from one strain.
The elastic constants of Fe 4 C 3 O 12 were calculated using the same strains used for MgCO 3 (R3c) and FeCO 3 (R3c), as it has similar R3c symmetry.
The anisotropy factor for V P ( A P ) is defined as where V P , max and V P , min are the maximum and minimum compressional wave velocities. The analogous polarization anisotropy where V S1 and V S2 are the shear wave velocities and V S the aggregate shear-wave velocity.

Pressure correction
It is well known that GGA exchange-correlation functionals, like PBE, overestimate pressure. In order, to account for this, an empirical pressure correction (Oganov et al., 2001), of a few GPa is applied to all of our results. A single pressure correction was   is shown as a grey line, with the dashed portion depicting its extrapolation to high pressure. Solid white circles indicate conditions where experimental studies report MgCO 3 (R3c) to be stable (Isshiki et al., 2004). Half-filled circles (Isshiki et al., 2004) and squares (Boulard et al., 2011) indicate conditions where a highpressure phase was observed. The black dashed line is a mantle geotherm (Stixrude and Lithgow-Bertelloni, 2011). calculated for all MgCO 3 phases and all FeCO 3 phases, based on experimental data for the R3c structure. For more details please see Supplementary Material Table S1.

Phase diagram of MgCO 3
Our computed phase diagram for MgCO 3 , based on the calculated Gibbs free energies of the R3c, P -1, C 2/m, P 2 1 /c and P 2 1 2 1 2 1 phases is shown in Fig. 1. In agreement with the work of Zhang et al. (2018), we find a direct transition from the R3c structure to the C 2/m structure at about 75 GPa at all temperatures above 1000 K, with the P -1 structure being stable over a short pressure range at lower temperature. The phase boundary at 75 GPa is in agreement with the predictions of Oganov et al. (2008) and observations of Boulard et al. (2011) and Maeda et al. (2017), but conflicts with the work of Isshiki et al. (2004), who found the R3c phase to be stable above 100 GPa at 2000 K. In contrast to Zhang et al. (2018), we find the P 2 1 /c structure to be stable in a small region at about 130 GPa, below about 500 K. However, the difference in the Gibbs free energies of the P 2 1 /c and C 2/m structures is only a few meV/atom at these conditions. The P 2 1 2 1 2 1 phase is found to be stable at pressures higher than those expected in the lower mantle, as reported by others (Pickard and Needs, 2015;Zhang et al., 2018).
Comparison of the sequence of phase transitions predicted at lower temperature, with that expected along a geotherm (Stixrude and Lithgow-Bertelloni, 2011) reveals some significant differences. In particular, the P -1 and P 2 1 /c structures, which are stable at 0 K, are unstable above about 700 K. This shows the importance of including temperature in studies of phase stability. Our calculations indicate that it is the C 2/m phase that is stable in the lowermost mantle. The C 2/m structure comprises (C 3 O 9 ) 6− rings, similar to those reported in high-pressure phases of iron-rich compositions (Boulard et al., 2015(Boulard et al., , 2012(Boulard et al., , 2011, discussed next.

Spin state of iron in FeCO 3 and Fe 4 C 3 O 12
In order to investigate the phase diagram of FeCO 3 it is essential to first determine the correct spin state of iron in the various phases being considered. In the present work, we consider FeCO 3 (R3c) and the Fe 4 C 3 O 12 phase reported by Cerantola et al. (2017). In addition, we consider FeCO 3 (C 2/m), as the iron end-member of the most stable MgCO 3 phase at lower mantle conditions.
Our calculated 0 K enthalpy differences between high-spin and low-spin FeCO 3 (R3c) are shown in Fig. 2. Calculations were performed with the PBE exchange-correlation functional (Perdew et al., 1996) and with the addition of a Hubbard U correction (Anisimov et al., 1997;Dudarev et al., 1998) (DFT+U), using effective U values of 2 eV and 4 eV. At 0 GPa, both PBE, and DFT+U (U = 2 eV) predict a ferromagnetic ground state, while DFT+U (U = 4 eV) predicts an antiferromagnetic ground state. The true ground state is antiferromagnetic (Badaut et al., 2010). To obtain an antiferromagnetic ground state using standard DFT functionals, spin-orbit coupling (SOC) must be included in the calculations (Badaut et al., 2010). Since the enthalpy difference between the ferromagnetic and antiferromagnetic states is small (3 meV/atom at 0 GPa) and physical properties of the ferromagnetic and antiferromagnetic structures are almost identical (Supplementary Material Table S2), neglect of SOC should not affect our results and is therefore not included in our calculations.
Using PBE we calculate a high-spin to low-spin transition pressure of about 20 GPa (Fig. 2). Inclusion of U shifts the spin transition to higher pressures, about 40 GPa for U = 2 eV and 70 GPa for U = 4 eV, as expected. The difference in the spin transition pressures calculated for ferromagnetic and antiferromagnetic structures is, at most, 2 GPa. Experimental and theoretical studies of FeCO 3 (R3c), report a spin transition pressure of about 45 GPa (Cerantola et al., 2015;Farfan et al., 2012;Hsu and Huang, 2016;Lavina et al., 2010bLavina et al., , 2009Mattila et al., 2007;Weis et al., 2017). Studies of a wide range of compositions of Mg 1-x Fe x CO 3 (x = 0.05-1) report similar spin transition pressures, in the range 40-52 GPa (Fu et al., 2017;Hsu and Huang, 2016;Lavina et al., 2010a;Lin et al., 2012;Liu et al., 2014;Merlini and Hanfland, 2013;Spivak et al., 2014), confirming suggestions that the concentration of iron has little effect on spin transition pressure. Shi et al. (2008) report a spin transition pressure of 28 GPa for FeCO 3 (R3c), using DFT+U (U = 4 eV), which is much smaller than our value of 70 GPa. The reason for this difference is not clear. Farfan et al. (2012) report a spin transition pressure of 44 GPa using PBE, which is much higher than our value of 20 GPa. This could be because they allowed spin to vary freely during their structural relaxations, allowing the magnetic moment to decrease gradually as the spin transition pressure is reached, before undergoing magnetic collapse. Temperature shifts the spin transition to higher pressure and leads to a mixed-spin region, where high-spin and low-spin iron co-exist (Hsu and Huang, 2016;Liu et al., 2014) (Supplementary Material Fig. S7). Our calculated spin transition pressure at 300 K (defined as the pressure at which n = 0.5, where n = fraction of low-spin iron) is about 44 GPa, in good agreement with previous studies (Cerantola et al., 2015;Farfan et al., 2012;Hsu and Huang, 2016;Lavina et al., 2010bLavina et al., , 2009Mattila et al., 2007;Weis et al., 2017). The width of the mixed-spin region (defined as the pressure range over which 0.05 < n < 0.95) is about 4 GPa at 300 K and about 13 GPa at 1200 K, in good agreement with Liu et al. (2014), who report a width of about 4 GPa at 300 K and about 10 GPa at 1200 K. Since DFT+U (U = 2 eV) leads to good agreement with experiment, we adopt this value for all calculations of thermodynamic and elastic and seismic properties.
Calculated 0 K enthalpies for the high-spin and low-spin states of Fe 4 C 3 O 12 [+ C (diamond)] and FeCO 3 (C 2/m) (Supplementary Material Fig. S8 and Table S3), indicate that they are in a highspin state at lower mantle pressures. Since spin transition pressure only increases with temperature, no spin transition is expected in these phases at lower mantle conditions. This supports the results of Cerantola et al. (2017) for Fe 4 C 3 O 12 , but disagrees with those of Liu et al. (2015). For Fe 4 C 3 O 12 , a mixed-spin state is more stable than the low-spin state at all pressures. In this mixed-spin state, iron atoms situated on the threefold axis are in a high-spin state, while all others are in a low-spin state.

Phase diagram of FeCO 3
Liu et al. (2015) observed a transition from FeCO 3 (R3c) structure to a high-pressure phase at about 50 GPa and 1400 K. They were able to resolve the space group of this new phase, but not atomic positions. However, more recent experimental work  indicates that the high pressure phase of FeCO 3 reported by Liu et al. (2015) is Fe 4 C 3 O 12 . Cerantola et al. (2017) report a complex series of decomposition products for FeCO 3 (R3c) between 70 GPa and 110 GPa, as temperature is increased. The first of these is Fe 4 C 3 O 12 , observed as a single phase at about 1400 K and with other phases up to about 2250 K. The formation of Fe 4 C 3 O 12 can be described by the reaction: FeCO 3 (R3c) = Fe 4 C 3 O 12 + C (diamond), which we investigate here.
Calculated 0 K enthalpies for the high-spin and low-spin states of FeCO 3 (R3c), Fe 4 C 3 O 12 + C (diamond) and FeCO 3 (C 2/m) (Supplementary Material Fig. S8 and Table S3), indicate that at 0 K low-spin FeCO 3 (R3c) is the most stable phase, at lower mantle pressures. Below about 20 GPa the Fe 4 C 3 O 12 structure reported by Cerantola et al. (2017) is found to undergo an iso-symmetric phase transition to a different structure identified by using FIND-SYM (Stokes and Hatch, 2005) (Supplementary Material Table S4), but it is less stable than high-spin FeCO 3 (R3c) in this pressure range. The situation is different at lower mantle temperatures.
Our calculated phase boundary between FeCO 3 (R3c) and Fe 4 C 3 O 12 + C (diamond), indicates that Fe 4 C 3 O 12 + C (diamond) is more stable than FeCO 3 (R3c) at lowermost mantle conditions (Fig. 3). The phase boundary was calculated using both the DFT+U method and a modified version of the HSE06 exchange-correlation functional to compute internal energies. The latter should be more reliable for reactions where there is a change in oxidation state, as it does not require specification of a U value, which may vary with oxidation state. The phase boundary calculated using the modified version of the HSE06 exchange-correlation functional fits the experimental data well, marking the edge of the pressuretemperature space where high-pressure phases have been ob- served Liu et al., 2015). The good agreement with the phase boundary of Liu et al. (2015), supports the idea that the phase they observed was Fe 4 C 3 O 12. The phase boundary calculated using DFT+U, also match reasonably well, but is about 500 K lower than experiments. The spin transition in FeCO 3 (R3c), calculated using the modified HSE06 exchange-correlation functional, is about 10 GPa lower than that calculated with DFT+U and experiment (Cerantola et al., 2015;Farfan et al., 2012;Lavina et al., 2010bLavina et al., , 2009Mattila et al., 2007;Weis et al., 2017), highlighting that it is not a perfect approximation. However, both methods predict the formation of Fe 4 C 3 O 12 + C (diamond) at about 50 GPa, along the mantle geotherm. The unusual shape of the phase boundary and the driving force for the formation of Fe 4 C 3 O 12 is magnetic entropy (Fig. 4). At pressures above the spin transition, FeCO 3 (R3c) is in a low-spin state, which means that its magnetic entropy is zero regardless of the temperature. On the other hand, Fe 4 C 3 O 12 is in a high-spin state at all pressures considered, which means its magnetic entropy increases with temperature, eventually stabilizing it with respect to FeCO 3 (R3c), even though it has a higher lattice enthalpy. This is seen at 120 GPa, where the formation of Fe 4 C 3 O 12 + C (diamond) occurs at about 1500 K if magnetic entropy is accounted for, but does not occur up to 3000 K if it is neglected (Fig. 4(a)). The same thing is also observed for FeCO 3 (C 2/m). The magnetic entropy of FeCO 3 (C 2/m) increases faster than that of Fe 4 C 3 O 12 + C (diamond), since it contains ferrous iron rather than ferric iron, and it becomes more stable above 3000 K. The phase boundary between FeCO 3 (R3c) and FeCO 3 (C 2/m) has a similar form to that of Fe 4 C 3 O 12 , but is about 500 K higher. Following a mantle geotherm, Fe 4 C 3 O 12 + C (diamond) becomes stable following the onset of the spin transition in FeCO 3 (R3c), where it loses its magnetic entropy ( Fig. 4(b)). The unusual importance of magnetic entropy arises due to the high concentration of iron in the phases and the high temperatures in the mantle. Cerantola et al. (2017) report the co-existence of Fe 4 C 4 O 13 with Fe 4 C 3 O 12 above 70 GP and 1600 K, and as a single phase above  (Stixrude and Lithgow-Bertelloni, 2011). The dashed lines show values without the magnetic entropy contribution (shaded region). FeCO 3 (R3c) has no magnetic entropy at 120 GPa, since iron is in a low-spin state. Following a mantle geotherm, the magnetic entropy of FeCO 3 (R3c) decreases as the spin transition progresses, reaching zero at about 80 GPa. In contrast, the magnetic entropy of FeCO 3 (C 2/m) and Fe 4 C 3 O 12 + C (diamond) increases with temperature, even at high pressure, as iron remains in a high-spin state. This means that FeCO 3 (C 2/m) and Fe 4 C 3 O 12 + C (diamond) are more stable than FeCO 3 (R3c) at high temperature, at pressures above the spin transition.  (Isshiki et al., 2004). Half-filled circles (Isshiki et al., 2004) and squares (Boulard et al., 2011) indicate conditions where a high-pressure phase was observed. The black dashed line is a mantle geotherm (Stixrude and Lithgow-Bertelloni, 2011). 3000 K, but it was not observed by Liu et al. (2015), which suggests that, perhaps, different starting materials or experimental procedures can lead to the formation of different high-pressure phases. Determining the stability of the Fe 4 C 4 O 13 phase is out of the scope of the present work.

Phase diagram of Mg 1-x Fe x CO 3
Previous investigations of high-pressure phases of Mg 1-x Fe x CO 3 have tended to focus on iron-rich compositions. While these may exist in some parts of the mantle where iron enrichment has taken place, it is expected that in eclogite 0.16 < x < 0.21 and in peridotite x = 0.07 (Dasgupta et al., 2004;Dasgupta and Hirschmann, 2006;Sanchez-Valle et al., 2011). Using our calculated Gibbs free energies for the magnesium and iron end-members of the R3c and C 2/m phases, we determine the binary phase loop for the R3c to C 2/m transition for (Mg 0.93 Fe 0.07 )CO 3 (Supplementary Material Fig. S9) and (Mg 0.81 Fe 0.19 )CO 3 (Fig. 5). These show that iron has a small effect, reducing the R3c to C 2/m phase transition pressure by 3-5 GPa and opening up a binary phase loop of 3-5 GPa. One would, therefore, expect the phase transition in natural ironbearing carbonates, which have not undergone iron enrichment, to  , while solid white squares indicate conditions where Mg 0.6 Fe 0.4 CO 3 (R3c) is stable (Boulard et al., 2012). Half-filled circles indicate conditions where a high-pressure phase, that was later identified as Fe 4 C 3 O 12  is stable, while half-filled squares indicate conditions where a high-pressure phase of Mg 0.6 Fe 0.4 CO 3 consistent with a C 2/m structure is stable (Boulard et al., 2012). The black dashed line is a mantle geotherm (Stixrude and Lithgow-Bertelloni, 2011). occur at a similar depth in the lower mantle as that for MgCO 3 (Fig. 1).
To compare our results with experimental studies of iron-rich compositions, we calculate the Gibbs free energies of Mg 0.35 Fe 0.65 CO 3 (R3c) and Mg 0.35 Fe 0.65 CO 3 (C 2/m) and consider their potential transformation to (0.35 MgCO 3 (C 2/m) + 0.1625 (Fe 4 C 3 O 12 + C (diamond))). For simplicity, the phase loop between the two Mg 0.35 Fe 0.65 CO 3 phases, which is only about 5 GPa, is neglected. The calculated boundary between the phases is plotted in Fig. 6. Our results indicate that above about 60 GPa and 1800 K, Mg 0.35 Fe 0.65 CO 3 (R3c) transforms to Mg 0.35 Fe 0.65 CO 3 (C 2/m), whereas at lower temperatures it transforms to MgCO 3 (C 2/m), Fe 4 C 3 O 12 and C (diamond). In their study of an Mg 0.35 Fe 0.65 CO 3 composition, Liu et al. (2015) reported the presence of Siderite-II, later identified as Fe 4 C 3 O 12 , close to our predicted phase boundary. However, they did not observe MgCO 3 (C 2/m). In another experimental study of an Mg 0.6 Fe 0.4 CO 3 composition, Boulard et al. (2012) detected the formation of HP-(MgFe)CO 3 (as well as ferropericlase, a high-pressure mag-  (Fiquet et al., 2002), down-pointing triangles (Litasov et al., 2008), squares (Sanchez-Valle et al., 2011) and circles (Yang et al., 2014). netite phase and diamond), above about 75 GPa and 2600 K. HP-(MgFe)CO 3 is the phase reported in one of their earlier studies (Boulard et al., 2011), which is consistent with the C 2/m space group, and fits with our phase diagram. Our results indicate that Fe 4 C 3 O 12 will only be present along a mantle geotherm for Mg 1-x Fe x CO 3 (with x > 0.7), forming in an intermediate layer comprising (1-x) MgCO 3 (C 2/m) + (x/4) (Fe 4 C 3 O 12 + C (diamond)) (Supplementary Material Fig. S10).

Elastic and seismic properties
The elastic properties of carbonate phases can be used to predict the seismic properties of mantle assemblages, which are important for constraining the global carbon cycle and carbon storage. In particular, the seismic anisotropy of Mg 1-x Fe x CO 3 is reported to be significant (Marcondes et al., 2016;Sanchez-Valle et al., 2011;Stekiel et al., 2017;Yao et al., 2018) and could be used as a diagnostic tool in regions of subduction of lithosphere. The elastic constants of MgCO 3 (R3c) have been measured to about 14 GPa (Yang et al., 2014) and calculated to 150 GPa at 0 K (Marcondes et al., 2016;Stekiel et al., 2017) and to 90 GPa at high temperature (Yao et al., 2018). Calculated values have also been reported MgCO 3 (C 2/m) to 150 GPa at 0 K (Marcondes et al., 2016). In contrast, much less is known about the elastic properties of ironbearing carbonate phases. The elastic constants of Mg 0.35 Fe 0.65 CO 3 have been measured to 70 GPa (Fu et al., 2017), while those of the iron end-member FeCO 3 have only been measured at ambient conditions (Sanchez-Valle et al., 2011). Calculated values for FeCO 3 have been reported to 60 GPa, with measurements of c 33 and c 44 (Stekiel et al., 2017). These are insufficient for discussions on their seismic detectability in the mantle.
Our calculated 0 K elastic constants of MgCO 3 (R3c) are plotted in Fig. 7 and listed in Supplementary Material Table S5, Table S6, are consistent with measurements for a Mg 0.35 Fe 0.65 CO 3 composition (Fu et al., 2017) and calculations and measurements of FeCO 3 (Stekiel et al., 2017). c 12 and c 13 soften to such an extent that they become negative above 40 GPa. c 11 and c 33 are also both significantly softened. In contrast, c 44 increase by 65% and c 14 increases by 19%. The elastic constants of MgCO 3 (C 2/m) are listed in Supplementary Material Table S5, with those of Fe 4 C 3 O 12 in Supplementary Material Table S6. Fig. 9 shows a comparison of the V P , V S , and density of MgCO 3 (R3c), MgCO 3 (C 2/m), FeCO 3 (R3c), Fe 4 C 3 O 12 and other major mantle minerals as a function of pressure, at 0 K. At pressures corresponding to the upper mantle and transition zone, the seismic velocities of MgCO 3 (R3c) are very similar to those of forsterite, although its density is much lower. Inclusion of iron reduces its seismic velocities and increases density, meaning that the seismic velocities of FeCO 3 (R3c) are much lower and its density much higher than those of forsterite, wadsleyite and ringwoodite. In the lower mantle, the seismic velocities and density of MgCO 3 (C 2/m) are most similar to those calculated for periclase and bridgmanite. Inclusion of iron reduces velocities and increases density.
Using calculations based on low pressure elastic constants of MgCO 3 (R3c) and estimated volume fractions, previous investigations (Sanchez-Valle et al., 2011;Yang et al., 2014) showed that the differences in velocities of carbonated and non-carbonate eclogite and peridotite are about 1% or less, making the detection of carbonated lithologies from seismic velocities difficult. However, they did not consider a possible high-pressure phase transition. Our calculations predict a phase transition from MgCO 3 (R3c) to MgCO 3 (C 2/m) at 75 GPa. If we make a first-order approximation and neglect temperature, based on our 0 K elastic constants, seismic discontinuities of about +4% in V P and +9% in V S are expected across the phase transition. However, given that the volume fraction of MgCO 3 in carbonated eclogite and peridotite is expected to be less than 15 vol% (Dasgupta et al., 2004; (Kiefer et al., 1999), wadsleyite (Kiefer et al., 2001), MgO (Karki et al., 1997) and MgSiO 3 (Kiefer, 2002). Experimental results for MgCO 3 (R3c) are shown as solid white squares (Sanchez-Valle et al., 2011) and circles (Yang et al., 2014) with a red edge, while those for FeCO 3 (R3c) are shown as solid white squares (Sanchez-Valle et al., 2011) with a black edge and those for Mg 0.35 Fe 0.65 CO 3 (R3c) as solid white circles (Fu et al., 2017) with a black edge. and Hirschmann, 2006;Sanchez-Valle et al., 2011), it is likely that these discontinuities are undetectable.
From our calculated 0 K elastic constants for FeCO 3 (R3c), we predict velocity anomalies across the spin transition of 45% in V P and 25% in V S . Their magnitude will decrease with iron-content and their sharpness decrease with temperature (Hsu and Huang, 2016). At 47 GPa and 0 K, V P of FeCO 3 (R3c) and Fe 4 C 3 O 12 differ by about 6%, while their V S differ by about 14%.
Using our 0 K elastic constants (Supplementary Material Table  S5 and S6), we estimate the anisotropy factor for V P ( A P ) and polarization anisotropy factor for V S ( A S ) for MgCO 3 (R3c), MgCO 3 (C /2m), low-spin FeCO 3 (R3c) and Fe 4 C 3 O 12 (Fig. 10). The centre of each figure is the direction perpendicular to the ab-plane, which is orthogonal to the page. We note that our 0 K elastic constants do not account for temperature. Yang et al. (2014) reported that the influence of temperature on the seismic anisotropy of MgCO 3 is weak, but their study was performed over a limited temperature range and should be confirmed.
We find that at 75 GPa, A S has a maximum value of 52% for MgCO 3 (R3c). This supports previous low pressure studies (Sanchez-Valle et al., 2011;Yang et al., 2014), which report that MgCO 3 (R3c) exhibits greater anisotropy than other major mantle minerals. In contrast, for the high-pressure phase, MgCO 3 (C 2/m), A S has a maximum value of 27%. This suggests that the R3c to C 2/m phase transition leads to a significant decrease in seismic anisotropy. For low-spin FeCO 3 (R3c), A S has a maximum value of 48% at 75 GPa, while that of the high-pressure phase, Fe 4 C 3 O 12 , Fig. 10. Seismic anisotropy of MgCO 3 (R3c), MgCO 3 (C /2m), FeCO 3 (R3c), and Fe 4 C 3 O 12 at 75 GPa. The centre of each figure is the direction perpendicular to the ab-plane, which is orthogonal to the page. Maximum and minimum values are marked by solid black squares and circles. Low-pressure phases (MgCO 3 (R3c) and FeCO 3 (R3c)) exhibit a very high degree of seismic anisotropy, compared to the high-pressure phases (MgCO 3 (C /2m) and Fe 4 C 3 O 12 ). Made using the MATLAB Seismic Anisotropy Toolkit (MSAT) (Walker and Wookey, 2012). is 11%. This implies that the transformation of FeCO 3 (R3c) to Fe 4 C 3 O 12 (diamond), will also result in a decrease in anisotropy.
Previous studies (Sanchez-Valle et al., 2011;Yang et al., 2014) argued that carbonates could develop a high degree of lattice preferred orientation and the large anisotropy of MgCO 3 (R3c) and FeCO 3 (R3c), could be used as a diagnostic feature to detect highly carbonated regions in the upper mantle. Our results support this, showing that, the large anisotropy of MgCO 3 (R3c) and FeCO 3 (R3c), persists at lower mantle pressures. In contrast, our predicted high-pressure phases, MgCO 3 (C 2m) and Fe 4 C 3 O 12 , exhibit much weaker anisotropy, making it likely impossible to detect carbonated minerals in the lowermost mantle.

Outstanding issues
Experimental studies show that carbonates react with SiO 2 to produce silicates + C + O 2 (Drewitt et al., 2019). The conditions at which this reaction takes place is debated. One study indicates that the reaction will not occur along a very cold slab geotherm (Maeda et al., 2017), permitting the subduction of carbonates down to the core-mantle boundary. In contrast, another investigation suggests that reaction should occur above 1500 km, regardless of slab geotherm (Drewitt et al., 2019), meaning no carbonates will persist below this depth. If the former study is correct, Mg 1-x Fe x CO 3 (C 2/m) and Fe 4 C 3 O 12 could be present in the lowermost mantle, with formation of the latter providing an alternative explanation for super-deep diamonds (Wirth et al., 2014). If the latter study is correct, Mg 1-x Fe x CO 3 (C 2/m) is unlikely to present in the lowermost mantle. The reaction between Fe 4 C 3 O 12 with SiO 2 remains to be investigated and it may be stable to greater depths.
Despite many studies of iron-rich carbonates (Boulard et al., 2015(Boulard et al., , 2012(Boulard et al., , 2011Cerantola et al., 2017;Liu et al., 2015;Merlini et al., 2015), none have demonstrated that iron-enrichment will occur, since the starting materials for the experiments did not contain other phases. It has been postulated that the spin transition in Mg 1-x Fe x CO 3 could induce iron enrichment (Cerantola et al., 2015;Liu et al., 2015), in a similar manner to that observed during the spin transition in ferropericlase (Muir and Brodholt, 2016). According to our calculations the degree of iron enrichment needed for Fe 4 C 3 O 12 to form in the lower mantle, is Fe/(Mg+Fe) > 0.7, which is much greater than the maximum of Fe/(Mg+Fe) = 0.35 calculated for Mg 1-x Fe x O (Muir and Brodholt, 2016). Further work is required to investigate iron partitioning to determine if such a degree of iron enrichment is possible, and explore the reason why the Fe 4 C 4 O 13 phase reported by Cerantola et al. (2017) was not observed in other studies (Boulard et al., 2012;Liu et al., 2015).
Our study has only considered the MgCO 3 -FeCO 3 system and neglected the role of CaCO 3 . There have been several studies of the structure and stability of CaCO 3 at lower mantle pressures (Oganov et al., 2008;Pickard and Needs, 2015;Santos et al., 2019;Zhang et al., 2018). The most recent investigation (Zhang et al., 2018) indicates that CaCO 3 readily reacts with SiO 2 and so is unlikely to be a major host of carbon in the lower mantle. However, more work is needed to investigate the MgCO 3 -CaCO 3 -FeCO 3 system.

Conclusions
To summarise, we have performed ab initio calculations to determine the thermodynamic and elastic properties of various phases of MgCO 3 and FeCO 3 , Fe 4 C 3 O 12 and C (diamond), in order to establish the stable phases in the lower mantle and their possible seismic detectability.
Based on our calculations, we predict that Mg 1-x Fe x CO 3 (with x < 0.7) undergoes a transition from R3c to C 2/m structure at conditions corresponding to a depth of about 1800 km in the lower mantle. Inclusion of iron at these concentrations leads to a narrow binary phase loop of only a few gigapascals. For higher iron concentrations, Mg 1-x Fe x CO 3 (with x > 0.7), we predict an additional intermediate layer comprising (1-x) MgCO 3 (C 2/m) + (x/4) (Fe 4 C 3 O 12 + C (diamond)), the thickness of which increases with iron content. In agreement with recent experimental work , we show that FeCO 3 (R3c) undergoes self-oxidation-reduction at conditions corresponding to a depth of about 1300 km in the lower mantle, to form Fe 4 C 3 O 12 + C (diamond). The unusual shape of the phase boundary for the reaction is governed by the oxidation and spin state of iron in the two phases, which leads to different magnetic entropies at high temperature. Similar to a number of previous studies (Marcondes et al., 2016;Sanchez-Valle et al., 2011;Stekiel et al., 2017;Yao et al., 2018), our calculations indicate that the seismic anisotropy of MgCO 3 (R3c) and FeCO 3 (R3c) is extremely high, meaning it might be possible to use it as a diagnostic tool in regions of subduction of lithosphere. In contrast, we find the seismic anisotropy of the corresponding high-pressure phases, MgCO 3 (C 2/m) and Fe 4 C 3 O 12 (R3c), to be considerably weaker, meaning that if carbonates do persist into the lowermost mantle, they are likely not seismically detectable.
Our work highlights the importance of accounting for magnetic entropy in calculations of iron-rich phases in planetary interiors and complex chemistry of iron-bearing minerals, arising from their unpaired d-electrons. Magnetic entropy can also influence phase transitions at lower temperatures (Zhou et al., 2014). Our original investigation began using crystal structure prediction software to predict the high-pressure phase of FeCO 3 . This did not take magnetic entropy into account and could not find a structure with a lower enthalpy than FeCO 3 (R3c). It may be possible to include an expression, similar to Eq. (3), in crystal structure prediction calculations, to act as a first-order approximation for magnetic entropy, when comparing the stability of phases with different magnetic states.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.