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

HAL
open science

# Metallization of warm dense SiO$_2$ studied by XANES spectroscopy

A. Denoeud, A. Benuzzi-Mounaix, A. Ravasio, F. Dorchies, P. M. Leguay, J.

Gaudin, F. Guyot, E. Brambrink, M. Koenig, S. Le Pape, et al.

► To cite this version:

A. Denoeud, A. Benuzzi-Mounaix, A. Ravasio, F. Dorchies, P. M. Leguay, et al. Metallization of warm dense SiO$_2$ studied by XANES spectroscopy. Physical Review Letters, 2014, 113 (11), pp.116404. 10.1103/PhysRevLett.113.116404. [hal-01561842](https://hal.science/hal-01561842)

HAL Id: hal-01561842

https://hal.science/hal-01561842

Submitted on 12 Jan 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Metallization of Warm Dense SiO$_2$ Studied by XANES Spectroscopy

A. Denoeud,<sup>1,*</sup> A. Benuzzi-Mounaix,<sup>1,3</sup> A. Ravasio,<sup>1,3</sup> F. Dorchies,<sup>2</sup> P. M. Leguay,<sup>2</sup> J. Gaudin,<sup>2</sup> F. Guyot,<sup>5</sup>
E. Brambrink,<sup>1</sup> M. Koenig,<sup>1</sup> S. Le Pape,<sup>6</sup> and S. Mazevet<sup>3,4</sup>

<sup>1</sup>Laboratoire pour l'Utilisation des Lasers Intenses (LULI), Ecole Polytechnique, CNRS, CEA, UPMC, 91128 Palaiseau, France

<sup>2</sup>Centre Lasers Intenses et Applications (CELIA), CNRS, CEA, Université Bordeaux 1, 33405 Talence, France

<sup>3</sup>LUTH, Observatoire de Paris, CNRS, Université Paris Diderot, 92195 Meudon, France

<sup>4</sup>Département de Physique Théorique et Appliquée, CEA, 91680 Bruyère-le-Chatel, France

<sup>5</sup>Institut de Minéralogie et de Physique des Milieux Condensés (IMPMC), MNHN, CNRS, UPMC, IRD, Sorbonne Universités, 75005 Paris, France

<sup>6</sup>Lawrence Livermore National Laboratory, Livermore, California 94550, USA

(Received 21 February 2014; published 12 September 2014)

We investigate the evolution of the electronic structure of fused silica in a dense plasma regime using time-resolved x-ray absorption spectroscopy. We use a nanosecond (ns) laser beam to generate a strong uniform shock wave in the sample and a picosecond (ps) pulse to produce a broadband x-ray source near the Si K edge. By varying the delay between the two laser beams and the intensity of the ns beam, we explore a large thermodynamical domain with densities varying from 1 to 5 g/cm³ and temperatures up to 5 eV. In contrast to normal conditions where silica is a well-known insulator with a wide band gap of 8.9 eV, we find that shocked silica exhibits a pseudogap as a semimetal throughout this thermodynamical domain. This is in quantitative agreement with density functional theory predictions performed using the generalized gradient approximation.

DOI: 10.1103/PhysRevLett.113.116404

PACS numbers: 71.30.+h, 31.15.ae, 52.25.Os, 52.50.Lp

Understanding the behavior of matter under extreme conditions of density and temperature is a major issue for various fields of physics, such as the design of new materials, inertial confinement fusion, or the modeling of planetary interiors [1–3]. Planetary science is witnessing a revolution with the discoveries of hundreds of planets orbiting nearby stars. Characterizing such astrophysical objects requires the physical properties of key constituents at multi-Mbar pressures and few-eV temperatures [4,5].

It is thus necessary to establish their equations of state and to identify structural changes, phase transitions, metallization, and dissociation processes at planetary interior conditions. The relevant thermodynamical conditions largely extend beyond those anticipated for the solar system planets and require expanding our knowledge of the physical properties at more extreme conditions. It currently presents a challenge for high-pressure physics, as such conditions are not directly accessible experimentally.

In this context, current planetary models rely almost exclusively on physical properties obtained using first principles simulations based on density functional theory (DFT) predictions [6–9]. This approach has known limitations such as band gap underestimation, underbinding, or overbinding of molecular bounds that all potentially affect the predictions at extreme conditions [10]. It is thus important to validate the basic underlying mechanisms occurring for key planetary constituents as pressure and temperature both increase.

Among these constituents, complex silicates [i.e., (Mg,Fe)SiO$_3$ and (Mg,Fe)$_2$SiO$_4$] and the products of

their dissociation (SiO$_2$, MgO) are of major interest. Their metallization and dissociation as pressure and temperature increase are a central issue for the modeling of terrestrial planets or the cores of giant and icy giant planets [11–14]. These two effects influence strongly the equation of state (EOS) and the evolution of the associated transport properties. To increase our confidence in the predictions at planetary interior conditions, we present here an experimental investigation of the metallization of silica in pressure and temperature.

Silica is an insulator with a large band gap of 8.9 eV at normal density-temperature conditions [15]. It turns metallic as density and temperature both increase. This was recently investigated using laser shock experiments by performing reflectivity, temperature, and shock velocity measurements along the principal Hugoniot [16,17]. These studies suggested a connection between the loss of correlation between the Si and O atoms, the closing of the band gap, and the increase in reflectivity.

To go beyond such indirect measurements, we present here an experimental and theoretical study of the properties of fused silica at Mbar pressures using time-resolved X-ray Absorption Near Edge Spectroscopy (XANES). We show that shocked silica is not an insulator with a closing band gap as pressure and temperature increase but exhibits a pseudogap and should thus be considered as a semimetal [18].

The experiment was performed both on the LULI2000 laser facility at the Laboratoire Pour l'Utilisation des Lasers Intenses at the Ecole Polytechnique and on the TITAN laser of the Jupiter Laser Facility at the Lawrence Livermore
<figure><img src="image_3.png" /><figcaption>FIG. 1 (color online). Experiment setup, shocked target geometry, and corresponding SOP typical image. The emission of the four transparent layers can be clearly followed.</figcaption></figure>

National Laboratory. The experimental scheme shown in Fig. 1 follows a previous setup used for aluminum to study electronic structure modifications along the shock Hugoniot [19] and on unloaded states where the nonmetal metal transition was evidenced [20].

A long-pulse laser (500 ps) was focused on a target to launch a shock wave to bring the sample to high pressures. The beam was spatially smoothed by a hybrid phase plate (HPP) producing a flat top focal spot of 600 µm diameter to remove large-scale intensity modulations and obtain a uniform shock. The laser energy at 532 nm was varied between 5 and 120 J, corresponding to intensities between 0.2 and $7 \times 10^{13}$ W/cm$^2$ to change the shock strength and to probe the $SiO_2$ sample over a large domain of density and temperature conditions. The main target consisted in 2-µm fused silica coated by a plasma-enhanced chemical vapor deposition (PECVD) process on 5 µm of free-standing CVD diamond, surrounded by 5 µm of plastic coated by CVD. The accuracy on the thicknesses was obtained by ellipsometry (±0.03 µm on $SiO_2$) and by interferometry (±0.05 µm on plastic and diamond). A layer of graphite (1 µm) was coated on the ablator to avoid direct transmission of the laser light through the rear-side diagnostics.

A second laser beam that delivered 15 J at $2\omega$ in 5 ps with a 50-µm diameter focal spot on an Ytterbium back-lighter target created a spectrally flat x-ray emission from the M band at around 1.85 keV [21]. This corresponds to the Si K edge energy. The short x-ray source duration (5 ps) allowed us to probe $SiO_2$ in a quasistatic state. Varying the delay between the two laser beams (uncertainty <50 ps) allowed us to probe $SiO_2$ under different thermodynamical conditions on the expansion adiabat.

We used two independent KAP conical crystals as x-ray spectrometers to obtain the time-resolved XANES spectra on image plates detectors as described in [22]. The measured spectra result of the convolution of the absorption spectra with an instrumental spectral response of 3-eV full-width at half maximum (FWHM) estimated by using an Si K edge spectrum measurement at ambient pressure and

temperature. The accuracy on the shift measurement is limited by the spectral shot-to-shot stability due to the detector replacement and was estimated at one pixel, corresponding spectrally to 0.4 eV. The configuration of the two detectors, which were placed out of the crystal focal plane, provided a spatial resolution of 25 µm, which allows us to resolve the shocked area with high precision.

Two rear-side diagnostics were used to control the thermodynamical state of the sample at the probed time. The VISARs system (velocity interferometer system for any reflector) [23] can possibly measure the shock velocity in the $SiO_2$ sample and in the last plastic (CH) layer when the generated shock pressure exceeds the metallization threshold that occurs around 1 Mbar in these two materials [17,24,25]. For pressures below 1 Mbar, the shock velocity was deduced from SOP data (Fig. 1). This SOP system (streak optical pyrometer) collected the time resolved optical radiation emitted from the transparent layers at the shock transition. This provided further constraints for analysis and the consistency between the two rear-side diagnostics is checked above 1 Mbar. Since the $SiO_2$ EOS is well known at these conditions [16,17], these measurements coupled to hydrodynamic simulations (multi) allowed us to deduce the thermodynamical conditions probed. We used the 7361 SESAME EOS table [26,27], which is in a very good agreement with our *ab initio* calculations in the range of densities and temperatures investigated here.

Figure 2 summarizes the thermodynamical conditions probed. This represents a large domain of temperatures (from 0.5 to 5 eV), densities (from 1 to 5 g/cm$^3$), and pressures (up to 5 Mbars). The error bars shown in Fig. 2 take into account the numerical laser intensity used to

<figure><img src="image_4.png" /><figcaption>FIG. 2 (color online). On the left, density and temperature conditions inferred from the shock timing and velocity SOP/VISAR data coupled to the hydrodynamic simulations. The different colors refer to different laser intensities. The pink dots represent the conditions at which the *ab initio* simulations were performed. The orange shaded areas represent the conditions at which the measurements and calculations were compared. The right-hand figures show a typical $\rho$ and T spatial profile along the shock direction of 2-µm probed $SiO_2$ as given by the Lagrangian code MULTI.</figcaption></figure>
reproduce the measured shock timings (15%) and the precision of the timing probe (50 ps). This large number of data allowed us to compare the experimental XANES spectra with the *ab initio* simulations for varying temperatures and along two different isochores, at 2 g/cm³ near solid density ($\rho_0 = 2.2$ g/cm³) and at 5 g/cm³.

We performed *ab initio* calculations at conditions representative of the region investigated experimentally, as indicated by pink dots in Fig. 2. These calculations consist in quantum molecular dynamics simulations based on DFT within the generalized gradient approximation (GGA) as parameterized by the PBE (Perdew, Burke, and Ernzerhof) functional [28]. We used the ABINIT electronic structure code [29] and PAW pseudopotentials (Projected Augmented Wave [30]), generated using the Atompaw generator [31]. They consisted in four outer electrons (3s²3p²) for Si and six (2s²2p⁴) for O and a cutoff radius of 1.5a$_B$ and 1.1a$_B$, respectively. We checked that the increase of the number of outer electrons and the decrease of the cutoff radius did not affect the results of the calculations at the highest densities investigated.

The simulations were performed using 162 atoms in the simulation cell. The resulting spectra were converged in terms of particles number by performing simulations using larger cells. The initial configuration was propagated for up to 0.5 ps using time steps of 1 fs in the isokinetics ensemble, where the velocities are rescaled at every time step to keep the temperature constant. This ensemble produces identical results for the hydrodynamical quantities calculated here than the canonical ensemble and is used in this work for convenience, as it is computationally more efficient [32,33]. We obtained the XANES spectra at a given condition by performing calculations within the impurity model for 10 to 15 snapshots along each trajectory as described in [34,35].

As pointed out in [35], we also need to go beyond the frozen-core approximation and obtain the variation of the 1s orbital energy with density and temperature to calculate the variation of the *K* edge. In the current work, we obtain this variation by building an all-electron pseudopotential for the absorbing Si atom and using it within a molecular dynamics run for each of the conditions investigated. This correction was then applied to each spectrum to adjust the position of the *K* edge with respect to the cold experimental XANES spectrum. Furthermore, each calculated spectrum was convoluted by a same Gaussian function of 3 eV corresponding to the experimental resolution. As shown in Fig. 2, we performed calculations on a regular density or temperature grid along the 2, 4, and 5 g/cm³ isochores to interpret the experimental data. This regular grid further indicates marginal variation of the XANES spectra across the density/temperature range corresponding to the quoted experimental error bars at each condition.

In Fig. 3, we compare the calculated and experimental spectra along the 2 g/cm³ isochore and for varying

temperatures. We first note a good agreement at standard density-temperature conditions (black curves). The slope, shape, and the maxima that are linked to the Si-O bounds in the tetrahedral structure are well reproduced by the impurity model [36]. We also note that it is essential to use the impurity model to obtain the right magnitude for the first maximum. This result is also in good agreement with synchrotron data and validates both the experimental and theoretical approaches used here to study the evolution of the XANES spectrum in density and temperature [37].

Figure 3 also shows a clear redshift of about 6 eV in the experimental spectrum as temperature increases to 0.75 eV. It is also found in the calculated spectrum but is underestimated by 1.5 eV. To understand this redshift, we show, in Fig. 4, the corresponding densities of states. The origin of the horizontal axis was set at the Fermi level. Figure 4 shows that starting from the lowest temperature measured, the large band gap found for cold SiO$_2$, is replaced by a pseudogap with a significant metallic component. Calculations indicate that the appearance of this metallic component is gradual and continuous starting around the melting temperature. It follows the gradual loss of correlations between the Si and O atoms [38]. Hence, we find

FIG. 4 (color online). On the left, total densities of states at 2 g/cm³: the gap of the insulator at 0 K disappears and becomes a pseudogap with the temperature. On the right, densities of state at 5 g/cm³: occupied (full line) and total (dashed line). The populated states above the Fermi level increase with the temperature, like in a metal.

<figure><img src="image_5.png" /><figcaption>FIG. 3 (color online). Si *K* edge XANES spectra measured (left) and calculated (right) near solid density.</figcaption></figure>

<figure><img src="image_6.png" /></figure>
that metallization of SiO₂ is due to a pseudogap forming with the increase in temperature. This contrasts with the usual picture of a band gap reducing as pressure and temperature increase. Indeed, it is well known for an atomic system like helium [39,40] that the band gap gradually decreases as pressure and temperature increase but remains finite up to metallization. It can thus be considered as a semiconductor through this density/temperature domain. In contrast, we find here for silica that a broad pseudogap forms as temperature increases. As such, silica in the regime investigated here should be considered as a semimetal rather than a semiconductor since a pseudogap remains present and broad up to the highest temperatures and densities (Fig. 4).

We can now interpret the behavior of the XANES spectra as temperature increases. We recall that the XANES spectrum near the *K* edge is proportional to the unoccupied partial *p*-projected density of states (*p*-DOS). In this regime, we further check that the partial *p*-DOS is the dominant contribution around and above the Fermi level and within the pseudogap. In temperature, the *K* edge energy corresponds to the transition from the Si 1 s orbital to the Fermi level in the pseudogap instead of the conducting band in the insulator at standard conditions. The redshift observed in Fig. 3 corresponds directly to the energy difference between these two states. Its experimental value is also less important than the 8.9-eV experimental band gap of cold silica. Indeed, the Fermi level of the semimetal does not correspond to the valence band of the insulator. In addition, the difference of 1.5 eV between the experimental and calculated spectra can be traced back to the well-known underestimation of the band gap by DFT approaches [10]. Indeed, as observed in Fig. 4, GGA calculations give a band gap value of 5.3 eV at standard conditions, that is underestimated by 3.6 eV. This is an issue only at normal conditions as a metallic component appears, filling the gap in temperature rather than a gap gradually closing.

Moreover, as temperature further increases to 2.2 eV, the theoretical predictions indicate a decrease of the slope and a reduction of the first maximum above the edge. In contrast, the experimental spectrum remains only slightly affected by this additional increase in temperature. We further check that neither the slight mismatch in densities nor the impurity model can explain this disagreement at the highest temperature.

We can now turn to the behavior at still higher temperatures obtained experimentally at higher densities. As indicated in Fig. 2, we compared the theoretical spectra along the 5 g/cm³ isochore with measurements obtained at densities varying between 4 and 5 g/cm³. Figure 4 shows the variation of the corresponding occupied density of states as temperature increases up to 6 eV along this isochore. Similar to what we found at lower density, silica also exhibits a pseudogap with a significant metallic

contribution. Its behavior in temperature is dominated by the metallic contribution and is driven by the pseudogap occurring near the Fermi level. The occupied density of states follows what is expected for a simple metal: the increase in temperature results in an increase of the density of states around the Fermi energy and in the population and depopulation above and below the Fermi level, respectively. We recover this behavior in the measured XANES spectra.

Figure 5 shows the XANES spectra obtained when increasing the temperature. On this isochore, we see that the maxima are shifted toward the blue side when compared with the result at lower densities, with some differences between the calculated and measured spectra. This discrepancy could come from some limitations in the calculations at reproducing the changes in the DOS as density increases (this point will be the subject of a future article). Despite this discrepancy, the simulations recover the overall behavior of the XANES spectrum observed experimentally. As can be anticipated from the occupied density of states, the XANES spectrum exhibits a broadening of the edge as temperature increases. This behavior in temperature is similar to what is found for a simple metal like aluminum [19]. This metallic component and the associated pseudogap also drive the transport properties. The good agreement between the calculations and the experimental data also brings confidence in the transport properties calculated using a ground-state approach as used in the current work.

In summary, we studied the electronic structure of silica in a wide thermodynamical range and a yet unexplored region of the phase diagram. The comparison of the measured XANES spectra with *ab initio* simulations provides an overall satisfying agreement. It demonstrates the appearance of a pseudogap with rising temperature and as soon as the Si and O atoms lose correlation. This shows that silica in this regime is not a system with a band gap reducing with density and temperature with a conduction band populated by temperature but rather behaves as a semimetal.

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

FIG. 5 (color online). Variation in temperature of the Si *K* edge XANES spectrum around the 5 g/cm³ isochore given by the calculations (right) and observed experimentally (left). The spectra are normalized to the x-ray absorption value far above the *K* edge and offset for clarity.
A first implication of this result is that the shortcomings of standard DFT methods regarding band gap underestimation and metallization pressure do not apply for warm dense silica. In contrast to what has been speculated until now, GW or TDDFT calculations commonly used to correct DFT approaches for the optical property of semiconductor in the solid state will not be of much use here. This also probably explains the good agreement previously found for the reflectivity data of shocked silica calculated using standard DFT approaches [38]. This situation probably occurs in other molecular systems in the warm dense matter regime and should be systematically explored. In conclusion, the investigation of the gap closure mechanism has important implications on the reliability of ab initio calculations to predict metallization onset conditions.

We acknowledge the expert support from the technical groups at LULI2000 and JUPITER facilities, from the GENCI program for providing computational time, and from the french RENATECH network and its FEMTO-ST technological facility for target fabrication. This work was supported by the ANR Project PLANETLAB (ANR-12-BS04-0015). The CELIA team thanks the ANR OEDYP (ANR-09-BLAN-0206-01) for equipment and travel support.

<footnote>*adrien.denoeud@polytechnique.edu</footnote>

<footnote>[1] G. Collins et al., Science **281**, 1178 (1998).</footnote>

<footnote>[2] M. D. Knudson, D. L. Hanson, J. E. Bailey, C. A. Hall, J. R. Asay, and W. W. Anderson, Phys. Rev. Lett. **87**, 225501 (2001).</footnote>

<footnote>[3] M. Koenig et al., Appl. Phys. Lett. **72**, 1033 (1998).</footnote>

<footnote>[4] G. Huser, M. Koenig, A. Benuzzi-Mounaix, E. Henry, T. Vinci, B. Faral, M. Tomasini, B. Telaro, and D. Batani, Phys. Plasmas **12**, 060701 (2005).</footnote>

<footnote>[5] J. Eggert, S. Brygoo, P. Loubeyre, R. S. McWilliams, P. M. Celliers, D. G. Hicks, T. R. Boehly, R. Jeanloz, and G. W. Collins, Phys. Rev. Lett. **100**, 124503 (2008).</footnote>

<footnote>[6] I. Baraffe et al., Protostars and Planets VI (University of Arizona Press, Tucson, Arizona, 2014).</footnote>

<footnote>[7] B. Militzer, Phys. Rev. B **87**, 014202 (2013).</footnote>

<footnote>[8] N. Nettelmann, A. Becker, B. Holst, and R. Redmer, Astrophys. J. **750**, 52 (2012).</footnote>

<footnote>[9] D. Alfè, G. Price, and M. Gillan, Phys. Rev. B **65**, 165118 (2002).</footnote>

<footnote>[10] S. Faleev, M. van Schilfgaarde, T. Kotani, F. Léonard, and M. Desjarlais, Phys. Rev. B **74**, 033101 (2006).</footnote>

<footnote>[11] K. Umemoto, Science **311**, 983 (2006).</footnote>

<footnote>[12] D. K. Spaulding, R. S. McWilliams, R. Jeanloz, J. H. Eggert, P. M. Celliers, D. G. Hicks, G. W. Collins, and R. F. Smith, Phys. Rev. Lett. **108**, 065701 (2012).</footnote>

<footnote>[13] J. Tsuchiya and T. Tsuchiya, Proc. Natl. Acad. Sci. U.S.A. **108**, 1252 (2011).</footnote>

<footnote>[14] R. S. McWilliams, D. K. Spaulding, J. H. Eggert, P. M. Celliers, D. G. Hicks, R. F. Smith, G. W. Collins, and R. Jeanloz Science **338**, 1330 (2012).</footnote>

<footnote>[15] T. H. DiStefano and D. E. Eastman, Solid State Commun. **9**, 2259 (1971).</footnote>

<footnote>[16] M. Knudson and M. P. Desjarlais, Phys. Rev. Lett. **103**, 225501 (2009).</footnote>

<footnote>[17] D. G. Hicks, T. R. Boehly, J. H. Eggert, J. E. Miller, P. M. Celliers, and G. W. Collins, Phys. Rev. Lett. **97**, 025502 (2006).</footnote>

<footnote>[18] E. J. Reed, M. R. Manaa, L. E. Fried, K. R. Glaesemann, and J. D. Joannopoulos, Nat. Phys. **4**, 72 (2008).</footnote>

<footnote>[19] A. Benuzzi-Mounaix, F. Dorchies, V. Recoules, F. Festa, O. Peyrusse, A. Levy, A. Ravasio, T. Hall, M. Koenig, N. Amadou, E. Brambrink, and S. Mazevet, Phys. Rev. Lett. **107**, 165006 (2011).</footnote>

<footnote>[20] A. Lévy et al., Phys. Rev. Lett. **108**, 055002 (2012).</footnote>

<footnote>[21] M. Harmand, F. Dorchies, O. Peyrusse, D. Descamps, C. Fourment, S. Hulin, S. Petit, and J. J. Santos, Phys. Plasmas **16**, 063301 (2009).</footnote>

<footnote>[22] A. Levy, F. Dorchies, C. Fourment, M. Harmand, S. Hulin, J. J. Santos, D. Descamps, S. Petit, and R. Bouillaud, Rev. Sci. Instrum. **81**, 063107 (2010).</footnote>

<footnote>[23] L. M. Barker and R. E. Hollenbach, J. Appl. Phys. **43**, 4669 (1972).</footnote>

<footnote>[24] M. Koenig, F. Philippe, A. Benuzzi-Mounaix, D. Batani, M. Tomasini, E. Henry, and T. Hall, Phys. Plasmas **10**, 3026 (2003).</footnote>

<footnote>[25] C. Wang, X.-T. He, and P. Zhang, Phys. Plasmas **18**, 082707 (2011).</footnote>

<footnote>[26] G. I. Kerley, *Equations of State for Composite Materials* (Kerley Albuquerque: NM, 1999).</footnote>

<footnote>[27] SESAME, Report No. LA-UR-92-3407 (1992).</footnote>

<footnote>[28] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. **77**, 3865 (1996).</footnote>

<footnote>[29] X. Gonze et al., Comput. Mater. Sci. **25**, 478 (2002).</footnote>

<footnote>[30] P. E. Blöchl, Phys. Rev. B **41**, 5414 (1990).</footnote>

<footnote>[31] N. A. W. Holzwarth, A. R. Tackett, and G. E. Matthews, Comput. Phys. Commun. **135**, 329 (2001).</footnote>

<footnote>[32] M. P. Allen and D. J. Tildesley, *Computer Simulation of Liquids* (Clarenton Press, Oxford, 1987).</footnote>

<footnote>[33] D. Frenkel and B. Smit, *Understanding Molecular Simulation: From Algorithms to Applications* (Academic Press, San Diego, 2001).</footnote>

<footnote>[34] S. Mazevet and G. Zerah, Phys. Rev. Lett. **101**, 155001 (2008).</footnote>

<footnote>[35] V. Recoules and S. Mazevet, Phys. Rev. B **80**, 064110 (2009).</footnote>

<footnote>[36] D. Li, G. M. Bancroft, M. Kasrai, M. E. Fleet, R. A. Secco, X. H. Feng, K. H. Tan, and B. X. Yang, Am. Mineral. **79**, 622 (1994).</footnote>

<footnote>[37] G. S. Henderson and J. C. St-Amour, Chem. Geol. **213**, 31 (2004).</footnote>

<footnote>[38] Y. Laudernet, J. Clérouin, and S. Mazevet, Phys. Rev. B **70**, 165108 (2004).</footnote>

<footnote>[39] P. Kowalski, S. Mazevet, D. Saumon, and M. Challacombe, Phys. Rev. B **76**, 075112 (2007).</footnote>

<footnote>[40] F. Soubiran, S. Mazevet, C. Winisdoerffer, and G. Chabrier, Phys. Rev. B **86**, 115102 (2012).</footnote>