# Critical point and supercritical regime of MgO

T. F. J. Bögels, Razvan Caracas

► To cite this version:

T. F. J. Bögels, Razvan Caracas. Critical point and supercritical regime of MgO. Physical Review B, 2022, 105, 10.1103/PhysRevB.105.064105. insu-03643043

HAL Id: insu-03643043

https://insu.hal.science/insu-03643043

Submitted on 6 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Critical point and supercritical regime of MgO

T. F. J. Bögels<sup>1,*,</sup> and R. Caracas<sup>1,2,3,†</sup>

<sup>1</sup>Centre National de la Recherche Scientifique, Laboratoire de Géologie de Lyon, LGLTPE UMR 5276, Centre Blaise Pascal, Ecole Normale Supérieure de Lyon, F-69364 Lyon Cedex 07, France

<sup>2</sup>The Center for Earth Evolution and Dynamics (CEED), University of Oslo, 0371 Oslo, Norway

<sup>3</sup>Université de Paris, Institut de Physique du Globe de Paris, CNRS, Paris, France

<figure><img src="image_2.png" /></figure>

(Received 31 July 2021; revised 9 December 2021; accepted 13 January 2022; published 18 February 2022)

The position of the critical point determines the top of the liquid-vapor coexistence dome, and it is a physical parameter of fundamental importance in the study of high-energy shocks, including those associated with large planetary impacts. For most major planetary materials, such as oxides and silicates, the estimated position of the critical point is below 1 g/cm³ at temperatures above 5000 K. Here we compute the position of the critical point of one of the most ubiquitous materials: MgO. For this we perform first-principles molecular dynamics simulations. We find the critical density to be in the 0.45–0.6 g/cm³ range and the critical temperature in the 6500–7000 K range. We investigate in detail the behavior of MgO in the subcritical and supercritical regimes, and we provide insight into the structure and chemical speciation. We see a change in Mg-O speciation toward lower degrees of coordination as the temperature is increased from 4000 to 10 000 K. This change in speciation is less pronounced at higher densities. We observe the liquid-gas separation in nucleating nanobubbles at densities below the liquid spinodal. The majority of the chemical species forming the incipient gas phase consists of isolated Mg and O atoms and some MgO and O₂ molecules. We find that the ionization state of the atoms in the liquid phase is close to the nominal charge, but it almost vanishes close to the liquid-gas boundary and in the gas phase, which is consequently largely atomic.

DOI: 10.1103/PhysRevB.105.064105

## I. INTRODUCTION

Giant impacts are a characteristic of the early stages of the evolution of planets, when chaotic trajectories of planets and planetesimals oftentimes intersect. The impact processes can be so energetic that they can produce partial or even total melting and/or vaporization of the bodies involved. The resulting ejecta gather to form a disk. If the temperatures reached during the peak of the impacts exceed the conditions of the critical point (CP), the constituting materials of those celestial bodies become supercritical. The resulting disks would then be monophasic. Upon cooling of a supercritical disk, or in case the supercritical conditions are not reached in the disk during the impact, the system evolves toward a biphase system along the liquid-vapor equilibrium curve. As the maximum of the liquid-vapor equilibrium is fixed by the CP, the position of the CP itself is of fundamental importance in understanding large and giant planetary impacts. And even though giant impacts may be considered rare events, they can be responsible for the creation of planets and moons, or sometimes for their destruction.

Currently, the most widely accepted hypothesis for the formation of our Moon is one such giant impact: a Mars-sized impactor, called Theia, collided with the proto-Earth. This giant impact melted, vaporized, and rendered supercritical a significant portion if not all of the impactor and of the proto-Earth [1–5], creating a large accretion disk from which the Moon formed. Using constitutive equations of the materials involved in the impact, smoothed-particle hydrodynamics (SPH) simulations [6] can describe some of the complex aspects of these giant impacts, and predict the outcome of such planetary impacts as well as the formation of protolunar disks. These SPH simulations are large consumers of data from shock equations of state, like SESAME [7] or ANEOS [8]. Additional meaningful data consist of supercritical points, equations of state of supercritical fluids, and liquid-vaporization equilibria.

Magnesium oxide, MgO, is one of the fundamental building blocks of rocky planets, being a ubiquitous component that appears in most complex silicate minerals. It adds up to roughly 38% [9] of the Earth's and Moon's composition. Deep inside rocky telluric planets, a distinct layer may develop where (Mg,Fe)O forms a mineral by itself, namely magnesiowüstite. In the Earth this layer, containing also (Mg,Fe)SiO₃ bridgmanite, corresponds to the lower mantle, which is the largest part of our planet by volume.

Moreover, MgO is an prototype of countless ionic AB diatomic compounds, many of which are minerals or technological materials. Its face-centered-cubic B1 structure, stable at ambient conditions, represents the most symmetric and most common structure of AB materials.

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*tim.bogels@ens-lyon.fr</footnote>

<footnote>†caracas@ipgp.fr</footnote>
Because of its relevance for both planetary sciences and materials science, MgO has been studied extensively, both experimentally and theoretically, over a wide range of pressures and temperatures. At ambient pressure, the melting point lies at 3125 K and 3.6 g/cm³, and the boiling point is at 3870 K. The first numerical predictions of the melting line [10,11] overestimated the temperatures, while the first experimental results [12] underestimated them. Modern molecular-dynamics simulations predicted the melting of the B1 phase to occur at 3100 K and 0 GPa and at 9400 K and 240 GPa [13].

The high-pressure and high-temperature region was investigated extensively in the past [14–19]. A B1-B2 phase transition was confirmed in an experimental setting [16,17] and in more recent numerical studies [19,20]. The pressure-volume equation of state was measured up to 600 GPa, and the temperature and optical reflectivity to beyond 1400 GPa and 50 000 K in shock experiments [16]. These conditions, while not relevant for the Earth, can be reached inside super-Earth exoplanets, where this phase transition might induce further layering in the rocky mantles of these planets [18].

At the other side of the phase diagram, the liquid-vapor equilibrium line and the position of the critical point are so far almost unknown. This is due to experimental difficulties of sampling both low densities and high temperatures. Theoretical [21] and experimental (molecular beam epitaxy and vacuum thermogravity apparatus) [22,23] studies suggest that the vaporization of MgO crystals is a congruent process, the resulting gas obtained from heating MgO crystals being formed of a stoichiometric mixtures of Mg and atomic O gas. This behavior is expected because of the simple chemistry and stoichiometry of magnesia. But these studies stop short at relatively low temperature, and they fail to reach the CP and to explore the supercritical regime. For the large majority of rock-forming minerals, such as feldspars [24], MgSiO₃, and Mg₂SiO₄, the vaporization is incongruent [21,24–26].

Here we characterize the high-temperature low-density region of MgO using *ab initio* molecular dynamics (MD) in the density functional theory (DFT) framework. These conditions cover the conditions of the liquid spinodal and the critical point and expand into the supercritical regime. The paper is organized in four main parts. Following the introduction, Sec. II details the methodology, the simulations, and the postprocessing. Section III discusses in detail the results: the position of the critical point, the structure of the fluids, the transport properties, the vibrational spectrum, and the electronic atomic charges. The paper ends with a short discussion and conclusions.

## II. METHODOLOGY

### A. First-principles molecular dynamics

We study MgO over a broad range of thermodynamic conditions that cover the low pressures and high temperatures characteristic of the liquid side of the liquid-vapor equilibrium dome. For this, we perform first-principles molecular dynamics (MD) calculations based on density-functional theory (DFT) using the Vienna Ab initio Software Package (VASP) implementation [27–31].

We employ the planar augmented wave function (PAW) [32] flavor of the DFT, with standard PAW pseudopotentials for Mg and O, $3s^2$ with a core radius of 1.06 Å for Mg, and $2s^22p^4$ with a core radius of 0.80 Å for O. We use the Perdew-Burke-Ernzerhof formalism of the generalized gradient approximation [33] for the exchange-correlation functional. The mass of the thermostat was set such that the temperature fluctuations have approximately the same frequencies as the typical phonon-frequencies of MgO, and it was not adapted to the density. The temperature of the system is controlled with a Nosé-thermostat [34]. We employ a kinetic energy cutoff of the plane waves of 550 eV, and we sample the reciprocal space in the $\Gamma$ point. These parameters ensure a precision of the calculations for the energy on the order of 5 meV/at, and for the pressure on the order of 2 kbar.

As is customary in molecular-dynamics simulations of fluids, the systems are modeled using cubic simulation boxes, which are periodically repeated along the three directions of the space. We start the simulations from a $3 \times 3 \times 3$ supercell of MgO with B1 structure, with lattice parameter $a = 4.211$ Å [35]. These supercells contain 216 atoms, i.e., 108 formula units. We heat this supercell up to 5000 K using a heating rate of 0.5 K/fs. We monitor the diffusion of the atoms and find that at this temperature and density, MgO is in a fluid state, i.e., the self-diffusivity of atoms is finite and positive. After thermalization for 1 ps at 5000 K, this configuration constitutes the starting point of our simulations. We go from one temperature to another as needed to follow the different isotherms, increasing or decreasing the temperature in steps of 0.5 K/fs. We sample the density space by changing the unit-cell parameter in steps of 1 Å. The production runs at any given point in pressure and temperature are started after allowing for a thermalization period of 0.5–1 ps. The average duration of the production simulations is on the order of 20 ps for the lower temperatures and higher density systems, and it decreases to a minimum of 5 ps for the high-temperature low-density systems, which are considerably more computationally heavy.

To ascertain the magnitude of the finite-size effects, we ran tests on systems of several densities at 5000 K, with 64 and 512 atoms. The results of these tests can be found in the supplemental material. We find that already at system sizes of 64 atoms, the results start to converge. Our simulations with 216 atoms yield almost the same result, in terms of energy, pressure, and equation of states, as the simulations with 512 atoms. This is consistent with the analysis of finite-size effects in forsterite (Mg₂SiO₄) melts [26], where a system size of 56 atoms already approaches the critical temperature within the accuracy range presented in this research. Furthermore, employing the van der Waals correction for the dispersive forces at the smallest density investigated here caused no appreciable difference for bond lifetimes or pressure.

### B. Finding the position of the critical point

We compute the variation of the pressure as a function of density along various isotherms. We start with the high-density simulations, where the stable phase is the liquid. As we decrease the density, the pressure continuously decreases.
Under enough stretching, the pressure passes below the liquid-vapor equilibrium value; the liquid becomes metastable. The metastable region is where the vapor is thermodynamically stable, but the liquid is mechanically stable. Under further expansion, the pressure reaches a minimum, which marks the liquid spinodal. At lower densities, the liquid is unstable and a gas fraction spontaneously separates from the liquid. As the density continues to decrease, the pressure builds up due to the gas phase until the system reaches a local maximum of the pressure, which marks the gas spinodal [36–38]. Between the local pressure minimum and the local pressure maximum, the fluid is a mechanical mixture of gas and liquid, as they are both unstable as a single phase. Upon further expansion, at densities lower than that of the gas spinodal, the gas becomes metastable. The pressure starts to decrease, passes again the gas-liquid equilibrium pressure, and then asymptotically decreases to zero under infinite expansion. At densities lower than the gas-liquid equilibrium line, the stable phase is the gas. This behavior is best described using the van der Waals gas-liquid equilibrium model, which employs third-order polynomials. We approximate this model with a standard cubic polynomial least-squares fit relating the pressure to the density in accordance with other studies [3,39] (see the supplemental material for more details [40]). The cubic function allows us to quickly find the local minimum and maximum of the polynomial, which yields, respectively, the liquid and the gas spinodal points. The curvature of the fit is influenced by the finite-size effect [41]. Upon increasing the system size, the amplitude of the curves will decrease, and the pressures will all become positive. Figure 1 illustrates the aforementioned features of the polynomials and the construction of the spinodal lines. The liquid and the gas spinodal lines intersect in the critical point.

As the liquid is stable at high densities, the MD simulations allow for the determination of the liquid spinodals at all temperatures. But because of limited ergodicity, the simulations cannot be reliably run at very low densities. This prevents covering the gas spinodal points at low temperatures. Only close to the critical point we can extend the simulations over a density range that encompasses both the gas and the liquid spinodals. Consequently, we approximate the van der Waals model with a cubic function of the pressure as a function of density. Then, as stated above, the local minimum yields the liquid spinodal, and close to the critical point, the local maximum yields the gas spinodal. Hence the position of the critical point is bracketed in density between the gas and the liquid spinodal. In temperature, it lies between the last isotherm whose pressure shows a local minimum and a local maximum, and the first isotherm that shows a monotonous decrease of pressure. This procedure was applied with success on various other material in the same phase space [3, 24, 26, 42].

### C. Postprocessing

We perform all of the postprocessing of the *ab initio* MD runs using the Universal Molecular Dynamics (UMD) software package [43]. We are analyzing structural, transport, vibrational, thermodynamic, and electronic properties in the subcritical and supercritical regimes.

<figure><img src="image_3.png" /></figure>

FIG. 1. Construction of the phase stability fields and critical point from pressure-density relations along several isotherms. The maxima and minima of the isotherms coincide with the liquid and gas spinodal points, respectively. The spinodal lines connect these points. The Maxwell construction delimits regions of equal area between the pressure-density curve and the liquid-vapor equilibrium pressure at each temperature. The equal-area regions are represented with hashed fields. The line joining the densities of the vapor, $\rho_v^{eq}$, in equilibrium with the liquid, $\rho_l^{eq}$, at all temperatures is the binodal line. Both the spinodal and the binodal line have the critical point as maximum in common. In the area between these lines (dark gray), the liquid and gas phases are metastable. The light gray areas outside of the binodal lines are where the individual phases are stable, and the light gray area between the spinodal lines is the region where the liquid and gas coexistence is stable. The area above the isotherm of the critical temperature is where the supercritical fluid is stable. Figure adapted from Kobsch and Caracas [24].

#### a. Pair distribution function

The average interatomic bonding and the coordination environment are important structural properties that stem from the analysis of the pair distribution function, commonly referred to as $g(r)$. The pair distribution function describes the relative distribution of atoms as a function of distance. It zeros out at the start since that is the region where atoms are repulsing each other; the distance up to which the distribution remains null defines the interatomic exclusion zone. The first maximum is associated with the most common interatomic distance, oftentimes referred to as the average bond length. The first minimum of the pair distribution function represents the limit of the first coordination sphere. The second minimum is the limit of the second coordination sphere, and so on. In calculations where the liquid is approximated by a periodic box, the applicability of the distribution function is limited by half the size of the edge of the simulation box, in order to avoid artifacts related to the periodicity.

We take the radius of the first coordination sphere as the threshold value of the interatomic bonds: if two atoms lie closer than this radius, they are considered *bonded*. All the bonded ligands to a central atom define the coordination
<figure><img src="image_4.png" /><figcaption>FIG. 2. Variation of the pressure as a function of density for various isotherms (a). The solid lines are cubic function fits. Their local minima and maxima yield, respectively, the liquid and the gas spinodal points, represented with thick black crosses. The spinodal lines are represented with thin dashed black lines. The critical point lies between the liquid and gas spinodals in density, and between the last isotherm where the pressure still reveals local extremes, and the first isotherm where the pressure is monotonously decreasing the decreasing density. For MgO, this places the critical point at 0.45–0.6 g/cm³ in the density range and between 6500 and 7000 K in the temperature range. The corresponding pressures are on the order of 0.1–0.2 GPa (b).</figcaption></figure>

polyhedra. All the bonded atomic pairs define a connectivity graph, building polymers that describe the structure of the liquid.

b. Mean-squared displacement

The mean-squared displacement (MSD) is the square of the average distance that an atom or cluster of atoms travels as a function of time.

It is calculated using Eq. (1), where $N_α$ is the number of atoms of type α, T is the total time of the simulation, $N_{init}$ is the number of initial times (the number of displacements

measured), and τ is the width of the time window,

$$ \text{MSD}_{\alpha}(\tau) = \frac{1}{N_{\alpha}} \frac{1}{N_{\text{init}}} \sum_{i=1}^{N_{\alpha}} \sum_{t=0}^{T/2} [r_{\alpha,i}(\tau + t) - r_{\alpha,i}(t)]^2. \quad (1) $$

The slope of the MSD yields the self-diffusion coefficient,

$$ D_{\alpha} = \lim_{\tau \to \infty} \frac{1}{nt} \text{MSD}_{\alpha}, \quad (2) $$

where n = 2, 4, 6 for one, two, and three dimensions, respectively. A positive slope of the MSD is a clear

<table><caption>TABLE I. Computed pressure values for MgO obtained for each isotherm at various densities.</caption><thead><tr><th rowspan="2">Density</th><th colspan="8">Pressure (GPa)</th></tr><tr><th>4000 K</th><th>5000 K</th><th>6000 K</th><th>6500 K</th><th>7000 K</th><th>8000 K</th><th>9000 K</th><th>10 000 K</th></tr></thead><tbody><tr><td>3.29 g/cm<sup>3</sup></td><td>29.94</td><td>35.95</td><td>41.56</td><td></td><td>47.27</td><td>52.41</td><td>57.59</td><td>62.77</td></tr><tr><td>2.63 g/cm<sup>3</sup></td><td>9.34</td><td>13.20</td><td>17.22</td><td></td><td>21.20</td><td>24.95</td><td>28.55</td><td>32.16</td></tr><tr><td>2.14 g/cm<sup>3</sup></td><td>1.49</td><td>4.06</td><td>6.77</td><td></td><td>9.38</td><td>12.05</td><td>14.66</td><td>17.28</td></tr><tr><td>1.76 g/cm<sup>3</sup></td><td>-0.80</td><td>0.73</td><td>2.32</td><td>3.36</td><td>4.17</td><td>6.10</td><td>8.10</td><td>9.93</td></tr><tr><td>1.24 g/cm<sup>3</sup></td><td>-0.82</td><td>-0.53</td><td>0.13</td><td>0.50</td><td>1.02</td><td>1.96</td><td>3.07</td><td>4.10</td></tr><tr><td>0.90 g/cm<sup>3</sup></td><td>-0.49</td><td>-0.34</td><td>-0.08</td><td>0.14</td><td>0.43</td><td>0.96</td><td>1.57</td><td>2.27</td></tr><tr><td>0.78 g/cm<sup>3</sup></td><td>-0.33</td><td>-0.25</td><td>-0.04</td><td>0.12</td><td>0.32</td><td>0.79</td><td>1.29</td><td>1.85</td></tr><tr><td>0.68 g/cm<sup>3</sup></td><td></td><td>-0.17</td><td>-0.01</td><td>0.10</td><td>0.27</td><td>0.65</td><td>1.07</td><td>1.52</td></tr><tr><td>0.59 g/cm<sup>3</sup></td><td></td><td>-0.15</td><td>-0.01</td><td>0.13</td><td>0.26</td><td>0.56</td><td>0.92</td><td>1.30</td></tr><tr><td>0.52 g/cm<sup>3</sup></td><td></td><td></td><td>0.03</td><td>0.11</td><td>0.22</td><td>0.51</td><td>0.82</td><td>1.14</td></tr><tr><td>0.37 g/cm<sup>3</sup></td><td></td><td></td><td>0.07</td><td>0.12</td><td>0.20</td><td>0.40</td><td>0.60</td><td>0.84</td></tr></tbody></table>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 3. The pair distribution function for Mg-O (a),(c),(e) and O-O (b),(d),(f) at three isotherms: 4000 (a),(b), 6000 (c),(d), and 10 000 (e),(f) K, and several densities. The first maxima yield a good approximation of the average bond distances. The first minima yield the radius of the first coordination sphere. We use this radius further in the manuscript to define the threshold for interatomic bonding.</figcaption></figure>

indication of the fluid nature of the system studied in the simulation.

c. Velocity autocorrelation function

The general expression for a time correlation function, such as the velocity self-correlation function, is shown in

Eq. (3):

$$C(\tau) = \frac{1}{\tau} \sum_{t=0}^{T/2} A(\tau + t) * A(t), \quad (3)$$

where $\tau$ is the width of the time interval, $T$ is the total time of the simulation, and $A$ is a time-dependent variable. Here
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 4. The speciation of MgO$_{x}$ polyhedra at several densities and temperatures. The coordination polyhedra around each atom are obtained using the analysis of the pair distribution function (Fig. 3). The MgO fluid is dominated by MgO$_{5}$ and MgO$_{6}$ at high densities. The coordination decreases sharply as the density decreases toward and passes the spinodal density. Parts (a), (b), (c), and (d) correspond, respectively, to 2.63, 1.24, 0.90, and 0.52 g/cm$^3$ densities.</figcaption></figure>

we study the self-correlation function of the atomic velocities. Just like with the MSD, we can use the velocity correlation function to determine the diffusion coefficient by taking its integral, shown in Eq. (4), where again $n = 2, 4, 6$ for one, two, and three dimensions, respectively,

$$D_{\alpha} = \frac{1}{nm_{i}N} \int_{-\infty}^{\infty} C(t)dt. \quad (4)$$

### d. Bader charge analysis

We apply the atoms-in-molecule approach of the Bader analysis [44] to obtain the static atomic volumes and charges, using a postprocessing code from the Henkelman Group [45–48]. The procedure finds the saddle points of the total electronic charge distribution around each atom, which, when connected, build the zero flux surface of the charge. These surfaces delimit the parts of the volume of the structure that are assigned to each atom. The integrals of the electronic density inside the atomic volumes yield the total negative atomic charge. The atomic charges are obtained after subtract-

ing the positive charge of the nucleus from the total negative charges.

## III. RESULTS AND DISCUSSION

### A. Critical point

We monitor the variation of the pressure as a function of density at several isotherms in the 4000–10 000 K temperature range (see Table I). These temperatures are above the melting point at ambient pressure conditions, and they extend into the gas and supercritical domains. Depending on the isotherm, we cover the 0.37–3.29 g/cm$^3$ density range. We approximate the van der Waals model with a cubic function fit to the pressure-density points, whose local extrema yield the two spinodals.

Figure 2(a) shows the pressure-density relation at all the isotherms considered here. The liquid spinodals are also indicated on the diagrams. The last two isotherms at which a local minimum can be identified are the 6000 and 6500 K isotherms [Fig. 2(b)]. For these two isotherms, the calculations can reliably sample low enough densities, i.e., the local maxima in the pressure variation can be identified, corresponding to
the gas spinodals. In particular, the 6500 K isotherm shows a local maximum around 0.4 g/cm³ and a local minimum around 0.7 g/cm³ density for pressures of about 1 kbar. On the contrary, there is no local minimum or local maximum along the 7000 K isotherm, but only a monotonous decrease of the pressure with decreasing density. Consequently, the position of the critical point can be constrained in temperature by the isotherms 6500 and 7000 K, and in density by the gas and liquid spinodals at 6500 K, i.e., in the 0.45–0.6 g/cm³ density interval. The pressure range corresponding to these temperatures and densities is 1–2 kbars.

The computed position of the critical point of MgO lies at higher temperatures than any phase in the MgO-SiO₂ phase space, which is relevant for the bulk composition of rocky planets. The critical points calculated for MgSiO₃ [25], Mg₂SiO₄ [26], and SiO₂ [38] lie in the 6200–6500 K temperature range and around 0.50 g/cm³ density. The higher temperatures of MgO confirm its refractory character, while its smaller range of density corresponds to the lighter mass of MgO.

## B. Structure of the fluids

Figure 3 shows the pair distribution functions as a function of density at 4000, 6000, and 10 000 K. The average Mg-O bond length, i.e., the first peak of the distribution function, is on the order of 1.96 Å, weakly dependent on temperature and density. Under compression up to 63 GPa, the bond length decreases by only 0.02 Å along the same isotherm. The radius of the first coordination sphere, which is chosen as the *bonding* criterion, shows a larger variability with both density and temperature. Increasing the density reduces the bonding threshold. The range of the threshold increases considerably with increasing temperature. As a reference, at 30 GPa and 4000 K, the Mg-O bond length in the MgO fluid is slightly larger than the Mg-O bond length in pyrolite [49].

At all temperatures, at 2.63 g/cm³ density, the MgO fluid is dominated by MgO₅, with the second most abundant coordination being MgO₆. Decreasing the density changes the dominant species toward MgO₄ and MgO₃, and down to MgO₂ in the supercritical fluid at 0.5 g/cm³. This decrease in coordination is natural, as it accompanies the decompression of the fluid. Increasing the temperature broadens the distribution of the coordination polyhedra. In the liquid, the coordination number of Mg by O is similar to the one encountered in liquid pyrolite [49]. Figure 4 shows the chemical speciation in the MgO fluid at all densities at several temperatures.

The O-O bond distances are sensitive to both density and temperature. At subcritical temperatures, the first maximum of the pair distribution functions lies around 3.5 Å for densities below about 2 g/cm³. At higher densities, there is a clear decrease of the O-O bond distance, which can be directly related to the increase in coordination of the Mg-O polyhedra from MgO₂-3 to MgO₄-5. At high temperatures and at densities below 1.24 g/cm³, the O-O pair distribution function reveals the presence of a peak around 1.3 Å, which corresponds to the characteristic bond length of the O₂ molecule. As observed in various silicate systems, such as feldspars [24] or silica [38], oxygen molecules are present right below the critical temperature as well as in the supercritical fluid.

<figure><img src="image_15.png" /></figure>

FIG. 5. The speciation of fluid MgO as a function of temperature at 0.68 g/cm³. Each dot represents one MgₓOᵧ cluster, and the vertical axis indicates its size, i.e., x + y. At this density, increasing the temperature takes the system from inside the liquid-vapor dome to the supercritical state. The bimodal distribution of cluster sizes is characteristic for a gas + liquid mixture, while a continuous distribution characterizes the supercritical state. The gap between the cluster sizes closes as the system approaches supercritical temperatures; the most stable clusters are always found at the two extremes.

The analysis of the ··Mg-O-Mg-O-Mg·· polymerization in the fluid allows us to separate the gas phase from the liquid phase. Indeed, the fluid is characterized by largely connected [MgOₓ]n clusters, which represent branched polymers of alternating cations, i.e., Mg, and anions, i.e., O. The gas phase shows isolated clusters, of very limited size. Figure 5 shows the population distribution of all the [MgₓOᵧ] polymers and gas clusters, at a density of 0.68 g/cm³ as a function of temperature. At the lowest temperatures, i.e., 5000 and 6000 K, there is a clear separation between two groups of cluster sizes. The highest values of x + y build the liquid MgO phase, and the lowest values build the incipient gas species escaping from the fluid. As the temperature increases, the gap in the distribution of the cluster sizes closes up. This is indicative of the continuous character of the gas to liquid transition as the system reached the supercritical regime. The supercritical feldspars show a similar behavior [24].

The simulations suggest that the dominant species that form the incipient gas of our system are atomic Mg, atomic O, and some MgO and O₂ molecules. However, further quantification of the gas phase, in terms of relative amounts of stable component species, requires a better sampling of the configuration space, which requires both considerably longer and larger simulations and exploring considerably lower densities.

## C. Vibrational spectrum

At ambient conditions, the solid B1 phase of MgO has only one infrared active phonon mode, whose transverse optical (TO) component lies around 380 cm⁻¹, and the longitudinal optical (LO) one is around 700 cm⁻¹, with the bulk of the spectrum in the 300–500 cm⁻¹ frequency range. Infrared reflectivity of B1 [50] shows the presence of a shoulder around
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

FIG. 6. Total vibrational spectra at several densities of fluid MgO are obtained from the velocity-velocity self-correlation function. Only the values computed at the lowest and the highest isotherms are shown. The noisy data, shown in the background, were filtered using the Savitzky-Golay filter [51]. At high density and low temperature, the spectrum shows a broad peak around 300–500 cm⁻¹. The peak is smoothed out with decreasing density and increasing temperature. Parts (a), (b), (c), and (d) correspond, respectively, to 3.29, 2.63, 0.90, and 0.52 g/cm³ densities.

620 cm⁻¹, corresponding to the tail of the LO-TO splitting. Increasing temperature makes this shoulder disappear and shifts the entire region toward lower frequencies.

The computed vibrational spectrum for the fluid at 4000 K and a density of 3.29 g/cm³ (Fig. 6), corresponding to 30 GPa, shows similarities with the B1 phase. There is a broad peak around 300 cm⁻¹ with a broad shoulder at higher frequencies. Compared to the solid, the spectrum of the fluid is, as expected, more smoothed out, with no detailed features, because of the temperature and the variety of the local coordinations. Another important feature is the nonzero component at zero frequency, which is due to the diffusion.

Decreasing the density to 0.5 g/cm³ and increasing the temperature to 10 000 K smoothes out the main vibrationally active region, albeit shifted toward lower frequencies. The details in the spectrum become less and less pronounced. Eventually, at even lower densities and/or higher temperatures, the spectrum should asymptotically become featureless and approach that of an atomic gas.

### D. Transport properties

We determine the MSDs at all volume and temperature points. The resulting MSDs for oxygen and for magnesium at several isochores are shown in Fig. 7. At all conditions, the MSDs of both atomic types show positive slopes, the systems being in a fluid state.

The O atoms travel for longer distances than the Mg atoms over the same amount of time. The differences between the two atoms depend on both temperature and density, being much more pronounced at higher temperatures and densities. At 2.63 g/cm³, i.e., 9 GPa and 4000 K, the Mg and O atoms travel, respectively, 10 and 10 Å over 10 ps. At the same density and 6000 K they travel, respectively, 14 and 15 Å over 10 ps, and at 10 000 K they travel, respectively, 20 and 24 Å over 10 ps. Inside the liquid-gas dome, at 4000 K and 0.90 g/cm³, the Mg atoms travel 13 Å over 10 ps, and the O atoms travel 14 Å over the same amount of time. The net increase is due to the decompression associated with the opening of the nanobubbles. There is also a clear increase
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /><figcaption>FIG. 7. Mean-square displacements calculated at four isochores. (a),(b) 2.63 g/cm³; (c),(d) 1.24 g/cm³; (e),(f) 0.78 g/cm³; (g),(h) 0.52 g/cm³. The differences in lengths are the result of differences in the length of simulations, the shortest being 5 ps and the longest 20 ps. MgO is fluid at all conditions studied here.</figcaption></figure>
<table><caption>TABLE II. Comparison between the diffusion coefficients at 0.78 g/cm³ as a function of temperature, estimated from the slope of the mean-square displacements as a function of time (MSD) and obtained from the velocity autocorrelation function (VA).</caption><thead><tr><th rowspan="2">Temperature</th><th colspan="4">Diffusion coefficient (m²/s)</th></tr><tr><th>Mg (VA)</th><th>O (VA)</th><th>Mg (MSD)</th><th>O (MSD)</th></tr></thead><tbody><tr><td>4000</td><td>4.87 × 10<sup>-8</sup></td><td>4.34 × 10<sup>-8</sup></td><td>3.22 × 10<sup>-8</sup></td><td>3.25 × 10<sup>-8</sup></td></tr><tr><td>5000</td><td>6.36 × 10<sup>-8</sup></td><td>8.59 × 10<sup>-8</sup></td><td>6.69 × 10<sup>-8</sup></td><td>7.43 × 10<sup>-8</sup></td></tr><tr><td>6000</td><td>1.51 × 10<sup>-7</sup></td><td>1.35 × 10<sup>-7</sup></td><td>1.39 × 10<sup>-7</sup></td><td>1.29 × 10<sup>-7</sup></td></tr><tr><td>7000</td><td>2.24 × 10<sup>-7</sup></td><td>2.69 × 10<sup>-7</sup></td><td>1.85 × 10<sup>-7</sup></td><td>2.26 × 10<sup>-7</sup></td></tr><tr><td>8000</td><td>2.83 × 10<sup>-7</sup></td><td>3.32 × 10<sup>-7</sup></td><td>2.37 × 10<sup>-7</sup></td><td>2.50 × 10<sup>-7</sup></td></tr><tr><td>9000</td><td>2.92 × 10<sup>-7</sup></td><td>3.65 × 10<sup>-7</sup></td><td>2.52 × 10<sup>-7</sup></td><td>3.46 × 10<sup>-7</sup></td></tr><tr><td>10000</td><td>2.98 × 10<sup>-7</sup></td><td>4.09 × 10<sup>-7</sup></td><td>2.86 × 10<sup>-7</sup></td><td>3.85 × 10<sup>-7</sup></td></tr></tbody></table>

in distance traveled by the atoms from 6000 to 7000 K for 1.24 g/m³ and 0.78 g/cm³ [Figs. 7(b) and 7(c)], which corresponds to the passage to the supercritical fluid.

The slope of the MSD with respect to time yields the diffusion coefficients. The results obtained from integrating the velocity self-correlation function confirm these results. Table II lists the values obtained from both methods for comparison. Figure 8 shows the diffusion coefficients as a function of density or different isotherms, as obtained from the self-correlation function.

Below the critical temperature, the dependence of the diffusion coefficients displays a clear separation into two regimes; a linear trend on a log scale at higher densities, and a roughly constant diffusion at lower densities. The point of the slope change corresponds to the density at which the first nanobubbles start to nucleate in the system. As the volume of the simulation box increases, it is the density of the entire system, liquid + gas, that decreases. But the density of the liquid is roughly constant. As the liquid is the dominant phase in these systems at these conditions, the diffusion of both Mg and O atoms reflects their behavior in the liquid phase.

As the temperature increases above the supercritical point, the system is monophasic, so the diffusion coefficient reflects the behavior of the atoms in the total homogeneous system. Here the increase of the volume of the simulation box induces a decrease of the density of the entire system. As the atoms lie farther apart, their diffusion continues to increase with decreasing density.

Consequently, the temperature variation of the diffusion coefficients yields another way of quantifying the transition towards the liquid + vapor dome, and/or the passage to supercritical conditions.

### E. Atomic charges

Finally we analyze the atomic charges of all the atoms in our simulations using the atoms-in-molecule approach, as mentioned in Sec. II. We select several snapshots inside the liquid-vapor dome, which show atoms in both liquid and gas phases, and one snapshot from the supercritical phase.

The values of the Bader charges correlate with the coordination number for both Mg and O. This correlation is visible in Fig. 9. The trend indicates that decreasing the coordination numbers makes the atoms more neutral.

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /><figcaption>FIG. 8. Diffusion coefficients for Mg atoms (a) and O atoms (b) as calculated from the velocity autocorrelation function plotted in log scale against the density. Below the critical temperature, the diffusion coefficients exhibit a change of slope corresponding to the passage inside the liquid-vapor dome. Above the critical temperature, the diffusion changes monotonously with density.</figcaption></figure>

The isolated atoms in the gas phase all have charges approaching zero. This suggests that the gas is close to an atomic-gas model and is not ionized. Indeed, these temperatures are far below the first ionization energy of monatomic magnesium, while the monatomic oxygen does not carry a supplementary electronic charge. The atoms in the liquid, which lie on or close to the interface with the cavities, are ionized, but to a significantly lesser extent than the atoms in the bulk. This suggests that the surface of the bubbles tends to become neutral, and does not carry dipoles. Both Mg and O atoms that lie inside the bulk liquid phase have large charges, between 1 and 2 in absolute values, negative for O and positive for Mg. This suggests an ionic liquid.

The charges in the supercritical fluid show a smaller spread than for the subcritical conditions. Their values are close to the nominal values for $Mg^{2+}$ and $O^{2-}$. This suggests that the supercritical fluid preserves the ionic character of the homogeneous liquid.
<figure><img src="image_30.png" /><figcaption>FIG. 9. Atomic charges for all Mg (blue) and O (red) atoms for a representative snapshot inside the liquid-vapor dome (a). The order of the coordination polyhedron around each atom (b). Charge values close to the nominal correspond to highly coordinated atoms, which lie in the bulk liquid. Charge values close to zero correspond to low coordination numbers, as encountered on the liquid-gas interface and in the gas.</figcaption></figure>

IV. CONCLUSION

We explore an as-yet uncharted area of the phase space of MgO. We provide a thorough analysis of the behavior of fluid MgO at the low-density and high-temperature conditions typically occurring around the supercritical point, using *ab initio* molecular dynamics. We apply a wide range of postprocessing tools to describe different facets of the MgO system in the phase space.

We determine the critical point to be in the density range of 0.45–0.6 g/cm³ and between 6500 and 7000 K. That puts it in the refractory category when compared to other major rock-forming compounds. We characterize a series of transport and structural properties and find similarities to other Mg-rich natural fluids [49]. From the atomic charge analysis we ascertain that the bulk liquid and the supercritical fluid are ionic. For conditions inside the liquid-vapor dome, we find that the surface of the gas bubbles does not carry charge, and that any incipient gas is not charged.

Future studies should address the process of vaporization at much larger time and space scales, and they should determine the liquid-vapor equilibrium curve. For this, more computationally efficient methods will need to be applied, such as machine learning potentials, as they are beyond the scope of the present work. Comparison to other B1-B2 diatomic phases would be interesting as well. Our work contributes to the characterization of an archetypal material, which is also one of the primary constituents of all rocky planets, at conditions typically encountered during planetary formation.

ACKNOWLEDGMENTS

We acknowledge support from the European Research Council under EU Horizon 2020 research and innovation program (Grant Agreement No. 681818–IMPACT to R.C.), the Research Council of Norway, Project No. 223272, and access to supercomputing facilities via the eDARI stl2816 grants, the PRACE RA4947 grant, and the Uninet2 NN9697K grant.

[1] R. M. Canup, Dynamics of lunar formation, *Annu. Rev. Astron. Astrophys.* **42**, 441 (2004).

[2] S. Stewart, E. Davies, M. Duncan, S. Lock, S. Root, J. Townsend, R. Kraus, R. Caracas, and S. Jacobsen, in *The Shock Physics of Giant Impacts: Key Requirements for the Equations of State*, AIP Conf. Proc. No. 2272 (AIP, New York, 2020), p. 080003.

[3] Z. Li, R. Caracas, and F. Soubiran, Partial core vaporization during Giant Impacts inferred from the entropy and the

critical point of iron, *Earth Planet. Sci. Lett.* **547**, 116463 (2020).

[4] R. G. Kraus, S. T. Stewart, D. C. Swift, C. A. Bolme, R. F. Smith, S. Hamel, B. D. Hammel, D. K. Spaulding, D. G. Hicks, J. H. Eggert, and G. W. Collins, Shock vaporization of silica and the thermodynamics of planetary impact events, *J. Geophys. Res. E* **117**, E09009 (2012).

[5] R. G. Kraus, S. Root, R. W. Lemke, S. T. Stewart, S. B. Jacobsen, and T. R. Mattsson, Impact vaporization of
planetesimal cores in the late stages of planet formation, Nat. Geosci. **8**, 269 (2015).

[6] J. J. Monaghan, Smoothed particle hydrodynamics and its diverse applications, Annu. Rev. Fluid Mech. **44**, 323 (2011).

[7] J. D. McHardy, An introduction to the theory and use of SESAME equations of state, Tech. Rep., Report No. LA-14503 (Los Alamos National Laboratory, Los Alamos, NM, 2018).

[8] S. L. Thompson, ANEOS—Analytic equations of state for shock physics Codes—Input Manual Report SAND89-2951 (Sandia National Laboratories, 1990), https://www.osti.gov/biblio/6939284.

[9] A. E. Ringwood, Phase transformations and their bearing on the constitution and dynamics of the mantle, Geochim. Cosmochim. Acta **55**, 2083 (1991).

[10] I. Jackson and R. Liebermann, Melting and elastic shear instability of alkali halides, J. Phys. Chem. Solids **35**, 1115 (1974).

[11] E. Ohtani, Melting temperature distribution and fractionation in the lower mantle, Phys. Earth Planet. Inter. **33**, 12 (1983).

[12] A. Zerr and R. Boehler, Constraints on the melting temperature of the lower mantle from high-pressure experiments on MgO and magnesiioüstite, Nature (London) **371**, 506 (1994).

[13] T. Taniuchi and T. Tsuchiya, The melting points of MgO up to 4 TPa predicted based on ab initio thermodynamic integration molecular dynamics, J. Phys.: Condens. Matter **30**, 114003 (2018).

[14] A. J. Cohen and R. G. Gordon, Modified electron-gas study of the stability, elastic properties, and high-pressure behavior of MgO and CaO crystals, Phys. Rev. B **14**, 4593 (1976).

[15] M. J. Mehl, R. E. Cohen, and H. Krakauer, Linearized augmented plane wave electronic structure calculations for MgO and CaO, J. Geophys. Res. **93**, 8009 (1988).

[16] R. S. McWilliams, D. K. Spaulding, J. H. Eggert, P. M. Celliers, D. G. Hicks, R. F. Smith, G. W. Collins, and R. Jeanloz, Phase transformations and metallization of magnesium oxide at high pressure and temperature, Science **338**, 1330 (2012).

[17] F. Coppari, R. F. Smith, J. H. Eggert, J. Wang, J. R. Rygg, A. Lazicki, J. A. Hawreliak, G. W. Collins, and T. S. Duffy, Experimental evidence for a phase transition in magnesium oxide at exoplanet pressures, Nat. Geosci. **6**, 926 (2013).

[18] K. Umemoto, R. M. Wentzcovitch, S. Wu, M. Ji, C. Z. Wang, and K. M. Ho, Phase transitions in MgSiO₃ post-perovskite in super-Earth mantles, Earth Planet. Sci. Lett. **478**, 40 (2017).

[19] F. Soubiran and B. Militzer, Anharmonicity and Phase Diagram of Magnesium Oxide in the Megabar Regime, Phys. Rev. Lett. **125**, 175701 (2020).

[20] R. Musella, S. Mazevet, and F. Guyot, Physical properties of MgO at deep planetary conditions, Phys. Rev. B **99**, 064110 (2019).

[21] B. Fegley, N. S. Jacobson, K. B. Williams, J. M. C. Plane, L. Schaefer, and K. Lodders, Solubility of rock in steam atmospheres of planets, Astrophys. J. **824**, 103 (2016).

[22] J. Vassent, A. Marty, B. Gilles, and C. Chatillon, Thermodynamic analysis of molecular beam epitaxy of MgO(s) II. Epitaxial growth of MgO layers on Fe(001) substrates, J. Cryst. Growth **219**, 444 (2000).

[23] N. Jacobson, N. Ingersoll, and D. Myers, Vaporization coefficients of SiO₂ and MgO, J. Eur. Ceram. Soc. **37**, 2245 (2017).

[24] A. Kobsch and R. Caracas, The critical point and the supercritical state of alkali feldspars: Implications for the behavior

of the crust during impacts, J. Geophys. Res.: Planets **125**, e2020JE006412 (2020).

[25] B. Xiao and L. Stixrude, Critical vaporization of MgSiO₃, Proc. Natl. Acad. Sci. (USA) **115**, 5371 (2018).

[26] J. P. Townsend, G. Shohet, and K. R. Cochrane, Liquid-vapor coexistence and critical point of Mg₂SiO₄ from ab initio simulations, Geophys. Res. Lett. **47**, e2020GL089599 (2020).

[27] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B **47**, 558 (1993).

[28] G. Kresse and J. Hafner, Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements, J. Phys.: Condens. Matter **6**, 8245 (1994).

[29] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. **6**, 15 (1996).

[30] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B **54**, 11169 (1996).

[31] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B **59**, 1758 (1999).

[32] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B **50**, 17953 (1994).

[33] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett. **77**, 3865 (1996).

[34] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. **81**, 511 (1984).

[35] R. M. Hazen, Effects of temperature and pressure on the cell dimension and X-ray temperature factors of periclase, Am. Mineral. **61**, 266 (1976).

[36] P. Aursand, M. A. Gjennestad, E. Aursand, M. Hammer, and Ø. Wilhelmsen, The spinodal of single- and multi-component fluids and its role in the development of modern equations of state, Fluid Phase Equilib. **436**, 98 (2017).

[37] P. Carlès, A brief review of the thermophysical properties of supercritical fluids, J. Supercrit. Fluids **53**, 2 (2010).

[38] E. C. Green, E. Artacho, and J. A. Connolly, Bulk properties and near-critical behaviour of SiO₂ fluid, Earth Planet. Sci. Lett. **491**, 11 (2018).

[39] V. V. Vasisht, S. Saw, and S. Sastry, Liquid-liquid critical point in supercooled silicon, Nat. Phys. **7**, 549 (2011).

[40] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevB.105.064105 for more details on the finite-size effect, van der Waals interaction, and EOS fitting.

[41] K. Binder, B. J. Block, P. Virnau, and A. Tröster, Beyond the Van Der Waals loop: What can be learned from simulating Lennard-Jones fluids inside the region of phase coexistence, Am. J. Phys. **80**, 1099 (2012).

[42] G. Faussurier, C. Blancard, and P. L. Silvestrelli, Evaluation of aluminum critical point using an ab initio variational approach, Phys. Rev. B **79**, 134202 (2009).

[43] R. Caracas, A. Kobsch, N. V. Solomatova, Z. Li, F. Soubiran, and J.-A. Hernandez, Analyzing melts and fluids from ab initio molecular dynamics simulations with the UMD package, J. Visual Exp. (175), e61534 (2021).

[44] R. F. Bader, S. G. Anderson, and A. J. Duke, Quantum topology of molecular charge distributions, J. Am. Chem. Soc. **101**, 1389 (1979).
[45] W. Tang, E. Sanville, and G. Henkelman, A grid-based Bader analysis algorithm without lattice bias, *J. Phys.: Condens. Matter* **21**, 084204 (2009).

[46] E. Sanville, S. D. Kenny, R. Smith, and G. Henkelman, Improved grid-based algorithm for Bader charge allocation, *J. Comput. Chem.* **28**, 899 (2007).

[47] G. Henkelman, A. Arnaldsson, and H. Jónsson, A fast and robust algorithm for Bader decomposition of charge density, *Comput. Mater. Sci.* **36**, 354 (2006).

[48] M. Yu and D. R. Trinkle, Accurate and efficient algorithm for Bader charge integration, *J. Chem. Phys.* **134**, 064111 (2011).

[49] N. V. Solomatova and R. Caracas, Pressure-induced coordination changes in a pyrolitic silicate melt from ab initio molecular dynamics simulations, *J. Geophys. Res.: Solid Earth* **124**, 11232 (2019).

[50] J. R. Jasperse, A. Kahan, J. N. Plendl, and S. S. Mitra, Temperature dependence of infrared dispersion in ionic crystals LiF and MgO, *Phys. Rev.* **146**, 526 (1966).

[51] A. Savitzky and M. J. E. Golay, Smoothing and differentiation of data by simplified least squares procedures, *Anal. Chem* **36**, 1627 (1964).