Critical vaporization of MgSiO3.

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 ([Formula: see text] K) is much lower than assumed in hydrodynamic simulations, pointing to much more extensive supercritical fluid after the Moon-forming impact than previously thought. The structure of the near-critical liquid is very different from what has been studied previously and includes a significant fraction (10%) of molecular species SiO and O2.

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., – g 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, O2 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 SiO2 previous estimates of the critical point range from 5,000 K to 13,000 K and from 0.07 g to 1.1 g, 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 (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 .

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 O2, SiO, SiO2, 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.

Fig. 2.

Simulated mass density profiles at 4,000 K (Top, ) and 5,500 K (Bottom, ). 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. to the time-averaged density profile. Both simulations have = 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 g and K. The density of the liquid decreases monotonically with increasing temperature. At 5,500 K the density ( g) is half that at the ambient melting point of this system (31) (2.58 g at 1,830 K). The density of the coexisting vapor at 5,500 K is 0.14 0.05 g. The density of the vapor phase increases with increasing temperature. We find no systematic variation of liquid or vapor densities with computational parameters or , as expected for orthobaric two-phase coexistence (28).

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: , ; with values of the exponents fixed to those expected from renormalization group theory, , (30); and the other parameters , 2.36, 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 () (29). The pressure may be represented by with and with in megapascals. This fit yields a critical pressure MPa at 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 = 12.45 and = 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) which yields Å, K, and . The value of determined from the fit agrees well with the critical temperature determined from the phase diagram (Fig. 3).

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 aswhere is the particle number, angle brackets denote time averages, is the structure factor, is the isothermal compressibility, is the number density, is temperature, and is the Boltzmann constant. We compute the structure factor from the definition (34) with given by a sum over the positions of the atoms in the liquid phase, i.e., those for which Eq. is satisfied. To focus on the properties of the liquid phase, we select vectors from the plane parallel to the interfacewith and integers and the cell dimension in the plane of the interface. To find we fit the structure factor to a quadratic in in the range . 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 value, .

The compressibility shows a steep increase with increasing temperature from 0.1 GP at 4,000 K to 0.5 GP 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 GP) (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 O2 molecules in the liquid phase.

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 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 numberwhere is the mean number of atoms of type surrounding an atom of type , is the corresponding partial radial distribution function, is the concentration of atoms of type , and is the location of the first minimum in . We also track , the fraction of atoms surrounded by atoms. All liquid-phase structural quantities are computed from positions of atoms that satisfy Eq. .

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, and decrease to values near 3 close to the critical point. The proportion of O atoms that are members of an O2 molecule () increases steadily with increasing temperature, as does the number of SiO molecules (), 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 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.

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 (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 and an O atom are mutual nearest neighbors and (ii) another O atom has as its nearest neighbor and is the second nearest neighbor of . The concept of mutual nearest neighbors has been used previously to identify dimers in dense fluid hydrogen (37) and generalized to a consideration of th 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.

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 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. K in our study).

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 in SiO2 (13), vapor composition in CaMgSi2O6 (11), and the thermodynamic conditions of one point on the liquid–vapor coexistence line in SiO2 (12).

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 = 135. To initiate our = 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 where is a constant, and is the coarse-grained density field obtained by convolving atomic positions with a Gaussian. The proximity of atom to the surface iswhere is the position of atom , is the point on the surface nearest to atom , and is the surface normal in the direction of the density gradient at that point. We have found that this method yields robust location of the dividing interface even in those cases where the surface is time dependent or when multiple liquid–vapor surfaces appear as a result of splitting of the liquid slab. We also locate the Gibbs dividing surface and determine its properties by replacing the instantaneous interface with the time-averaged surface in the computation of (52).

We compute the mean mass density profilewhere is the mass of atom and is the dimension of the cell parallel to the interface. We find that the density profile follows the expected form (53)from which we extract our best-fitting values of and , the density of the liquid and vapor phases, respectively, and , the width of the dividing surface.

Computation of the allows us to identify the phase to which each atom belongs at each time step in the simulation. We assume that atoms belonging to the liquid and vapor phases are defined byandrespectively. This identification allows us to determine several properties of the liquid and vapor phases, including composition, speciation, and structure.