## Abstract

First-principles molecular-dynamics simulations show that over the pressure regime of Earth's mantle the mean silicon-oxygen coordination number of magnesium metasilicate liquid changes nearly linearly from 4 to 6. The density contrast between liquid and crystal decreases by a factor of nearly 5 over the mantle pressure regime and is 4% at the core-mantle boundary. The ab initio melting curve, obtained by integration of the Clausius-Clapeyron equation, yields a melting temperature at the core-mantle boundary of 5400 ± 600 kelvins.

The structure of silicate liquids over most of the mantle pressure regime is unknown. There is experimental evidence that silicate liquid structure changes substantially on compression and that these changes include an increase in Si-O coordination number from 4 toward 6 (*1*), although there are no experimental observations within the liquid stability field at lower mantle pressure.

Most simulations of silicate liquids have been based on atomistic models, which permit much faster computation but have the disadvantage of being based on empirical force fields, the forms of which are uncertain (*2*, *3*). First-principles molecular-dynamics (FPMD) simulations are considerably more costly because the electronic structure must be computed at each time step, but they are also more robust because they make no assumptions about the type of bonding or the shape of the charge density.

Here, we present FPMD simulations of MgSiO_{3} liquid over the pressure-temperature range of Earth's mantle, based on density functional theory in the local density approximation and the plane-wave pseudopotential method (*4*). The simulations are performed in the canonical ensemble with periodic boundary conditions and a Nosé thermostat (*5*). The primary cell is cubic, with 80 atoms (16 formula units). The initial condition is a pyroxene structure homogeneously strained to a cubic cell shape and the desired volume. The simulated strained crystal is melted at 6000 K and then cooled isochorically to 4000 K and 3000 K. Melting was confirmed by inspection of the radial distribution function and the mean square displacement. Total run durations are 3 ps, with the last 2.4 ps used for computing equilibrium averages. To ensure that the results are not sensitive to initial conditions, we performed simulations with twice the number of atoms (160 atoms initiated with a cubic majorite structure), twice the duration (6 ps), and different initial configurations (80-atom perovskite super-cell homogeneously strained to a cubic cell shape) and found that the computed properties of the liquid were unchanged within statistical uncertainties. We explore a range of volumes *V* = *V*_{X} to *V* = *V*_{X}/2, where *V*_{X} = 38.9 cm^{3} mol^{–1} is the experimental value for the liquid at the ambient-pressure melting point (1830 K) (*6*).

To explore freezing, we also performed simulations of orthorhombic MgSiO_{3} perovskite, following the same procedures as the liquid simulations except that at each volume the shape of the unit cell was taken to be that found by structural relaxation at static conditions. We found that the unstrained perovskite structure persisted (meta) stably at all volumes *V* < 0.7*V*_{X}. We are thus able to compare directly the simulated properties of liquid and perovskite over the entire lower mantle pressure-temperature regime.

Inspection of the equilibrated liquid structure shows that compression has a large influence on the atomic arrangement (Fig. 1). At the experimental zero-pressure volume, the liquid shows four-fold Si-O coordination, as does the stable crystalline phase at this volume, but without the elements of longer ranged order exhibited by the crystal, such as the infinitely long linear chains of tetrahedra that characterize pyroxenes. The liquid also shows a few non-four-fold-coordinated silicon atoms. At high pressure, the liquid structure is much more densely packed.

In the simulations, the average Si-O coordination number increased smoothly and nearly linearly with compression, with no discernible influence of temperature, reaching a value of 6 at a pressure corresponding to the base of Earth's mantle (Fig. 1). The coordination change does not take place over a narrow pressure interval, as has been suggested (*7*). Such a rapid structural change would be expected to cause changes in the compressibility of the liquid (*7*) and in the slope of the melting curve (*8*) over a correspondingly narrow pressure interval, anomalies for which we find no evidence in our simulations.

The liquid equation of state can be accurately described by the Mie-Gr”neisen form (1) where *P*_{c} is the reference isotherm, taken to be at *T*_{0} = 3000 K, γ is the Gr”neisen parameter, and *C*_{V} is the isochoric heat capacity. The latter quantities are computed from the simulations by noting that the internal energy *E* and pressure *P* vary linearly with temperature along isochores to within our uncertainty. The value of γ is then calculated at each volume from (2) and *C*_{V} = (*∂*E/∂*T*)_{V}. We find that *C*_{V} decreases by about 10% from 4.12 ± 0.11 *Nk* at *V* = *V*_{X} [the experimental value is 4.17 ± 0.18 *Nk* (*9*)] to 3.54 ± 0.11 *Nk* at *V* = *V*_{X} /2, where *N* is the number of atoms and *k* is the Boltzmann constant. γ increases by a factor of three over the same range of volume, reflecting the increase in thermal pressure. All known mantle crystalline phases show the opposite behavior: γ decreases on compression (*10*).

The unusual behavior of γ in the liquid can be understood on the basis of the pressure-induced change in liquid structure (Fig. 2). Experimental data on crystalline MgSiO_{3} phases show that, although compression reduces the value of γ in each individual phase, polymorphic phase transformations have an even larger and opposing effect, with the higher coordinated phase showing the larger value of γ. As the coordination number of the liquid increases on compression, the liquid adopts values of γ characteristic of the more highly coordinated state. The behavior of the liquid is consistent with theoretical analyses of the influence of pressure-induced coordination changes on γ (*11*). Previous theoretical and experimental studies have also found that γ increases on isothermal compression in liquids (*12*–*14*).

The volume of the liquid closely approaches, but does not fall below, that of perovskite over the mantle pressure regime (Fig. 3). The volume difference is 4% at 6000 K and 140 GPa, as compared with 18% at the ambient-pressure melting point (*6*). Studies based on dynamic compression also find that the volume of fusion is markedly reduced at deep lower mantle conditions (*15*). The enthalpy difference between liquid and perovskite increases nearly linearly with increasing temperature and changes little with compression over the pressure regime of perovskite stability (25 < *P* < 140 GPa) (Fig. 3). The quantity Δ*H*/*T* varies little throughout this regime and has a value near 1.5 *Nk*, significantly greater than the value *Nk* ln2 often assumed for the entropy of fusion in high-pressure melting studies, and based primarily on one-component, close-packed materials. The larger value in our simulations can be understood on the basis of liquid structure: Throughout the lower mantle pressure regime, the liquid shows a wide range of different coordination environments, which adds a configurational component to the entropy that is absent in structurally simpler liquids (*16*). A recent FPMD study of MgO also found values of the entropy of melting ∼1.5 *Nk* at high pressure (*17*).

We used our simulations to predict the melting curve of MgSiO_{3} by integrating the Clausius-Clapeyron equation (*18*) (3) The quantities on the right-hand side are computed from our simulations over the entire lower mantle regime. The integration constant is set by assuming a single fixed point along the melting curve. We choose this point as 25 GPa, 2900 K near the low-pressure limit of perovskite melting, where many different experiments are mutually consistent.

Our predicted melting curve lies between previous experimental determinations, which diverge widely at pressures above 30 GPa (Fig. 4). The differences between theory and experiment can be understood in detail. Because we find Δ*V* > 0, our melting slope and melting temperature must be greater than experimental results with a vanishingly small melting slope, which requires Δ*V* = 0 (*19*, *20*). Zerr and Boehler (*21*) argued that their much higher melting slope is consistent with the Lindemann law, a simple scaling argument that is valid if neither solid nor liquid undergo considerable structural change on compression (*22*). This picture of melting is incompatible with the large structural changes that we find in MgSiO_{3} liquid. Moreover, because the structural changes we see lead to more efficient packing, and therefore to reduced liquid-solid volume contrast, our melting slope and melting temperature should be correspondingly reduced. Some of the discrepancy among experiments may be due to differences in composition—Fe-free (*21*) versus ∼10% Fe (*8*, *19*, *20*)—but the wide divergence among diamond-anvil cell studies more likely reflects difficulties in detecting melt at high pressure (*23*). In this context, we note that a recent determination based on dynamic compression (*15*) agrees well with our predicted melting curve.

The marked decrease in the volume of fusion of MgSiO_{3} that we find means that liquids are likely to be denser than coexisting solids in the deep mantle. In a multicomponent system, the density contrast between liquid and coexisting solid is a function of the volume contrast and the difference in chemical composition. In the shallow mantle, the volume contrast dominates, and liquids are less dense despite being enriched in heavy major elements such as Fe and Ca. In the deep mantle, where the volume contrast is small, the difference in composition is likely to be the determining factor. Experiments show that Fe and Ca continue to be preferentially incorporated in the liquid phase at least up to 25 GPa (*24*). For example, an apparent Fe-Mg partition coefficient between solid and liquid *K* = [(*X*_{Fe}^{sol} *X*_{Mg}^{liq})/(*X*_{Fe}^{liq} *X*_{Mg}^{sol}) ≈ 0.4, where *X*_{A}^{α} is the mole fraction of element *A* in phase or phase assemblage α, would be sufficient to overcome the volume difference between MgSiO_{3} liquid and perovskite and make the liquid denser.

The mantle solidus has been estimated to lie 800 to 1300 K below that of pure MgSiO_{3} (*25*). Combining this estimate with our ab initio value of the melting point at 136 GPa of 5400 ± 600 K yields a mantle solidus temperature of 3500 to 5200 K, which overlaps a range of previous estimates of the temperature at the core-mantle boundary (*26*, *27*). Melting of normal mantle at its base cannot be ruled out and appears to be plausible.

The presence of neutrally or negatively buoyant melt in the deep mantle would have fundamental implications for our understanding of the Hadean Earth (*28*) and of the structure and chemical evolution of the coremantle boundary. If melting extended to the base of the mantle during accretion, perovskite may have floated in coexisting liquid as it crystallized out of the magma ocean, providing a possible means of large-scale chemical differentiation. In the present-day Earth, dense melt may provide a natural explanation for the ultralow velocity zone at the base of the mantle (*29*).