Critical vaporization of MgSiO3
See allHide authors and affiliations
Edited by Henry J. Melosh, Purdue University, West Lafayette, IN, and approved March 27, 2018 (received for review November 2, 2017)

Significance
During accretion, giant impacts may partially vaporize planets, with profound consequences for energy and mass transfer between the growing planet and its satellites. We find that giant impacts produce much more supercritical fluid than previously thought, with a much lower liquid–vapor critical temperature. The vapor is depleted in Mg with respect to coexisting liquid and has very different chemical bonding, being dominated by covalent molecules (O2 and SiO) up to the critical point. These molecules also appear in significant fractions within the near-critical liquid phase, making the liquid structure fundamentally different from that near the melting point.
Abstract
Inhomogeneous ab initio molecular dynamics simulations show that vaporization of MgSiO3 is incongruent and that the vapor phase is dominated by SiO and O2 molecules. The vapor is strongly depleted in Mg at low temperature and approaches the composition of the liquid near the critical point. We find that the liquid–vapor critical temperature (
Planets present fundamental challenges to our understanding of materials because of the immense range of pressure and temperature relevant to their evolution and formation. The highest temperatures, exceeding even those at the planetary center, are produced in giant impacts, which appear to be a generic part of planetary accretion. A substantial portion of the planet is vaporized in such events (1, 2); in the case of Earth, the Moon may have condensed from an extended fluid cloud (3, 4). Knowing how much vapor is produced and its composition and speciation is thus essential for understanding the initial stages of planetary evolution. Silicate vaporization is important in other planetary contexts as well, for example, in accreting and steadily evaporating rocky exoplanets (5, 6) and in the passage of meteorites through planetary atmospheres (7).
The conditions of planetary vaporization lie in the realm of warm dense matter, defined by conditions of density and temperature such that interactions among molecules are essential and the electronic structure is not yet fully ionized (e.g.,
Despite the importance of silicate vaporization, our knowledge is still limited because of the very high temperature involved; the liquid–vapor coexistence curve has never been experimentally measured for any silicate composition. The conditions of silicate vaporization extend from the solid–liquid–vapor triple point up to the critical point, which is unknown. In SiO2 previous estimates of the critical point range from 5,000 K to 13,000 K and from 0.07 g⋅
Silicate vaporization presents at least two important challenges to simulation methodology. First, the change in bonding on vaporization means that methods based on classical potentials that are typically tuned to one type of bonding will fail to capture the essential physics. Second is the incongruent nature of silicate vaporization: The chemical composition of the liquid and vapor coexisting in equilibrium may differ (20, 21). This means that methods of determining liquid–vapor coexistence based on the Maxwell equal area construction are incomplete as they assume that liquid and vapor have the same chemical composition. There has been one previous molecular dynamics study of silicate vaporization (SiO2 system) (22), which was based on classical pair potentials and which assumed congruent vaporization based on the Maxwell construction.
Here, we use inhomogeneous ab initio molecular dynamics to simulate silicate vaporization. We have chosen the MgSiO3 system as a prototypical binary silicate that has been extensively studied via first-principles molecular dynamics in the liquid phase (23⇓–25) and that is the dominant chemical component of telluric planets, making up, for example, more than 80% of Earth’s silicate fraction (26).
Our method simulates vaporization directly, mimicking a realizable experimental process (Fig. 1). A slab of liquid is placed in a vacuum. During the dynamical trajectory in the canonical ensemble, we see spontaneous evaporation, leading eventually to a steady-state chemical equilibrium between coexisting liquid and vapor. This approach does not assume congruent vaporization as did the previous molecular dynamics study of silicate vaporization (22). Indeed, we find that the chemical compositions of the vapor and the liquid differ significantly. Our simulations allow us to determine the orthobaric densities of coexisting liquid and vapor, their chemical compositions, the speciation of the vapor phase, and the properties of the interface via direct analysis of the trajectories. Previous two-phase molecular dynamics simulations of liquid–vapor coexistence have been based on classical potentials and have been limited to congruently vaporizing systems, including simple fluids and water (27⇓–29). Further details of our simulations are given in Materials and Methods.
Snapshot from a simulation at 5,500 K (270 atoms,
Fig. 2 shows the results of typical simulations. The density shows the expected spatial variation: from large density within the liquid slab to much lower density in the vapor region. The density of the liquid is reduced at higher temperature, whereas the density of the vapor is augmented. The width of the Gibbs dividing surface grows with increasing temperature. As expected for near critical fluids, density fluctuations in space and time are apparent in both phases and the magnitude of these fluctuations increases with increasing temperature. However, a steady state has been achieved: There are no secular variations in the mean density of either the liquid or the vapor phase.
Simulated mass density profiles at 4,000 K (Top,
The densities of coexisting liquid and vapor approach each other with increasing temperature (Fig. 3) and the phase diagram can be represented by the Wegner expansion (30) which yields the location of the critical point (technically the maxcondentherm) at
Liquid–vapor phase diagram from 270-atom (solid symbols) and 135-atom (open symbols) simulations, with Ω = 4 (circles) and 6 (squares). Red, liquid; blue, vapor. The green solid circle is the critical point of the liquid–vapor boundary assumed in recent simulations of the Moon-forming impact (4, 19). The black solid line is the best fit to the Wegner expansion:
The pressure along the coexistence line increases monotonically with increasing temperature (Fig. 3). We determine the vapor pressure as the component of the pressure tensor normal to the interface (
The width of the Gibbs dividing surface between liquid and vapor phases increases systematically with increasing temperature (Fig. 4). We have found that our results can be fitted to the expression (33)
The width of the interface (blue) and the compressibility (red, right-hand axis) as a function of temperature.
As temperature increases, spatial fluctuations in the density of the liquid phase increase in magnitude. We quantify density fluctuations and show that their magnitude is thermodynamically sensible (Fig. 4). Density fluctuations of the liquid may be related to the isothermal compressibility and to the structure factor in the limit of small q as
The compressibility shows a steep increase with increasing temperature from 0.1 GP
The structure of the liquid undergoes important changes with increasing temperature along the vapor–liquid coexistence line (Fig. 5). At 4,000 K, the Si-O, Mg-O, and O-O radial distribution functions are very similar to those seen in previous studies of the homogeneous liquid at similar conditions of temperature and density (23, 36). As temperature increases a new peak emerges in the O-O radial distribution function at a very small distance (1.3 Å). This peak corresponds to the appearance of O2 molecules in the liquid phase.
Cation–oxygen coordination numbers (Top) and the abundance of molecular species in the liquid (Bottom) shown as the proportion of O atoms that are members of O2 molecules and the proportion of Si atoms that are members of an SiO molecule. Bottom Inset shows the partial radial distribution functions of the liquid phase at 6,000 K: Si-O (blue), O-O (red), and Mg-O (gold).
We quantify liquid structure via the coordination number
We find that the Si-O (4) and Mg-O (5) coordination numbers at 4,000 K are very similar to those found in previous experimental and theoretical studies near the triple-point density (23, 36) (Fig. 5). With increasing temperature,
The vapor is significantly depleted in Mg compared with the liquid (Fig. 6): At 5,000 K, only 7% of the vapor is Mg, compared with 20% for bulk MgSiO3 composition. The vapor is enriched in O compared with the bulk composition. With increasing temperature, the O content in the vapor decreases, and the Mg content increases, so that the vapor composition approaches the bulk composition at the critical point.
Chemical composition of the vapor phase: atomic fractions of O (red), Si (blue), and Mg (gold). Also shown are the predictions of the MAGMA code (dashed lines). The arrows along the right-hand axis indicate the atomic fractions in the bulk (MgSiO3) composition.
The vapor phase in our simulations consists of several clearly defined and long-lived molecular species, the most important of these being SiO and O2. Because molecules interact frequently, particularly near the critical point, through either collisions or exchanges of atoms, molecules cannot be uniquely defined. Definitions based on bond-length criteria are unsuitable as vibrations are large in amplitude and highly anharmonic so that molecular proportions are sensitive to the choice of bond-length cutoffs. Instead, we use a definition of a molecule that does not rely on arbitrary bond-length distance cutoffs. We define a two-atom cluster as a molecule when the atoms are mutual nearest neighbors. We also permit three-atom molecules if (i) an Si atom i and an O atom j are mutual nearest neighbors and (ii) another O atom k has i as its nearest neighbor and k is the second nearest neighbor of i. The concept of mutual nearest neighbors has been used previously to identify dimers in dense fluid hydrogen (37) and generalized to a consideration of Nth mutual neighbors in the context of simple liquids (38).
We find that the vapor is dominated by SiO (33% at 4,500 K) and O2 (33%) (Fig. 7). The next most abundant species are O (12%), Mg (10%), and SiO2 (9%). With increasing temperature SiO and O2 remain the most abundant species, although their relative proportions decrease somewhat, while the fractions of atomic O and Mg increase. The proportion of SiO2 molecules depends nonmonotonically on temperature, first increasing with temperature, up to 5,500 K, and then decreasing on further heating.
Speciation in the vapor phase. In Top, Middle, and Bottom, monatomic species are indicated by open symbols and molecular species by solid symbols. Top Inset shows the lifetimes of the molecular species compared with their vibrational periods indicated by the tags on the right-hand axis and the vertical blue band joining the periods of the stretching and bending mode of SiO2 (39, 40).
Our results highlight the limitations of the concept of a molecule: At the temperatures under consideration molecular lifetimes may be less than one vibrational period (Fig. 7). While O2 and SiO molecules are long-lived, the lifetime of SiO2 is much less than the period of the bending mode and the lifetime of MgO is much less than the period of the stretching mode. The dynamics and geometry of the SiO2 and MgO “molecules,” and thus their energetics, may differ significantly from what would be expected on the basis of extrapolation from low-temperature experimental data.
As the liquid–vapor coexistence curve has never been determined before via experiment or first-principles simulation, there are few data with which to compare our results. The one previous simulation study of silicate vaporization (on SiO2) (22) used classical potentials and a homogeneous simulation approach, thereby assuming congruent vaporization. The assumption of congruency and the inability of classical potentials to capture the stabilization of covalent species in the vapor phase, in addition to the difference in system composition, likely account for the much higher critical temperature found in the previous study (12,000 K vs.
Experimental data on vapor phase composition in the MgSiO3 system do not exist. However, our lowest-temperature results are in reasonable accord with available experimental data in the vicinity of the triple point. The silica component appears to be more volatile than the magnesia component in MgSiO3 (21) as well as in a wide variety of multicomponent silicates (41), so that evaporation tends to produce vapor that is Mg depleted, as we find.
We may also compare our results with extrapolations of thermochemical data. For this purpose, we have used the MAGMA code (16), computing the composition of the vapor phase in equilibrium with MgSiO3 liquid. The MAGMA code also predicts vapor depleted in Mg, as we find, and vapor speciation dominated by SiO and O2 with O and SiO2 being next most abundant. We find the composition of the vapor changes more rapidly with increasing temperature than predicted by the MAGMA code, which shows more nearly constant vapor-phase composition (Fig. 6). These differences are due in part to the ideal gas law underlying the MAGMA code, which breaks down as the critical point is approached; our results could be used to expand the scope of such modeling efforts. The MAGMA code also makes a series of assumptions regarding the component species in the liquid phase and does not include SiO or O2 as liquid-phase species. The results of our study may be used to improve the liquid-phase model, for example, by providing information on these molecular components. A recent study of the silica system (15) argued that it is important to include SiO and O2 as liquid-phase species, particularly near the critical point.
The temperature of the liquid–vapor critical point that we find (6,600 K) is much less than what has been assumed in hydrodynamic simulations of planetary accretion (8,800 K) (Fig. 3). The much higher critical temperature assumed in hydrodynamic simulations reflects the properties of the underlying semianalytic thermodynamic model, which is based on the MANEOS approach (42), and the lack of experimental data on near-critical vaporization. Our results point toward the production of much greater amounts of supercritical fluid over a wider range of impact scenarios than previously assumed. The proto-Earth and its surrounding postimpact disk may have been supercritical out to greater distances, possibly encompassing the radius at which the Moon formed. The presence of such an extended one-phase structure may have important implications for understanding the similarity between the Earth and the Moon in O and other isotopes (43) and for isotopic fractionations between the two bodies such as those recently observed in potassium (3).
The incongruent nature of vaporization that we find highlights another shortcoming of the phase diagram used in hydrodynamic simulations: Silicate vaporization is assumed to be congruent. This assumption is made for simplicity, but also based on the experimental observation that forsterite Mg2SiO4 vaporizes congruently near the triple point. However, the Mg/Si ratio of forsterite is very different from that of the bulk silicate Earth, which is much closer to that of our system (Mg/Si = 1.3). Moreover, thermodynamic modeling of the Mg2SiO4–SiO2 system, based on experimental observations, suggests that Earth-like compositions vaporize incongruently to Mg-poor vapor at low temperature (21) as we find. We find that the degree of Mg depletion in the vapor depends strongly on temperature and that vapor and liquid have very similar compositions near the critical point.
Our results open avenues for studying silicates in the lower-density portion of the warm dense-matter regime. Our methods are in principle applicable to any silicate composition, including those that have been studied more extensively by experiment, such as SiO2, and multicomponent systems that capture even more of the richness of planetary vaporization. By the same token, current experimental technology is in principle capable of directly testing our current predictions: Shock and release studies have already measured fluid structure via XANES (X-ray absorption near edge structure) down to 0.9 g⋅
Materials and Methods
Our molecular dynamics simulations are based on density functional theory using the projector augmented wave (PAW) method (44) as implemented in VASP (45), with valence configurations Si (3s23p2), Mg (3s23p0), and O (2s22p4). We use the PBEsol approximation to the exchange-correlation potential (46) as we have previously found it to yield excellent agreement with experimental measurements of condensed-phase physical properties of oxides (47, 48). We sample the Brillouin zone at the Gamma point and use a plane-wave cutoff of 500 eV, which we find yields pressures and energies that are converged to within 0.01 GPa and 2 meV per atom, respectively. We assume thermal equilibrium between ions and electrons via the Mermin functional (49).
The initial condition is a liquid slab embedded in a vacuum. The simulations are performed in the canonical ensemble using the Nosé–Hoover thermostat (50) with a time step of 1 fs for durations of 4–15 ps. We explored a range of system sizes: 135 atoms or 270 atoms and initial vacuum-to-liquid volume ratios Ω = 4 or 6. Quantities are computed as time averages over the stationary portion of the trajectory, taken to the be the last 50%, and uncertainties are computed with the proper non-Gaussian statistics via the blocking method (51). Initial configurations are chosen from equilibrated homogeneous simulations near zero pressure from our previous work (25). This choice of initial condition helps to minimize behavior that we have observed which complicates analysis: a tendency for the liquid slab to split in two, particularly as the critical point is approached. The homogeneous simulations were performed with N = 135. To initiate our N = 270 simulations, we doubled the liquid slab in the direction normal to the interface. Because of the large cell sizes, these are very demanding simulations: The number of plane waves required is as much as 100 times greater than for homogeneous liquid-phase simulations.
To determine the physical properties and composition of the two coexisting phases, we locate the liquid–vapor interface using the method of ref. 52. Briefly, the interface is defined at every time step by the surface s such that
We compute the mean mass density profile
Computation of the
Acknowledgments
We thank Prof. Bruce Fegley (Washington University) for kindly providing the MAGMA code, Zhigang Zhang (Chinese Academy of Sciences) for helpful discussions, and K. Kurosawa and J. Connolly for insightful comments on the manuscript. This project is supported by the European Research Council under Advanced Grant 291432 MoltenEarth (FP7). Calculations were carried out using the IRIDIS cluster partly owned by University College London, ARCHER of the United Kingdom national high-performance computing service, and MARENOSTRUM at the Barcelona Supercomputing Center, Spain of the PRACE Consortium.
Footnotes
- ↵1To whom correspondence should be addressed. Email: l.stixrude{at}ucl.ac.uk.
Author contributions: B.X. and L.S. designed research, performed research, analyzed data, and wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
- Copyright © 2018 the Author(s). Published by PNAS.
This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).
References
- ↵
- Pahlevan K,
- Stevenson DJ
- ↵
- Nakajima M,
- Stevenson DJ
- ↵
- Wang K,
- Jacobsen SB
- ↵
- Lock SJ,
- Stewart ST
- ↵
- Visscher C,
- Fegley B
- ↵
- van Lieshout R, et al.
- ↵
- ↵
- Lee RW, et al.
- ↵
- Clerouin J, et al.
- ↵
- ↵
- Kurosawa K, et al.
- ↵
- Kraus RG, et al.
- ↵
- Denoeud A, et al.
- ↵
- ↵
- Connolly JAD
- ↵
- Fegley B,
- Cameron AGW
- ↵
- Roine A
- ↵
- Schaefer L,
- Fegley B
- ↵
- ↵
- Nakagawa H,
- Asano M,
- Kubo K
- ↵
- ↵
- ↵
- Stixrude L,
- Karki B
- ↵
- ↵
- Karki BB,
- Zhang J,
- Stixrude L
- ↵
- Workman RK,
- Hart SR
- ↵
- Rao M,
- Levesque D
- ↵
- Alejandre J,
- Tildesley DJ,
- Chapela GA
- ↵
- ↵
- ↵
- ↵
- Ahrens TJ,
- O’Keefe JD
- ↵
- Beysens D,
- Robert M
- ↵
- Hansen JP
- ↵
- Rivers ML,
- Carmichael ISE
- ↵
- ↵
- ↵
- Larsen RE,
- Stratt RM
- ↵
- Irikura KK
- ↵
- Pacansky J,
- Hermann K
- ↵
- ↵
- ↵
- Young ED, et al.
- ↵
- ↵
- Kresse G,
- Furthmuller J
- ↵
- ↵
- ↵
- Scipioni R,
- Stixrude L,
- Desjarlais MP
- ↵
- ↵
- ↵
- ↵
- ↵
Citation Manager Formats
Article Classifications
- Physical Sciences
- Earth, Atmospheric, and Planetary Sciences