• ab initio molecular dynamics simulation
• MgSiO₃
• vaporization
• giant impact
• density functional theory

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., 10−3–101 g⋅cm−3 and 0.5–20 eV) (8). Experimental interest has expanded considerably with the development of methods to produce and accurately characterize matter in this regime, for example on adiabatic release from either shocked or isochorically heated states and because of applications to understanding astrophysical objects, plasma production and characterization, and inertial confinement fusion (8, 9). In silicates most dynamic compression studies have focused on achieving the high pressures characteristic of planetary interiors (10), whereas less attention has focused on expanded states and vaporization (11⇓–13). The behavior of silicates under expansion may be particularly rich as experiments near the triple point suggest that a major change in bonding occurs: from dominantly ionically bonded solid and liquid phases to dominantly covalently bonded vapor, containing, for example, O₂ and SiO molecules (14).

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 SiO₂ previous estimates of the critical point range from 5,000 K to 13,000 K and from 0.07 g⋅cm−3 to 1.1 g⋅cm−3, illustrating the magnitude of uncertainty that remains (15). Predictions based on the extrapolation of experimental data (16⇓–18) are accurate near the triple point, but their range of validity at higher temperatures is unknown. Silicate vaporization in hydrodynamic simulations of giant impact events (4, 19) is based on a semianalytic model, which assumes that the liquid–vapor critical point occurs at 8,800 K and that vaporization is congruent, but the accuracy of these assumptions is untested experimentally near the critical point.

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 (SiO₂ 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 MgSiO₃ 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Snapshot from a simulation at 5,500 K (270 atoms, Ω= 6) showing Mg (gold), Si (blue), and O (red) atoms. The box outlines the periodically repeated simulation cell. Molecular species are apparent in the vapor phase, including O₂, SiO, SiO₂, and atomic species, Mg and O.

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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Simulated mass density profiles at 4,000 K (Top, Ω=4) and 5,500 K (Bottom, Ω=6). We represent temporal and spatial variations in the density with density profiles from snapshots taken from the simulation at times arranged by rainbow color from 2 ps (red) to 5 ps (blue). The density is coarse grained with a Gaussian of width 2.5 Å and the black line is the best fit of Eq. 6 to the time-averaged density profile. Both simulations have N = 270 atoms.

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 ρc=0.48±0.05 g⋅cm−3 and Tc=6,600±150 K. The density of the liquid decreases monotonically with increasing temperature. At 5,500 K the density (1.29±0.11 g⋅cm−3) is half that at the ambient melting point of this system (31) (2.58 g⋅cm−3 at 1,830 K). The density of the coexisting vapor at 5,500 K is 0.14 ± 0.05 g⋅cm−3. The density of the vapor phase increases with increasing temperature. We find no systematic variation of liquid or vapor densities with computational parameters N or Ω, as expected for orthobaric two-phase coexistence (28).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

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: ρl−ρv=C11−T/Tcβ+C21−T/Tcβ+Δ, 12(ρl+ρv)=ρc+C3(1−T/Tc); with values of the exponents fixed to those expected from renormalization group theory, β=0.325, Δ=12 (30); and the other parameters C1=1.32, C2= 2.36, C3= 1.74. Inset shows the pressure of the two-phase system from our simulations (red symbols), the best fit to the simulation results (red curve), the critical pressure (blue symbol), and the extrapolation of low pressure–temperature experimental data described in the text (black dashed curve).

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 (Pzz) (29). The pressure may be represented by P=exp(A−B/T) with A=11.8±2.0 and B=45,000±14,000 with P in megapascals. This fit yields a critical pressure Pc=140 MPa at Tc=6,600 K. Our pressure–temperature curve agrees remarkably well with an extrapolation from experimental observations at lower temperature: Fitting to the triple point from ref. 21 (0.2 Pa, 1,773 K) and the 1-bar boiling point from ref. 32 (0.1 MPa, 3,350 K) yields A = 12.45 and B = 49,420 and a curve that differs from the best fit to the simulation results by no more than a few megapascals. Incongruent vaporization requires that the vapor pressure line is actually a two-phase coexistence region of finite width. However, a thermochemical modeling study of the silica system (15) indicates that the width of the two-phase region is on the order of 10% of the critical pressure and therefore not resolvable in our simulations within our present uncertainties.

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) w=w01−T/Tc−ν which yields w0=1.4 ±0.4 Å, Tc=6,700±300 K, and ν=0.3±0.2. The value of Tc determined from the fit agrees well with the critical temperature determined from the phase diagram (Fig. 3).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

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⟨N2⟩−⟨N⟩2⟨N⟩=ρNkBTχT=limq→0S(q),[1]where N is the particle number, angle brackets denote time averages, S(q) is the structure factor, χT is the isothermal compressibility, ρN is the number density, T is temperature, and kB is the Boltzmann constant. We compute the structure factor from the definition (34) S(q)=⟨ρq→ρ−q→⟩/N with ρq→=∑i⁡exp(iq→⋅r→i) given by a sum over the positions r→i of the atoms in the liquid phase, i.e., those for which Eq. 7 is satisfied. To focus on the properties of the liquid phase, we select q→ vectors from the plane parallel to the interfaceq→=2πL(h,k,0)[2]with h and k integers and L the cell dimension in the plane of the interface. To find S(q→0) we fit the structure factor to a quadratic in q in the range 0<q<3Å−1. We have found that this procedure is robust for the relatively large compressibilities of our system: The extrapolated value of the structure factor does not differ significantly from that at the smallest q value, S(2π/L).

The compressibility shows a steep increase with increasing temperature from 0.1 GPa−1 at 4,000 K to 0.5 GPa−1 at 6,500 K (Fig. 4). The compressilbity at 4,000 K is two times larger than that measured experimentally near the triple point (0.05 GPa−1) (35).

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 O₂ molecules in the liquid phase.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

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 O₂ 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 numberZαβ=4πρNxβ∫0rminr2gαβ(r)dr,[3]where Zαβ is the mean number of atoms of type β surrounding an atom of type α, gαβ is the corresponding partial radial distribution function, xβ is the concentration of atoms of type β, and rmin is the location of the first minimum in gαβ. We also track ζαβ(Z), the fraction of α atoms surrounded by Z β atoms. All liquid-phase structural quantities are computed from positions of atoms that satisfy Eq. 7.

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, ZSiO and ZMgO decrease to values near 3 close to the critical point. The proportion of O atoms that are members of an O₂ molecule (ζOO(1)) increases steadily with increasing temperature, as does the number of SiO molecules (ζSiO(1)), both approaching 10% near the critical point (Fig. 5).

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 MgSiO₃ 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.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

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 (MgSiO₃) 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 O₂. 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 O₂ (33%) (Fig. 7). The next most abundant species are O (12%), Mg (10%), and SiO₂ (9%). With increasing temperature SiO and O₂ remain the most abundant species, although their relative proportions decrease somewhat, while the fractions of atomic O and Mg increase. The proportion of SiO₂ molecules depends nonmonotonically on temperature, first increasing with temperature, up to 5,500 K, and then decreasing on further heating.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 7.

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 SiO₂ (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 O₂ and SiO molecules are long-lived, the lifetime of SiO₂ 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 SiO₂ 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 SiO₂) (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. Tc=6,600 K in our study).

Experimental data on vapor phase composition in the MgSiO₃ 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 MgSiO₃ (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 MgSiO₃ liquid. The MAGMA code also predicts vapor depleted in Mg, as we find, and vapor speciation dominated by SiO and O₂ with O and SiO₂ 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 O₂ 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 O₂ 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 Mg₂SiO₄ 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 Mg₂SiO₄–SiO₂ 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 SiO₂, 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⋅cm−3 in SiO₂ (13), vapor composition in CaMgSi₂O₆ (11), and the thermodynamic conditions of one point on the liquid–vapor coexistence line in SiO₂ (12).

• Copyright © 2018 the Author(s). Published by PNAS.