Path integral Monte Carlo and density functional molecular dynamics simulations of warm dense ${\mathrm{MgSiO}}_3$

In order to provide a comprehensive theoretical description of ${\mathrm{MgSiO}}_3$ at extreme conditions, we combine results from path integral Monte Carlo and density functional molecular dynamics simulations and generate a consistent equation of state for this material. We consider a wide range of temperature and density conditions from ${10}^4$ to ${10}^8\mathrm{K}$ and from 0.321 to 64.2 $\mathrm{g}{\mathrm{cm}}^{-3}$ (0.1- to 20-fold the ambient density). We study how the L and K shell electrons are ionized with increasing temperature and pressure. We derive the shock Hugoniot curve and compare with experimental results. Our Hugoniot curve is in good agreement with the experiments, and we predict a broad compression maximum that is dominated by the K shell ionization of all three nuclei while the peak compression ratio of 4.70 is obtained when the Si and Mg nuclei are ionized. Finally we analyze the heat capacity and structural properties of the liquid.

Figure 1

Temperature-density conditions of our DFT-MD and PIMC simulations along with computed isobars, isentropes, and three principal shock Hugoniot curve that were derived for an initial density of ${\rho }_0=3.207911\mathrm{g}$ ${\mathrm{cm}}^{-3}$ $(V_0=51.965073\AA /\mathrm{f}.\mathrm{u}.)$.

Figure 2

Temperature-pressure conditions for the PIMC and DFT-MD calculations along isochores corresponding to the densities of 0.1-fold (0.321 g ${\mathrm{cm}}^{-3}$) to 20-fold (64.20 g ${\mathrm{cm}}^{-3}$). The blue dash-dotted line shows the principal Hugoniot curve of ${\mathrm{MgSiO}}_3$ obtained from our simulations, using an initial density of ${\rho }_0=3.207911\mathrm{g}$ ${\mathrm{cm}}^{-3}$ $(V_0=51.965073\AA /\mathrm{f}.\mathrm{u}.)$. The red dashed line corresponds to the Hugoniot curve from Ref. [85], calculated from DFT-MD simulations. Experimental measurement of the principal Hugoniot curve from Ref. [86], an isentrope derived from this experiment (solid green line), and the Hugoniot curve for ${\mathrm{MgSiO}}_3$ glass [28] (orange region) are shown for reference. The melting line of ${\mathrm{MgSiO}}_3$ derived from two-phase simulations [87] is shown in the dashed gray line, while the melting curve derived from shock experiments [86] is represented by the thick black line.

Figure 3

Excess pressure as a function of temperature relative to the ideal Fermi gas, computed with PIMC, DFT-MD, and the Debye-Hückel plasma model. The results are plotted for densities of (a) 6.4, (b) 3.651, (c) 7.582, and (d) 15.701 g ${\mathrm{cm}}^{-3}$.

Figure 4

Excess internal energy, relative to the ideal Fermi gas, computed with PIMC, DFT-MD, and the Debye-Hückel plasma model. As in the corresponding Fig. 3, the results are plotted for densities of (a) 6.4, (b) 3.651, (c) 7.582, and (d) 15.701 g ${\mathrm{cm}}^{-3}$ as a function of temperature.

Figure 5

Total internal energy as a function of density, computed with PIMC and DFT-MD.

Figure 6

Pressure as a function of density, computed with PIMC and DFT-MD.

Figure 7

Grüneisen parameter of ${\mathrm{MgSiO}}_3$ as a function of volume for different temperatures. The horizontal, dashed line represents the high-temperature limit (the Grüneisen parameter of the ideal gas), ${\gamma }_0=2/3$. Values of the Grüneisen parameter along the principal Hugoniot from shock compression experiments [86] (shown in gray circles with error bars) correspond to pressures between 230–380 GPa and temperatures between 6200–10 000 K.

Figure 8

Integrated nuclear-electron pair correlation, functions $N\left(r\right)$, computed with PIMC simulations. The three columns correspond to the Mg, Si, and O nuclei, while the four rows show results for the densities of 0.321 (0.1-fold compression), 3.21 (1.0-fold compression), 12.83 (4.0-fold compression), and 64.16 g ${\mathrm{cm}}^{-3}$ (20-fold compression). For given density, results are shown for the same set of temperatures for all three nuclei. However, these temperatures adjusted with increasing density because the degree of ionization is density dependent. $N\left(r\right)$ represent the average of number of electrons contained within a sphere of radius, $r$, around a given nucleus. For comparison we show the corresponding functions with thin dashed lines for isolated nuclei with double occupied $1\mathrm{s}$ core states that we computed with the GAMESS software [93].

Figure 9

$N\left(r\right)$ and $g\left(r\right)$ correlation functions for liquid ${\mathrm{MgSiO}}_3$ at $T=50523\mathrm{K}$ and $\rho =19.247$ $\mathrm{g}{\mathrm{cm}}^{-3}$ (sixfold compression). In the three lower panels, the $g\left(r\right)$ functions were split into the contributions from the $n\mathrm{th}$ nearest neighbors. The functions of the third, ninth and 21st neighbors were shaded because their locations respectively correspond to first maximum, first minimum, and second maximum of the Si-O $g\left(r\right)$ function. Their locations are also marked by the vertical dotted lines. To improve the readability, the scales of the $Y$ axes in the three lower panels were split into two separate linear parts, one from 0 and 1, and a compressed region from 1 and 2.5.

Figure 10

Nuclear radial distribution functions computed with DFT-MD simulations of liquid ${\mathrm{MgSiO}}_3$ at a fixed temperatures of $50523\mathrm{K}$ and $202095\mathrm{K}$. Functions are calculated in 65-atom cells and compared for densities of (from top to bottom) 6.41 (twofold, top), 16.04 (fivefold), 19.25 (sixfold), and 32.08 g ${\mathrm{cm}}^{-3}$ (10-fold).

Figure 11

Comparison of the principal shock Hugoniot curves of ${\mathrm{MgSiO}}_3$ with those of carbon [44], oxygen [92], neon [83], aluminum [18], and silicon [19, 45]. The two dash-dotted curves and shaded regions show Hugoniot curves without any electronic excitations and just without K shell ionization. Electronic excitations increase the shock compression. Without them, a maximum compression ratio of 4.7 could not be attained.

Figure 12

Comparison of principal and secondary shock Hugoniot curves with isotherms and isentropes in pressure-density space.

Figure 13

Left: particle velocity $u_p$ is shown as a function of compression ratio. Right: The deviations from a linear $u_p-u_s$ fit. The coefficients are $A=1.2764$ and $B=8.2315$ km/s.

Figure 14

Heat capacity and Grüneisen parameter temperature dependence, derived from DFT-MD and PIMC calculations at various densities. The horizontal dashed line in the upper panel represents the high-temperature limit of $C_V=\frac{3}{2}N{k}_B$ where $N=55$ is the total number of free particles per ${\mathrm{MgSiO}}_3$ formula unit, with five ions and 50 electrons. Values above this line mark the temperature region of ionization where the internal energy increases significantly. In the lower panel, the Grüneisen parameter of the ideal gas, ${\gamma }_0=2/3$, also represents the high-temperature limit.