Melting curve of SiO2 at multimegabar pressures: implications for gas giants and super-Earths

Ultrahigh-pressure phase boundary between solid and liquid SiO2 is still quite unclear. Here we present predictions of silica melting curve for the multimegabar pressure regime, as obtained from first principles molecular dynamics simulations. We calculate the melting temperatures from three high pressure phases of silica (pyrite-, cotunnite-, and Fe2P-type SiO2) at different pressures using the Z method. The computed melting curve is found to rise abruptly around 330 GPa, an increase not previously reported by any melting simulations. This is in close agreement with recent experiments reporting the α-PbO2–pyrite transition around this pressure. The predicted phase diagram indicates that silica could be one of the dominant components of the rocky cores of gas giants, as it remains solid at the core of our Solar System’s gas giants. These results are also relevant to model the interior structure and evolution of massive super-Earths.


S1. Z-METHOD
In atomistic computer simulation in the microcanonical ensemble it is relatively easy to overheat a crystal beyond its melting temperature T m . This puts the system in a metastable state (the superheated solid state) which survives until the final equilibrium temperature reaches a limit T LS known as the critical superheating limit 2 . Once beyond this limit, the superheated state is unstable and melting will spontaneously ocurr during the simulation (if long enough times are used).
For a temperature of the superheated phase T + LS (slightly above T LS ) the equilibrium temperature after melting will drop to T + m (slightly above T m ). A good estimation of T m can thus be obtained if one searches for the lowest T that triggers melting during the simulation and computes the final equilibrium temperature. As it has been observed 1-19 , the system can transit from the superheated state at T + LS to the liquid state at T + m at constant energy and volume; thus it is verified where E S and E L are the energies of the solid and liquid branches, respectively, as a function of temperature and volume. We denote the common energy in Eq. S1 as the energy of melting E S (V ) for a given volume V .
The procedure for the Z method computation of the melting point is then as follows: at a fixed volume V , the (E, T (E)) points from different simulations draw a "Z" shape, as shown in figure S2, hence the name of the method.
Here the states at T LS and T m have the same energy Note that one also can draw the isochore in the pressure-temperature plane (P (E), T (E)), as is shown in figure S3. In these Z-shaped curves the sharp inflection at the higher temperature corresponds to T LS (V ) and the one at the lower temperature to T m (V ). Thus, by computing the temperature T m and pressure P m of the lower inflection point for different volumes one obtains an estimate of the melting curve T m (P m ) for a particular range of pressures. In practice, several microcanonical simulations are needed, each with a different initial kinetic energy applied to the ideal crystalline structure with a fixed volume. Each simulation at energy E i will add a new (E i , T i ) point to the energy-temperature isochore at volume V (equivalently to the pressure-temperature isochore). With enough points, we obtain isochores in the T -P and T -E diagram which are Z-shaped, as shown in figure S3. However, we need to make sure that the point (T m , P m ) is in fact the lowest in temperature for the liquid branch, and this can take a considerable number of runs. One possible test is to verify that our estimate of T m is aligned in energy with our highest estimate of T LS .
It is possible to take advantage of advanced statistical methods to improve the estimations of T LS , T m and E LS from a number of microcanonical simulations, i.e., statistically reconstruct the most probable isochore given a set of points. For a small number of points, simple least squares fitting discards considerable information which can be included by instead implementing a Bayesian estimation procedure, described below.

S2. AB INITIO MD SIMULATION
The FPMD simulations are in the framework of Kohn-Sham density functional theory. We worked under the NVEensemble, using the Born-Oppenheimer molecular dynamics (BOMD) method, as implemented in the VASP code 28 , employing PAW pseudopotential 29 , and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 30 . For the plane wave expasion of the wavefunctions, we used a cutoff energy of 900 eV, and a 1 × 1 × 2 k-point grid to sample the Brillouin zone in the case of cotunnite-type SiO 2 , and Γ-point only for the pyrite and Fe 2 P structures, as used in previous studies 5 .
We checked sistem size dependence, repeating calculations for the Fe 2 P structure in a 243 atoms cells, which we show in Fig. S4. We observe no difference in the predictions of the melting temperature, validating the 72-atoms structure as large enough.
Time convergence is also an important issue that must be addressed. In Fig. S5, we show the evolution of temperature as a function of time for different initial conditions close to the melting temperature, using the pyrite-type structure. The initial velocities are randomly assigned, so that they correspond to the a given temperature, leading to different values of instantaneous temperatures in the simulation. Homogeneous melting occures during the first 200 steps (40 fs) for the sample initiated at 26000 K , but the structure remains as a superheated solid when initiated at 23000 K (Fig. S5, left). The mean temperatures of the superheated solid and the final liquid are T sol = 9821 K and T liq = 8815 K, while the mean pressures are P sol = 375 GPa and P liq = 391 GPa. We observe that, whether we use 5000 steps (1 ps) or 30000 steps (6 ps), the mean temperature does not change, showing that 5000 steps of 0.2 fs is converged enough for our purposes. FIG. S5: Four independent simulation runs for a volume V = 5.06Å 3 /atom (ρ = 6.58 g/cm 3 ) of the pyrite structure: two initiated at 23000 K (left) and two initiated at 26000 K (right).

S3. BAYESIAN ESTIMATION FOR THE Z-METHOD
We describe an isochoric curve (a Z-curve, in the following) by the four-tuple (T LS , T m , C S v , C L v ), where the last two parameters are the specific heat values for the solid and liquid branches, respectively. The energy of melting E LS is not an independent parameter, as can be determined from where Φ 0 is the potential energy of the ideal crystalline structure.
In Bayesian estimation 31 we propose a probability distribution for the unknown parameters of a Z-curve given the actual (E i , T i ) points obtained from n different microcanonical simulations. According to Bayes' theorem, this probability is given by Here P 0 represents the prior probability for the Z-curve and η is just a normalization constant. This prior probability encodes what is known about the possible values of the parameters of the Z-curve before any actual data points are taken into account in the analysis. In our case, the prior is constructed simply by including consistency information that rejects impossible combinations of parameters which generate invalid Z-curves (for instance, we impose that T LS > T m ).
The estimations of T m with their associated error bars, were computed from this probability distribution in Eq. S3 using Monte Carlo Metropolis sampling 32 .
The melting curve thus obtained, and fitted accordingly 33 , corresponds to the blue line in Fig. 1 of the main manuscript.