# Superionic-Superionic Phase Transitions in Body-Centered Cubic H2O Ice

Jean-Alexis Hernandez, Razvan Caracas

► To cite this version:

Jean-Alexis Hernandez, Razvan Caracas. Superionic-Superionic Phase Transitions in Body-Centered Cubic H2O Ice. Physical Review Letters, 2016, 117 (13), 10.1103/PhysRevLett.117.135503. hal-02333289

HAL Id: hal-02333289

https://univ-lyon1.hal.science/hal-02333289

Submitted on 24 Nov 2021

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
Superionic-Superionic Phase Transitions in Body-Centered Cubic H₂O Ice

Jean-Alexis Hernandez\* and Razvan Caracas†

Laboratoire de Géologie de Lyon, UMR CNRS 5276 (CNRS, ENS, Université Lyon1),
École Normale Supérieure de Lyon, 69364 Lyon Cedex 07, France

(Received 1 May 2016; revised manuscript received 23 June 2016; published 21 September 2016)

From first-principles molecular dynamics, we investigate the relation between the superionic proton conduction and the behavior of the O-H···O bond (ice VII' to ice X transition) in body-centered-cubic (bcc) H₂O ice between 1300 and 2000 K and up to 300 GPa. We bring evidence that there are three distinct phases in the superionic bcc stability field. A first superionic phase characterized by extremely fast diffusion of highly delocalized protons (denoted VII'' hereinafter) is stable at low pressures. A first-order transition separates this phase from a superionic VII', characterized by a finite degree of localization of protons along the nonsymmetric O-H···O bonds. The transition is identified in structural, energetic, and elastic analysis. Upon further compression a second-order phase transition leads to the superionic ice X with symmetric O-H-O bonds.

DOI: 10.1103/PhysRevLett.117.135503

The thermodynamic conditions spanned by the planetary environments hosting H₂O ices are so diverse that almost all stable structures from the phase diagram have a chance to exist somewhere. These phases exhibit remarkable diversity up to about 3 GPa with various ordered and disordered molecular structures. Above this pressure the phase diagram is essentially dominated by structures that are built on a body-centered cubic (bcc) sublattice formed by oxygen atoms. At room temperature, under compression, these ices present continuous transition from a molecular crystal (ice VII) to an ionic crystal (ice X) [1–13].

In ice VII, water molecules preserve their identity and are arranged in a semirandom fashion obeying the Bernal-Fowler-Pauling ice rules [14,15], resulting in O-H···O bonds along the half-diagonals of the bcc sublattice. These bonds are characterized by an underlying double-well proton transfer potential. Increasing pressure and/or temperature activates the translation of the proton between the two potential wells, eventually resulting in a distinct H translationally disordered regime called ice VII'. This phenomenon was first predicted by Holzapfel [10] and then confirmed experimentally and computationally [1,3,4,8,11,12,16,17]. The potential barrier is crossed by proton tunneling and thermal fluctuations, which tend to decrease the transition pressure [3,17,18]. Morrone et al. [19] have shown that above 1000 K the proton tunneling becomes negligible in comparison to the thermal fluctuations, allowing us to model the transition classically. At higher pressure, ice VII' transforms progressively in a symmetric ionic state, called ice X, which is stable up to about 400 GPa [20].

This transition sequence—orientational disorder (ice VII) to dynamical disorder (ice VII') to ionic ice (ice X) —was studied at ambient temperature, but few studies

mention it at high temperature [21,22]. With increasing temperatures, the mobility of the protons is enhanced such that H₂O ice reaches a superionic state [23] while the sublattice of oxygen atoms remains fixed. At several thousand kelvin and over several megabars superionic a bcc-face-centered cubic (fcc) transition affects the oxygen sublattice [24–26]. At 0 K and higher pressure, multiple ground-state structures are predicted to be stable [27–29]. Sun et al. [25] have recently shown that some of these structures become superionic at high temperature (the P₂₁/c-SI phase).

Proton diffusion has been quantified in superionic H₂O with a bcc oxygen sublattice, and different bonding regimes have been noticed at 2000 K, suggesting a behavior close to the ice VII to ice X transition [22,30]. However, no clear relation between proton diffusion and these bonding regimes has been shown.

Here, we show from multiple analysis that several distinct phases populate the superionic regime of the bcc arrangement of the oxygen atoms. For this we study the superionic regime and its transition to the nondiffusive regime up to 300 GPa pressure range in the 1300–2200 K temperature range using first-principles molecular dynamics as implemented in the VASP package [31–34]. We employ the projector augmented-wave method approach [35] with the generalized-gradient approximation in the Perdew-Burke-Ernzerhof formalism [36] for the exchange correlation. We checked and confirmed that the correction for van der Waals interactions does not influence the physical properties of ice at the pressures encountered in our simulations. We model ice in a 4 × 4 × 4 supercell containing 128 water molecules. At every volume-temperature point we generate the starting ice VII configuration independently to remove any possible memory effect. We allow 1–2 ps for thermalization and at least
2–8 ps for data production (more details in S1.A of the Supplemental Material [37]). All simulations are performed in an isokinetic NVT ensemble; i.e., the velocities of all the atoms are rescaled at each time step to match the desired temperature.

First, we calculate the diffusion coefficients ($D_i$) of the oxygen and hydrogen atoms from their mean square displacement (details in S1.B [37]). We distinguish the fluid phase ($D_O$ and $D_H > 0$), the superionic regime ($D_O = 0$ and $D_H > 0$), and the nondiffusive regime ($D_O = 0$ and $D_H = 0$). Only the results corresponding to a fixed bcc sublattice of oxygen atoms ($D_O = 0$) are presented here. Figure 1(a) shows the evolution of hydrogen diffusion coefficient as a function of the shortest

<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Diffusion coefficients of the hydrogen atoms ($D_H$) in solid ice ($D_O = 0$) as function of the OO distance. Superionic diffusion regime corresponds to positive $D_H$. The transition pressure to the superionic regime increases with increasing temperature. (b) Variation of the position of H in the O-H···O bond under compression. Black line represents OH = OO/2, corresponding to ice X. (c) Internal energy ($U$) under compression. The lattice parameter $a$ is related to the OO distance as $a = OO \times 2/\sqrt{3}$. The abrupt decrease in $D_H$ in the superionic regime corresponds to the jump in the internal energy and to the beginning of the O-H···O bond symmetrization. This marks the phase transition between the superionic ice VII'' and the superionic ice VII'.</figcaption></figure>

oxygen-oxygen (OO) distance (the corresponding densities and pressures are given in Table S1 [37]). In measurements, electrical conductivities larger than $1 \text{ S} \cdot \text{cm}^{-1}$ are typically considered to characterize the superionic phases [50]. Here we consider that finite and positive $D_H$ indicates the superionic regime. All our values yield conductivity larger than the threshold. The superionic regime expands to higher pressures as the temperature increases. Along all isotherms, the superionic regime is characterized by a transition from high values of $D_H$ ($3.6 \times 10^{-8} \text{ m}^2\text{s}^{-1}$ at 2000 K), corresponding to the low compression, to small values of $D_H$ (less than $1.0 \times 10^{-8} \text{ m}^2\text{s}^{-1}$), corresponding to higher compression. We note that this transition is sharper at lower temperatures. The values of $D_H$ are coherent with the ones found by French et al. [51] in the fluid phase and with the ones of Goldman et al. [22] in the superionic regime at 2000 K in their simulations starting from a solid initial configuration. At higher pressure (i.e., a smaller OO distance), H atoms diffusion stops altogether, marking the normal, nondiffusive regime. As expected, the limit between the superionic regime and the nondiffusive one moves upwards in pressure with increasing temperature.

Next we look at the relative position of the H atom in the O-H···O triplet, using the pair distribution functions and a three-body analysis (details in S2.A [37]). As shown in Fig. 1(b), there are several cases. For large OO distances, the OH distance remains almost constant under compression: this is ice VII''. As the OO distance decreases, the OH distance starts to increase: this is ice VII'. In this range the O-H···O bond is symmetrizing; i.e., the H continuously changes position approaching the fully symmetric state, halfway between the two O atoms. A fully symmetric O-H-O bond is reached for OO distances between 2.31 and 2.35 Å: this is ice X. The beginning of the O-H···O symmetrization, namely, the transition between ice VII' and ice VII'', corresponds to the abrupt decrease of the hydrogen diffusion, marked by dotted lines in Fig. 1(a). The comparison of the degree of symmetrization of the O-H···O bond with the diffusion rate of the H atoms indicates that the transition between ice VII' and ice X occurs in the superionic regime at 1800 K and above.

We perform one further step in our analysis and look at the individual trajectory of each H atom during the entire simulation time to identify specifically the translational and the rotational component of the H diffusion. This analysis shows that in the temperature range of our simulations H translations in between the two neighboring O atoms are always present. The superionic phases are associated and, thus, can be identified by the presence of both rotations of H around the O atoms and translations along the O-H···O triplets. Moreover, the transition to the superionic ice VII' is associated with a brutal increase of the translation rate and a strong decrease of the rotation rate. The actual ratio of translations vs rotations is characteristic to each phase of
<figure><img src="image_3.png" /><figcaption>FIG. 2. Behavior of the H and O atoms along the 1300 K isotherm. (a) Representation of the atomic trajectories in one supercell representing the simulation box during a 2.5-ps long run in ice X ($P = 300$ GPa), in ice VII' (57.3 GPa), in superionic ice VII' (34.5 GPa), and in ice VII'' (34.2 GPa) just after the transition. (b) Projection of the above trajectories along one O—H···O bond. In ice X the H-atom distribution is unimodal and presents a maximum at the center of the O—H···O bond, whereas in ice VII' the distribution is bimodal. In ice VII'', the delocalization of the H atoms is too large to clearly distinguish the O—H···O bond.</figcaption></figure>

ice (S2.B [37]). Ice VII'' is, as expected, characterized by the highest rate of rotations and a linear increase of translation rate as a function of pressure. Superionic ice VII' and superionic ice X still present a non-null rotation rate, whereas nondiffusive ice VII' and X are only affected by translation jumps.

These structural observations are supported by a jump in the internal energy [Fig. 1(c)] associated with the beginning of the O—H···O symmetrization and with the start of the abrupt decrease in protonic diffusion. The jump is sharp at 1300 K and it smooths out at higher temperatures. This jump is a mark of a first-order phase transition, which takes place between a superionic phase with largely delocalized hydrogen atoms (found at lower pressures, hereinafter called ice VII'') and a superionic phase with H atoms in the O—H···O triplets (the superionic regime of ice VII').

Figure 2(a) shows qualitative analysis of the H delocalization by comparing nuclei trajectories at 1300 K in the ice X, the superionic ice VII', and the superionic VII'' phases. In the superionic ice VII', H atoms remain along 〈111〉 directions (i.e., O—H···O bonds) between the diffusive jumps, whereas the H atoms are delocalized all around a given O atoms in the VII'' phase. In the former case the interstitial space is quasiunoccupied by H, but in the latter there is a non-null population of H traversing it. Another way of looking at the H delocalization, via the O—H bond orientation, shows a progressive delocalization of the H atoms over decompression from the ices X and VII' to the superionic ice VII' to the superionic ice VII'' (S2.C [37]). Figure 2(b) represents the projection along one O—H···O bond of all atomic positions sampled during a given simulation. It shows the symmetrization of H distribution when compressing ice VII' up to ice X, and also the inadequacy of a double-well potential model for the VII'' phase—this is because of the too-broad distribution of the

O—H bond orientations away from the O—H···O line, due to the higher degree of delocalization. The different pattern of H delocalization between ices VII' and VII'' supports the distinction of two phases, as observed by the bonding and energetic analysis.

Based on the structural, diffusion, and energetic analysis we propose modifications to the phase diagram as shown in Fig. 3. It displays the proposed stability fields of the four different phases (fluid, VII'', VII', and X) and the extent of the superionicity in the solids. A first-order phase transition separates the superionic phase VII'' from the superionic ice VII'. The stability of the VII'' phase is artificially extended at low pressures towards the fluid phase; this comes from classical overheating effects observed in molecular-dynamics simulations [52]. However the VII'' phase has its own distinct thermodynamic stability field between the melting curve [30,53], shown with the red line in Fig. 3, and the stability field of ice VII'. As such, it should be observed experimentally at low pressures and high temperatures before reaching the melting line. Vibrational spectroscopy, like Raman and/or infrared, could be used to trace these transitions from the change of slope in pressure of the stretching modes of the O—H bonds. This is visible in the analysis of the vibrational density of states (VDOS); for example, the softening of the OH-stretching during the transition between ice VII' and ice X (see S3.A [37]). Moreover, the finite component of the VDOS at zero frequency is related to the diffusion coefficients, present only for the superionic phases. We find that the transition to ice X is quasi-temperature-independent in the nonsuper-ionic regime; it takes place around 100 GPa, the same transition pressure as obtained at low temperature when the H atoms are treated classically [4,18]. In the superionic regime the transition pressure increases with increasing temperature.
<figure><img src="image_4.png" /></figure>

FIG. 3. Proposed phase diagram obtained from our structural, diffusion, and energetic analysis. Circles indicate the thermodynamic conditions sampled in this study. The colors of the circles correspond to the various phases and/or regimes identified in the phase diagram: light green, ice VII"; cyan, ice VII' in the superionic regime; blue, ice VII' in the nondiffusive regime; magenta, ice X in the superionic regime; and white, ice X in the nondiffusive regime. Black dotted lines correspond to the phase transitions lines. The fluid–ice VII" and ice VII"–ice VII' transitions are first order. The pale blue background area shows the thermodynamic conditions in which we found the superionic regime. The upper red line is the melting curve observed by Schwager et al. [53], as confirmed by calculations [30]; the lower red line is the melting curve as obtained by spectroscopic measurements [6,54,55].

The existence of a VII" phase and the transition between ice VII' and ice VII" are obvious when computing the pressure variation of the elastic properties. Figure 4 shows the elastic constants at several temperatures as a function of pressure. In ices X and VII', all elastic constants decrease with decreasing pressure (increasing OO distance). The slopes are different for each constant, and they are not affected by the development of the superionic regime. $C_{12}$ approaches $C_{11}$ consistently with the development of elastic instabilities observed at 0 K [20,56]. The transition between the VII" phase and ice VII' appears very clearly in the diagram: there is a dip in the pressure variation of the $C_{11}$ and $C_{12}$ elastic constants, exactly at the transition pressure, at all temperatures. The dip shows a clear change in the elastic behavior of H$_2$O ice during the transition from the VII" phase into ice VII'. Similar behavior is observed in other elastic phase transitions, like stishovite and hollandite [57,58]. In the transformations of H$_2$O ice, because the oxygen sublattice is not affected, the elastic transition must be only induced by a change of the localization degree of the H atoms in the superionic regime. As the OO distance increases, the H becomes more and more delocalized; they are less confined to a linear O–H···O triplet. They can move more freely around the O atoms while still preserving the weak H···O bonds that hold the structure together. This is the mark of the proposed ice VII" phase.

Therefore, we propose an addition to the phase diagram of H$_2$O ice at low pressures and high temperatures. We

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

FIG. 4. Elastic constants calculated along different isotherms. $C_{11}$, $C_{12}$, and $C_{44}$ are, respectively, represented in blue, red, and green. In the left panels, the stability fields of the different phases are separated with the black vertical lines, and the pale blue area corresponds to the superionic regime. Right panels are zooms of the transition region between the VII" phase and ice VII'. Black arrows show the conditions in which the Born criterion is not respected.

show that between the melting line and the superionic ice VII' there is another superionic phase with distinct physical properties. This new phase, denoted here as ice VII", is characterized by a large delocalization of the H atoms around the O atoms, which are no longer confined along the O–H···O triplet in between the diffusion jumps. This results in a higher rate of protonic diffusion and a constant OH distance under compression. The transition between ice VII' and ice VII" is marked by a change in the H diffusion rate, by a change in the relation between OH and OO
distances, by a jump in energy, and by a dip in the pressure variation of the elastic constants. This transition could be investigated experimentally by either conductivity or Brillouin measurements. Finally, at high pressures we show that ice VII' and ice X can host H superionic diffusion. Thus the complete symmetrization of the O-H···O bond can occur in the superionic regime if the temperature is high enough.

The simulations were carried out on the Occigen machine at the CINES supercomputing facility (Montpellier, France) and on the Curie machine at the CEA TGCC facility (Bruyères-le-Châtel, France) using computational Grant No. eDARIx2015106368.

jeanalexis.hernandez@ens-lyon.fr

razvan.caracas@ens-lyon.fr

[1] K. Aoki, H. Yamawaki, M. Sakashita, and H. Fujihisa, Phys. Rev. B **54**, 15673 (1996).

[2] M. Benoit, M. Bernasconi, P. Focher, and M. Parrinello, Phys. Rev. Lett. **76**, 2934 (1996).

[3] M. Benoit, D. Marx, and M. Parrinello, Nature (London) **392**, 258 (1998).

[4] M. Benoit, A. H. Romero, and D. Marx, Phys. Rev. Lett. **89**, 145501 (2002).

[5] M. Bernasconi, P. L. Silvestrelli, and M. Parrinello, Phys. Rev. Lett. **81**, 1235 (1998).

[6] A. F. Goncharov, V. V. Struzhkin, M. S. Somayazulu, R. J. Hemley, and H.-K. Mao, Science **273**, 218 (1996).

[7] A. F. Goncharov, V. V. Struzhkin, H.-K. Mao, and R. J. Hemley, Phys. Rev. Lett. **83**, 1998 (1999).

[8] R. J. Hemley, A. P. Jephcoat, H.-K. Mao, C. S. Zha, L. W. Finger, and D. E. Cox, Nature (London) **330**, 737 (1987).

[9] R. J. Hemley and P. Dera, Rev. Mineral. Geochem. **41**, 335 (2000).

[10] W. B. Holzapfel, J. Chem. Phys. **56**, 712 (1972).

[11] P. Pruzan, J. C. Chervin, E. Wolanin, B. Canny, M. Gauthier, and M. Hanfland, J. Raman Spectrosc. **34**, 591 (2003).

[12] M. Song, Phys. Rev. B **60**, 12644 (1999).

[13] K. S. Schweizer and F. H. Stillinger, J. Chem. Phys. **80**, 1230 (1984).

[14] J. D. Bernal and R. H. Fowler, J. Chem. Phys. **1**, 515 (1933).

[15] L. Pauling, J. Am. Chem. Soc. **57**, 2680 (1935).

[16] R. J. Nelmes, J. S. Loveday, W. G. Marshall, G. Hamel, J. M. Besson, and S. Klotz, Phys. Rev. Lett. **81**, 2719 (1998).

[17] M. Benoit, D. Marx, and M. Parrinello, Comput. Mater. Sci. **10**, 88 (1998).

[18] Y. Bronstein, P. Depondt, F. Finocchi, and A. M. Saitta, Phys. Rev. B **89**, 214101 (2014).

[19] J. A. Morrone, L. Lin, and R. Car, J. Chem. Phys. **130**, 204511 (2009).

[20] R. Caracas, Phys. Rev. Lett. **101**, 085502 (2008).

[21] M. French and R. Redmer, Phys. Rev. B **91**, 014308 (2015).

[22] N. Goldman, L. E. Fried, I.-F. W. Kuo, and C. J. Mundy, Phys. Rev. Lett. **94**, 217801 (2005).

[23] C. Cavazzoni, Science **283**, 44 (1999).

[24] H. F. Wilson, M. L. Wong, and B. Militzer, Phys. Rev. Lett. **110**, 151102 (2013).

[25] J. Sun, B. K. Clark, S. Torquato, and R. Car, Nat. Commun. **6**, 8156 (2015).

[26] M. French, M. P. Desjarlais, and R. Redmer, Phys. Rev. E **93**, 022140 (2016).

[27] J. M. McMahon, Phys. Rev. B **84**, 220104 (2011).

[28] M. Ji, K. Umemoto, C.-Z. Wang, K.-M. Ho, and R. M. Wentzcovitch, Phys. Rev. B **84**, 220105 (2011).

[29] C. J. Pickard and M. Martinez-Canales, Phys. Rev. Lett. **110**, 245701 (2013).

[30] E. Schwegler, M. Sharma, F. Gygi, and G. Galli, Proc. Natl. Acad. Sci. U. S. A. **105**, 14779 (2008).

[31] G. Kresse and J. Hafner, Phys. Rev. B **47**, 558 (1993).

[32] G. Kresse and J. Hafner, Phys. Rev. B **49**, 14251 (1994).

[33] G. Kresse and J. Furthmiiller, Comput. Mater. Sci. **6**, 15 (1996).

[34] G. Kresse and J. Furthmüller, Phys. Rev. B **54**, 11169 (1996).

[35] P. E. Blöchl, Phys. Rev. B **50**, 17953 (1994).

[36] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. **77**, 3865 (1996).

[37] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.117.135503, which includes Refs. [38–49], for detailed information about first-principles simulations and post-processing of the molecular dynamics trajectories, additional supporting results and tests.

[38] P. Hohenberg and W. Kohn, Phys. Rev. **136**, B864 (1964).

[39] W. Kohn and L. J. Sham, Phys. Rev. **140**, A1133 (1965).

[40] G. Kresse and D. Joubert, Phys. Rev. B **59**, 1758 (1999).

[41] S. Nose, J. Chem. Phys. **81**, 511 (1984).

[42] W. G. Hoover, Phys. Rev. A **31**, 1695 (1985).

[43] M. Parrinello and A. Rahman, Phys. Rev. Lett. **45**, 1196 (1980).

[44] M. Parrinello and A. Rahman, J. Appl. Phys. **52**, 7182 (1981).

[45] J. F. Nye, *Physical Properties of Crystals: Their Representation by Tensors and Matrices* (Oxford University Press, New York, 1985).

[46] M. Born, Math. Proc. Cambridge Philos. Soc. **36**, 160 (1940).

[47] F. Mouhat and F. X. Coudert, Phys. Rev. B **90**, 224104 (2014).

[48] J. G. Kirkwood, J. Chem. Phys. **3**, 300 (1935).

[49] J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B **83**, 195131 (2011).

[50] J. Boyce and B. Huberman, Phys. Rep. **51**, 189 (1979).

[51] M. French, T. R. Mattsson, and R. Redmer, Phys. Rev. B **82**, 174108 (2010).

[52] A. B. Belonoshko, Am. Mineral. **86**, 193 (2001).

[53] B. Schwager, L. Chudinovskikh, A. Gavriliuk, and R. Boehler, J. Phys. Condens. Matter **16**, S1177 (2004).

[54] T. Kimura, Y. Kuwayama, and T. Yagi, J. Chem. Phys. **140**, 074501 (2014).

[55] J.-F. Lin, Geophys. Res. Lett. **32**, L11306 (2005).

[56] B. Journaux, R. Caracas, P. Carrez, K. Gouriet, P. Cordier, and I. Daniel, Phys. Earth Planet. Inter. **236**, 10 (2014).

[57] M. A. Carpenter, R. J. Hemley, and H.-K. Mao, J. Geophys. Res. **105**, 807 (2000).

[58] T. Boffa Ballaran, J. Liu, L. S. Dubrovinsky, R. Caracas, and W. Crichton, Phys. Rev. B **80**, 214104 (2009).