Thermal and nuclear quantum effects in the hydrogen bond dynamical symmetrization phase transition of \delta-AlOOH

Thermal and nuclear quantum effects in the hydrogen bond dynamical symmetrization phase transition of delta-AlOOH


We conducted ab initio molecular dynamics simulations of the phase of the hydrous mineral AlOOH at ambient temperature and high pressure. Nuclear quantum effects were included through a Langevin dynamics in a bath of quantum harmonic oscillators. We confirm that under increasing pressure -AlOOH undergoes a phase transition from a structure with asymmetric and disordered O-H bonds to a stiffer phase with symmetric hydrogen bonds, which should be stable within the pressure and temperature ranges typical for the Earth’s mantle. The transition is initially triggered by proton tunneling, which makes the mean proton position to coincide with the midpoint of the O-O distance, at pressures as low as 10 GPa. However, only at much larger pressures, around 30 GPa as previously found by other calculations, the phase with symmetric hydrogen bonds is stable from the classical point of view. The transition is also characterized through the analysis of the H-O stretching modes, which soften considerably and fade out around 10 GPa in the structure, when thermal and nuclear quantum effects are taken into account in the simulations. At variance, the harmonic picture is not adequate to describe the highly anharmonic effective potential that is seen by the protons at the transition. Finally, we propose that the picture of a dynamical transition to the high-symmetry and proton-centered phase, which is brought about by the onset of proton tunneling, could be confirmed by quasi-elastic neutron scattering and vibrational spectroscopy under pressure.

I Introduction

Water transport in the Earth deep lower mantle is a cenral issue in geology geological studies (Ohtani et al., 2001). In particular, hydrous minerals are believed to play a major role as water vectors (Panero & Stixrude, 2004; Sano et al., 2004, 2008; Ohira et al., 2014). Among the various dense hydrous minerals, aluminium oxide hydroxide (AlOOH) has received much attention recently, in particular its phase, which is a high-pressure polymorph of the more known diaspore (-AlOOH) and boehmite (-AlOOH). -AlOOH has been synthesized for the first time in 2000 at GPa and C (Suzuki et al., 2000). Since then, it has been shown to be stable in a wide range of pressure and temperature (Ohtani et al., 2001; Sano et al., 2004, 2008), making it a candidate as a hydrogen and water carrier within the deep mantle. Its structure is known through diffraction experiments (Suzuki et al., 2000; Ohtani et al., 2001; Vanpeteghem et al., 2002; Sano et al., 2004, 2008; Komatsu et al., 2006; Kuribayashi et al., 2014) and theoretical calculations (Tsuchiya et al., 2002; Panero & Stixrude, 2004), except for the positions of the hydrogen atoms, which are not easily accessible experimentally. Moreover, isotope effects deriving from the differences between -AlOOH and -AlOOD (Sano-Furukawa et al., 2009) indicate that nuclear quantum effects (NQE) play a role, as it happens in other minerals(Umemoto et al., 2015; Hermann & Mookerjee, 2016).

Several studies suggested that the ionocovalent and hydrogen bonds in -AlOOH become symmetric under increasing pressure, in contrast to -AlOOH (Winkler et al., 2001; Friedrich et al., 2007). Recent synchrotron X-ray diffraction measurements indicated that there is a phase transition between and GPa, which is characterized by a modification of the compression properties of -AlOOH (Kuribayashi et al., 2014). However, this technique was unable to specify whether the post-transition phase consists in disordered and asymmetric hydrogen bonds, or if they are indeed symmetric.

From the theoretical viewpoint, ab initio calculations predicted that -AlOOH undergoes a second-order phase transition from an asymmetric (hydrogen off-centered - HOC) with group symmetry to a symmetric (hydrogen centered - HC) structure with group symmetry at about GPa (Tsuchiya et al., 2002; Tsuchiya & Tsuchiya, 2009). Ab initio phonon dispersion curves computed at K showed that there are strong modifications of the modes involving protons at the transition (Tsuchiya et al., 2008). However, the previous theoretical studies included neither thermal effects nor proton tunneling, which would be likely to reduce the transition pressure as in the case of the VII-X transition of ice under high pressure (Benoit et al., 1998; Bronstein et al., 2014).

However, the scenario for the transition in -AlOOH seems to be more involved: in contrast with ice VII, it is not clear whether before the transition the protons are ordered or disordered (Komatsu et al., 2006; Sano-Furukawa et al., 2008; Tsuchiya et al., 2008); moreover, -AlOOH belongs to the space-group, a much lower symmetry than ice VII, with two inequivalent protons per unit cell in a less symmetric atomic environment. It is not a priori obvious that the same symmetric double-well picture stands for -AlOOH as well. In particular, contrarily to the high-pressure ice structure, the oxygen-oxygen midpoint is not a symmetry center in -AlOOH, whereas it is in the space-group; finally the experimental transition pressure (Kuribayashi et al., 2014) ( 10 GPa) is much lower than that predicted by static density functional theory calculations (Tsuchiya et al., 2008) ( 30 GPa): dynamical effects in the symmetrization of the H-bonds may be relevant to characterize the transition at the atomistic scale.

We carry out, for the first time in -AlOOH to our knowledge, ab initio molecular dynamics simulations at ambient temperature, in the - GPa pressure range. The quantum nature of the nuclei is taken into account by coupling their dynamics to a ”Quantum Thermal Bath” (QTB).(Dammak et al., 2009) From these simulations, we extract structural and dynamical quantities in order to investigate the hydrogen behavior in -AlOOH and disentangle thermal from quantum effects as pressure increases. We anticipate that both thermal fluctuations and NQE play a crucial role in the modification of hydrogen bonds in -AlOOH, contributing to a substantial downshift of the theoretical transition pressure, in agreement with the the experimental observations.

After presenting in section II the simulation techniques, including the QTB method, in section III, we confront our simulation results with known structural data from mainly X-ray scattering. In section IV, we compute vibrational properties that yield a clear signature of the phase transition at the transition pressure and discuss the link with experiments.

Ii Methods

ii.1 The Quantum Thermal Bath

For our simulations, we use the Quantum Espresso package (Giannozzi et al., 2009), in which we implemented the QTB method. We perform molecular dynamics (MD) simulations in which a random force is added as well as a friction term according to a Langevin-type equation. When the nuclei are treated as classical particles (corresponding to standard MD), the random force simulates the thermal motion of the system (Kubo et al., 1978), while in the QTB method, it simulates both thermal and quantum fluctuations (Callen & Welton, 1951; Dammak et al., 2009). The equations of motion for the nuclei are in both cases:


where is the position of a nucleus at time , its mass and the interatomic force applied on nucleus and computed via the Density Functional Theory (DFT). is the random force, while is a friction coefficient. The frequency distribution of the random force is related to the friction coefficient through its power spectrum (Callen & Welton, 1951), where denotes the Fourier transform of and the function characterizes the kind of bath, either classical (standard Langevin MD) or quantum (QTB MD):


In practice, the only difference between standard Langevin MD and QTB MD is that in the standard case, has a Gaussian distribution with standard deviation ( is the temperature and the Boltzmann constant), while in the QTB case, the random force is a colored noise. Therefore, from the comparison between standard Langevin MD and QTB trajectories, we are able to disentangle thermal from purely quantum effects. From now on, we denote with ”standard Langevin MD” classical Langevin ab initio (DFT-based) molecular dynamics (thus using a white noise), while ”QTB MD” denotes Langevin ab initio molecular dynamics using the colored noise in Equation 3 and thus including NQE.

The QTB is exact for a system of harmonic oscillators, but also provides a good approximation for anharmonic systems (Dammak et al., 2012; Calvo et al., 2012; Savin et al., 2012; Bronstein et al., 2014, 2016). It has recently been the object of many studies that partly assessed its advantages and drawbacks (Calvo et al., 2012; Basire et al., 2013; Bedoya-Martinez et al., 2014; Hernandez-Rojas et al., 2015; Brieuc et al., 2016). In systems with few degrees of freedom, the probability distributions obtained via QTB capture most of the relevant quantum effects. Indeed, the computed probabilities along the atomic trajectories within the QTB approach are good approximations of the corresponding quantum distributions (Calvo et al., 2012; Dammak et al., 2012; Bronstein et al., 2014), although a perfect quantitative agreement with exact methods can only be reached for harmonic solids (Dammak et al., 2009; Basire et al., 2013; Brieuc et al., 2016).

With respect to other methods that include NQE such as path-integral-based molecular dynamics (Berne & Thirumalai, 1986; Marx & Parrinello, 1996), the QTB has the advantage of presenting no additional computational cost compared to standard MD and to give access to vibrational information about the system, which is particularly suited for the present case and for more complex mineral structures. Vibrational spectra computed from QTB compare well with experiments and exact results both in molecules, clusters (Calvo et al., 2012), and ice under pressure (Bronstein et al., 2014, 2016). Indeed, the spectra computed from Langevin MD, whether with a classical or a quantum noise, are a convolution of the real physical spectra with a Lorentzian of width , the friction coefficient (Izaguirre et al., 2001). As long as is chosen small enough, the vibrational spectra are reliable: hence, we set THz, which is well below significant O-H frequencies in our system.

ii.2 Technical details

The random force is generated at the beginning of the simulation using a standard numerical technique to generate a noise with a prescribed correlation function (Maradudin et al., 1990). In particular, no knowledge of the eigenfrequencies of the system is required before the simulation. The interatomic forces are computed via the DFT within the generalized gradient approximation (GGA) (Perdew et al., 1996), through the Becke-Lee-Yang-Parr functional (Becke, 1988; Lee et al., 1988). The interaction between the ionic cores and the valence electrons is described through ultra-soft pseudopotentials with non-linear core corrections. The Kohn-Sham orbitals are expanded in plane-waves with Ry energy cutoff. Our simulation cell is a supercell containing AlOOH units. The corresponding Brillouin Zone is sampled by a Monkhorst-Pack grid with a offset (half a grid step in each direction). We run QTB and standard Langevin MD simulations at ambient temperature ( K) for several volumes. The time length of each simulation is about ps with a fs integration time step. The instantaneous pressures are calculated via the stress theorem (Nielsen & Martin, 1985). From the atomic trajectories, we compute the mean pressures, and extract probability distributions and vibrational spectra.

Iii Structural properties

iii.1 versus structures in

First, we address from a static point of view and with no NQE, the apparent contradiction between theoretical (Tsuchiya et al., 2002; Tsuchiya & Tsuchiya, 2009) computations that predict a structure for pressures up to 30 GPa and X-ray experiments (Kuribayashi et al., 2014) which yield a symmetry above 10 GPa.

Static lattice parameter optimization and experimental data
-HOC; Theory -HC; Theory -HOC; Exp.
this work this work Kuribayashi et al. (2014)
Pressure (GPa) 0 0 0.0001
4.7882 4.2833 2.8617 4.7426 4.1944 2.8659 4.7127 4.2250 2.8310
  Al 0 0.029 0 0 0 0 0 0.02433 0
Al 1/2 0.471 1/2 1/2 1/2 1/2 1/2 0.47567 1/2
O 0.341 0.250 0 0.355 0.240 0 0.34246 0.24860 0
O 0.841 0.250 1/2 0.855 0.260 1/2 0.84246 0.25140 1/2
O 0.642 0.747 0 0.645 0.760 0 0.6458 0.74918 0
O 0.142 0.753 1/2 0.145 0.740 1/2 0.1458 0.75082 1/2
H 0.516 -0.059 0 1/2 0 0
H 0.016 0.559 1/2 0 1/2 1/2
-HOC; Theory -HC; Theory -HC; Exp.
this work this work Kuribayashi et al. (2014)
Pressure (GPa) 10 10 8.2
4.6829 4.1835 2.8201 4.6658 4.1404 2.8211 4.6379 4.1342 2.7990
  Al 0 0.020 0 0 0 0 0 0 0
Al 1/2 0.480 1/2 1/2 1/2 1/2 1/2 1/2 1/2
O 0.343 0.246 0 0.353 0.240 0 0.3510 0.2431 0
O 0.843 0.260 1/2 0.853 0.260 1/2 0.8510 0.2569 1/2
O 0.645 0.753 0 0.647 0.760 0 0.6490 0.7369 0
O 0.145 0.747 1/2 0.147 0.740 1/2 0.1490 0.7531 1/2
H 0.513 -0.042 0 1/2 0 0
H 0.013 0.542 1/2 0 1/2 1/2
-HOC; Theory -HC; Theory
this work this work
Pressure (GPa) 30 30
4.5557 4.0641 2.7475 4.5527 4.0572 2.7484
  Al 0 0.008 0 0 0 0
Al 1/2 0.492 1/2 1/2 1/2 1/2
O 0.347 0.242 0 0.350 0.240 0
O 0.847 0.258 1/2 0.850 0.260 1/2
O 0.648 0.759 0 0.650 0.760 0
O 0.148 0.741 1/2 0.150 0.740 1/2
H 0.506 -0.017 0 1/2 0 0
H 0.006 0.547 1/2 0 1/2 1/2
Table 1: Equilibrium lattice parameters (Å) and atomic positions in units of lattice vectors of -AlOOH, computed within GGA-BLYP, through static optimization at 3 pressures( (top), 10 (middle) and 30 GPa (bottom), compared with available experimental data from Kuribayashi et al. (2014).The irreducible atomic positions are in boldface characters.

The relation between the atomic positions in and space groups is not completely straightforward. From direct inspection to the International Tables of Crystallography, the two sets of positions can be continuously changed into one another, by shifting the origin of the atomic positions, all of the kind, in (space group 31) by with respect to the conventional origin.

After the shift, one can easily see that the symmetry group can be continously tranformed into the one (see e.g. table 1: atomic positions from static structure optimization and corresponding experimental data). More specifically, concerning the proton positions which are not easily detectable by X-ray scattering and subject to relevant quantum effects such as tunneling, we denote two kinds of structures, as previously done in the literature (Tsuchiya et al., 2002; Tsuchiya & Tsuchiya, 2009): off-centered Hydrogen (HOC) and centered Hydrogen (HC). While Hydrogen in can only exist as HC, one can build two virtual structures: HOC and HC. In the former, the atomic positions are , while in the latter the Hydrogen positions are identical to the ones in . Nevertheless, as a result of the loss of inversion symmetry, protons in -HC are not at the same distance from the two inequivalent O and O and their energy landscape is not symmetric in the order parameter . For instance, for GPa (see table 1, top) the computed distances are: 1.03 Å and 1.57 Å (-HOC); 1.22 Å and 1.22 Å (-HC).

At this stage, table 1 confirms with our static optimization data the previous ab-initio results by Tsuchiya et al. (2002) and Tsuchiya & Tsuchiya (2009) and obtain quantitative agreement with experimental crystallographic data by Kuribayashi et al. (2014). In addition, consecutive dynamical matrix calculations provide several negative eigenvalues and thus imaginary frequencies for the -HC structure (figure 1) for GPa confirming that the HC structure is unstable up to that pressure in contradiction with experimental data: this means that static ab-initio optimization is unable to account for the experimentally obtained transition pressure and we can suspect that thermal effects and/or nuclear quantum effects play a significant role in this transition. Furthermore, that the two structures can be transformed continuously into one another, as also shown for the high pressure data at 30 GPa in table 1, indeed allows dynamical effects to play a role in the phase transition.

iii.2 Equation of state

Synchrotron X-ray experiments showed a change of compressibility of -AlOOH around GPa (Sano-Furukawa et al., 2009). Moreover, theoretical calculations conducted in the asymmetric ( - HOC) and symmetric ( - HC) phases yielded two distinct bulk moduli (Tsuchiya et al., 2002; Tsuchiya & Tsuchiya, 2009). In the asymmetric structure, up to GPa, our results from standard MD runs agree with those from Tsuchiya et al. (2002). The same procedure was carried out within the QTB frame: we fitted our computed points at K through two distinct Vinet equations of state (Vinet et al., 1987), one adapted for the low-pressure (LP) structure, the other for the high-pressure (HP) structure (see Figure 2), as done in previous studies (Tsuchiya et al., 2002; Sano-Furukawa et al., 2009). Within this procedure, the two equations of state cross between and GPa; accordingly, the bulk modulus of -AlOOH changes from () GPa at low pressure to () GPa at high pressure (in both cases, ). The small discrepancies with the experiments (Sano-Furukawa et al., 2009) are mainly due to the GGA functional, which slightly overestimates the cell parameters with respect to the experimental values. Our dynamical simulations are therefore consistent with a change of compressibility in -AlOOH with an increase of the bulk modulus by about between and GPa at  K, which is compatible with a phase transition towards stiffer O-H-O bonds as indicated by experiments. However, given the statistical error bars, we cannot demonstrate that a single equation of state cannot fit our computed points and we therefore resort to more precise methods, such as the analysis of the interatomic distances and the vibrational spectra, to establish the transition pressure in the following.

iii.3 Interatomic distances

The O-O average distance is a relevant parameter for any transition from asymmetric O-HO to symmetric O-H-O bonds; indeed, in pure ice, proton tunneling and the onset of symmetrization happen at O-O distances around Å  (Benoit & Marx, 2005), and this is generally accepted as a good criterion for the onset of symmetrization. The average O-O distances that were computed along standard Langevin and QTB simulations at  K are reported in Figure 3. They are very similar, apart from larger fluctuations within the QTB framework, which are mainly due to zero-point motion. Our results are in agreement with experimental data (Komatsu et al., 2006; Sano-Furukawa et al., 2008; Kuribayashi et al., 2014), though we slightly overestimate the O-O distance as pressure increases to about GPa with respect to very recent X-ray data (Kuribayashi et al., 2014). One may note that both experimental results and ours for the O-O distance, remain significantly above the 2.4 Å  threshold until over 25 GPa.

In order to investigate the hydrogen bonds in -AlOOH, we compute the oxygen-hydrogen (O-H) pair correlation function (PCF), averaged over the trajectories obtained from QTB or standard Langevin MD (see Figure 4). At low pressure, the PCF displays two peaks indicating two distinct bond lengths: around Å  for ionocovalent O-H bonds and around Å  for hydrogen HO bonds. These distances are in agreement with experimental measurements on -AlOOD obtained through neutron scattering around GPa (Sano-Furukawa et al., 2008). Since NQE do not seem to affect significantly the positions of the two peaks in the PCF, the distances for hydrogen and deuterium are expected to be similar, even though quantum effects are weaker for deuterium than for protons.

Upon compression, the O-O mean distance decreases (see Figure 3); this is correlated with a dilatation of the O-H bond and a concomitant contraction of the HO length. As pressure further increases, the two peaks merge until they become hardly distinguishable from one another; at high pressures, the ionocovalent and hydrogen bonds cannot be differentiated anymore. However, the PCF from the QTB simulations differ from those obtained via standard Langevin MD: the peaks are broader than their classical counterparts, and they merge at a lower pressure, even though the positions of the peaks are similar. Moreover, the minimum of the PCF between the two peaks is much higher when NQE are taken into account, suggesting that either disorder or proton tunneling is occuring at low pressure. The PCF thus shows a merging of the two peaks, therefore indicating a structural change within the explored pressure range, however actual symmetrization does not occur: at GPa, the oxygen-oxygen distance is larger than 2.4Å  (Figure 3), while the O-H distance at the same pressure is less than 1.1Å, that is, less than half the O-O distance. This means that if a structural phase transition occurs, the PCF never becomes symmetric within the explored pressure range.

iii.4 Proton tunneling

The analysis of the PCF in section III.3 does not allow to distinguish between the intrinsic width of the two peaks on the one hand and proton hopping on the other. The determination of a precise transition pressure from the computed distributions is therefore problematic. The usual relevant proton transfer coordinate in symmetrization transition is where d and d are the O-H distances between a given proton H and its two nearest-neighbor oxygen atoms (Benoit & Marx, 2005; Morrone et al., 2009; Lin et al., 2011). Thus, when the proton is covalently bonded to one oxygen, we have while indicates that the proton is located at the midpoint of the O-O segment.

Figure 5 shows the distribution of the proton transfer coordinate computed in two different ways: from the time-averaged atomic positions and from the instantaneous atomic positions. The time-averaged results show a transition to a symmetric state () at approximately 10 GPa, which is consistent with space group, while the full distribution is not peaked at . This means that the average symmetry is that observed experimentally, but the system displays a large degree of disorder and the local instantaneous configurations may have lower symmetry than the average ones. Moreover, at low pressure, the iono-covalent O-H bond is longer and the hydrogen bond shorter within QTB MD with respect to classical Langevin dynamics, in agreement with the results on the position of the peaks in the PCF (see Figure 4). We note also that our from standard Langevin MD i.e. with virtually classical protons, are in good agreement with experimental measurements on -AlOOD (Sano-Furukawa et al., 2009), the deuterium nucleus being heavy enough to be treated as a classical particle. Moreover, we also compute the evolution of mean with pressure in -AlOOH (diaspore), where no symmetrization occurs in the pressure range studied here (up to GPa), which is in agreement with experimental results (Winkler et al., 2001; Friedrich et al., 2007). In addition, the distinction between results from QTB MD and standard Langevin MD is negligible: in diaspore, nuclear quantum effects do not play any major role in the hydrogen bond behavior, in contrast with the phase.

A closer look to the full proton transfer coordinate distribution that is computed from the instantaneous atomic positions allows to obtain more detailed information on the role of thermal and nuclear quantum effects along the symmetrization of the hydrogen bonds in -AlOOH. Within the classical frame, shows a well defined maximum at for pressures up to GPa, which shifts from Å at GPa to Å at GPa. At larger pressures, the proton transfer coordinate distribution becomes flat in a distance range which contracts when increases. Classical Langevin dynamics thus suggests a transition from a proton-ordered asymmetric hydrogen bond to a dynamically disordered symmetric hydrogen bond taking place around 15-20 GPa, where the variance of the distribution is largest. The picture that is obtained by the dynamics with the quantum thermal bath is in striking contrast with the classical one: the proton transfer coordinate distribution is quite flat even at low pressures, with , which evidences the presence of a small but non null proton tunneling and a large variance - the maximum being Å at GPa, close to its classical counterpart. The distribution flattens and its variance lessens for increasing pressures. Already at GPa, the maximum is hardly distinguishable and .

Therefore, according to the quantum frame, -AlOOH undergoes a transition from a proton-disordered asymmetric hydrogen bond to a proton-disordered symmetric hydrogen bond. By direct comparison between the two frames, quantum effects enhance the proton disorder even at low pressures and slightly downshift the transition pressure with respect to purely thermal effects. The situation is thus somewhat complex: from X-ray scattering, in which proton disorder should be barely visible, only the averaged symmetry can be seen at high pressure. Proton tunneling at the time-scale of the experiment thus triggers the dynamical symmetrization of the hydrogen bond that leads to the hydrogen-centered (HC) phase, which is observed experimentally (Sano-Furukawa et al., 2009; Kuribayashi et al., 2014) for larger pressures than about 8 GPa. Its symmetrization from a classical viewpoint, that is when the potential felt by the protons only has one minimum, occurs at a higher pressure, around 30 GPa as found by Tsuchiya et al. (Tsuchiya & Tsuchiya, 2009), as we will show in the next section that is devoted to the analysis of vibrational modes in -AlOOH.

Iv Vibrational properties

As previously pointed out, the proton distribution is difficult to establish through X-ray diffraction experiments, and pair-correlation functions as obtained from simulations are not always absolutely conclusive. On the other hand, vibrational spectroscopic measurements are an excellent means to detect experimentally phase transitions through mode softening (Aoki et al., 1996; Goncharov et al., 1996; Struzhkin et al., 1997; Goncharov et al., 1999), and it turns out these are also easily accessed through time-dependent simulations such as QTB molecular dynamics (Bronstein et al., 2014; Brieuc et al., 2016).

The vibrational spectra of -AlOOH are rather complex (Tsuchiya & Tsuchiya, 2009). Therefore, in order to determine which phonon modes are relevant to the hydrogen bond symmetrization in -AlOOH, we conduct a preliminary mode analysis via the dynamical matrix, which is computed through density functional perturbation theory (Baroni et al., 2001) at K and without any quantum or anharmonic effects (section IV.1), below and above the transition. Some of the relevant modes from the simulated vibrational spectra (at ambient temperature and including NQE) are then discussed in section IV.2.

iv.1 Phonons: harmonic approximation at K

In the asymmetric hydrogen bond configuration (, HOC) (Tsuchiya et al., 2008), the vibrations at the Brillouin zone center consist in two high-frequency O-H stretching modes (around cm at GPa), and four O-HO bending modes between and cm. The two low-frequency bending modes vibrate in the direction, while the two bending modes at higher frequencies are in the plane, which contains the O-HO bonds. When pressure increases, the O-H stretching mode frequencies decrease abruptly, while the bending mode frequencies increase slightly.

In the symmetric phase (, HC) (Tsuchiya et al., 2008) at , the two highest O-H bending modes (around cm at ) are in the plane and the two lowest ones along vibrate in the - cm range. In contrast, the O-H stretching modes at have imaginary frequencies for pressures below 30 GPa (see figure 1), their modulus tends to zero with increasing . This clearly indicates that the HC phase should be unstable below that pressure. At 30 GPa, the stretching frequencies at the Brillouin zone center become real and then increase with pressure. At GPa they mix with bending modes, in the 1400 - 1700 cm range and harden with further increasing pressure.

As already shown by Tsuchiya et al. (2002, 2008), this is consistent with a second-order transition from the HOC to the HC structure in -AlOOH, which is accompanied by a strong softening of the O-H stretching modes. Therefore, our dynamical matrix analysis at K and without NQE confirms a transition at GPa, which, however is about 20GPa above the experimental transition pressure from Kuribayashi et al. (2014), a far from negligible difference. Moreover, the imaginary frequency O-H stretching mode also show motion of the Aluminium atoms (figure 1) which means that the symmetry of the position thereof with respect to the hydrogen environment is unstable: this is reminiscent of the local symmetry as discussed in section III.4. Finally, for some pressures around 30GPa, the bending and stretching modes can mix. Hence, we cannot simply rely on the description of the vibrations in -AlOOH via the dynamical matrix approach close to the HOC - HC transition. We can indeed expect both thermal and quantum effects to allow the system to explore the anharmonic regions of phase space, as suggested by the issue of local versus average symmetric positions as done in section III.4 for protons. Thus, we turn to the study of the atomic trajectories in order to take into account the inherent anharmonicity of the system, as induced by thermal and quantum effects.

Nevertheless, the dynamical matrix approach allows us to choose a basis of eigenvectors - the phonon directions , - that is a key to the interpretation of the hydrogen vibrational spectra. In the configuration, is parallel to the O-O direction, while is orthogonal to in the plane. In the configuration, the O-H stretching mode eigenvector forms a small angle with the O-O direction, and is orthogonal to in the plane. In both configurations, is parallel to the direction. The eigenvectors from the two phases are quite similar: in particular, and . Therefore, the use of these eigenvectors to obtain mode-specific density of states (see Equation 5) appears to be reliable in both phases. At the transition, the OH modes can mix and become a linear combination of harmonic ones; apart from the modes here reported, there are no others involving proton displacements in the harmonic case. Therefore, the vectors form a minimal basis set to analyze the behavior with pressure of modes involving protons.

iv.2 Vibrational spectra from molecular dynamics trajectories

In order to take into account both thermal and quantum effects, as well as the anharmonicities of the system, the vibrational spectrum of hydrogen in -AlOOH is calculated directly from our MD simulations, through the Fourier Transform of the autocorrelation function involving the different nuclei (Marshall & Lovesey, 1971). For each pressure, the hydrogen spectrum density is proportional to:


where is the shortest O-H distance (corresponding to the O-H ionocovalent bond) and the sum runs over all H atoms. As contains information about all hydrogen vibrations, we also compute the partial spectrum densities in order to unravel the distinct mode contributions:


where is the projection of the hydrogen position vector onto one of the vectors defined in the previous section. Each partial spectrum then corresponds to hydrogen vibrations along the eigenvector . By comparing the full spectrum to the peaks of , we determine the frequencies of the distinct hydrogen modes: the O-H stretching modes and the O-HO bending modes. The evolution of these frequencies with increasing pressure is shown in Figure 6, as well as the width of the corresponding peaks (indicated by the vertical bars).

At low pressure, we distinguish two stretching modes at approximately and cm and two bending modes around cm (in the plane) and cm (in the direction). These results are in very good agreement with Raman measurements (Ohtani et al., 2001) at ambient pressure. Due to the large widths of the peaks in the calculated spectrum, we could distinguish only two high-frequency O-H stretching modes, while Ohtani and coworkers observed four broad bands (Ohtani et al., 2001). Tsuchiya and coworkers suggested that the presence of these multiple bands, instead of two sharp peaks as calculated in the dynamical matrix approach (see section IV.1), is due to hydrogen disorder at low pressure (Tsuchiya et al., 2008), which is consistent with our results on the O-H pair correlation function. Upon compression, the O-H stretching modes soften gradually, their frequencies dropping to approximately and cm. Moreover, the width of the high-frequency mode increases, which implies that the two peaks of merge at the transition and are no longer distinguishable. In contrast, the bending mode frequencies increase slightly with pressure. Above a critical pressure ( GPa), the O-H stretching mode frequency ( cm) does not vary appreciably with pressure, at least up to GPa, while the bending mode frequencies increase continuously. The second O-H bending mode in the plane becomes distinguishable only above GPa. The previous trends agree with those found by Tsuchiya and coworkers (Tsuchiya et al., 2008).

The softening of the O-H stretching mode up to about GPa, where bending and stretching modes mix up, is consistent with a second-order phase transition as also guessed by the discontinuity of the compressibility. This is also quite close to the high-pressure behavior of the O-H vibrations during the transition from asymmetric phase VII to symmetric phase X in ice (Aoki et al., 1996; Goncharov et al., 1996; Struzhkin et al., 1997; Goncharov et al., 1999; Bronstein et al., 2014, 2016). Above this critical pressure, there is a continuous hardening of the O-H-O bending modes, while the stretching mode frequencies do not vary in our pressure range (up to GPa). By extrapolating the trends shown in Figure 6, it appears that the bending modes show higher frequencies than the stretching modes above GPa, as also predicted by the dynamical matrix analysis in section IV.1.

All the previous observations conspire to indicate that GPa is close to the critical pressure for the transition from a proton-disordered 2nm phase with asymmetric hydrogen bonds to a proton-disordered phase with symmetric hydrogen bonds. The proton disorder decreases with pressure and a proton-ordered phase with symmetric hydrogen bonds is recovered at pressures as high as GPa.

V Conclusion

We study the symmetrization of the Hydrogen bond in the high-pressure phase of AlOOH by molecular dynamics simulations, thus including thermal effects. For the first time to our knowledge, we also include nuclear quantum effects in the description of the transition through a recent though relatively simple method - the Quantum Thermal Bath (QTB) - in which semi-classical atomic trajectories fulfilling the Bose-Einstein statistics are generated via the contact with a stochastic bath of quantum harmonic oscillators. From a detailed comparison between three simulation frames (static calculations, classical molecular dynamics and quantum molecular dynamics simulations) we disentangle thermal and quantum effects.

The transition can be schematically depicted as proceeding from a low-pressure hydrogen off-centered (HOC) phase with space group to a more symmetric, hydrogen-centered (HC) phase with space group at high pressure. The HC- phase is harder than the HOC- one, with important consequences on the stability and mechanical properties of -AlOOH in the lower Earth mantle.

The static calculations fully agree with previous ones (Tsuchiya et al., 2002; Li et al., 2006; Tsuchiya et al., 2008) indicating that a pressure as high as GPa is needed to stabilize the proton-ordered HC- phase. The apparent contradiction with recent X-ray diffraction data (Kuribayashi et al., 2014) providing experimental evidence that the transition to the high-symmetry phase happens at much lower pressures, around 8 GPa, is solved by analyzing the role of thermal and nuclear quantum effects in the transition. The onset of proton tunneling appears at rather low pressures, below 10 GPa, which is much lower than that needed to shorten O-O distances down to Å - the usual threshold for the hydrogen-bond symmetrization (Benoit & Marx, 2005; Morrone et al., 2009; Lin et al., 2011). Proton tunneling implies that the time-averaged proton site coincides with the O-O midpoint, although at GPa this is still not an equilibrium position, neither in the classical nor in the quantum frames. However, the average structure fully coincides with the HC- phase. Any structural measurement that is recorded on larger time lapses than the proton tunneling time would therefore provides such a picture. Nevertheless, only at larger pressures (about 20 GPa when NQE are included or around 25 GPa in classical simulations, at room temperature) the symmetric hydrogen bond corresponds to a maximum of the proton distribution. We propose that the measurement of the Debye-Waller factors through quasi-elastic neutron scattering (Teixeira J. & Luzar A., 1999) could distinguish between the two regimes.

As in other order-disorder phase transitions, the analysis of the O-H modes is crucial to point out the transition. From the trends of H-related modes with pressure, which we extracted from the dynamical simulations including NQE as well as anharmonic effects, the transition should happen at pressures as low as 10 GPa, where the O-H stretching modes considerably weaken and eventually fade out. Those trends are fully consistent with the analysis of the structural observables (average positions and pair-correlation functions) but they yield a much clearer signature of the transition as phonon modes are directly linked with the change of character of the proton dynamics. We suggest that high-resolution vibrational spectroscopy could complement the present knowledge and contribute to better identify the precise transition pressure. In order to point out the role of nuclear quantum effects, the experiments could be ideally conducted on both AlOOD and AlOOH.

Our results might also be relevant for other hydrous minerals, such as guyanaite (Jahn et al., 2012) or dense hydrous magnesium silicates (Ohira et al., 2014; Tsuchiya, 2013), where nuclear quantum effects, such as tunneling, could be at work.

We thank E. Balan for useful discussions and a critical reading of the manuscript. Y.B. acknowledges financial support from the Conseil Régional d’Île-de-France (DIM Oxymore). This work was granted access to the HPC resources of CINES under the allocation 2016096719 made by GENCI.

Figure 1: Schematic representation of the atomic motion of -AlOOH, for the imaginary frequency O-H stretching mode in the (HC) structure at GPa. Color code: H in light blue, O in red, Al in violet. The crystallographic conventions used here are the same as in Reference (Tsuchiya et al., 2008).
Figure 2: Equation of state for -AlOOH. The points computed through our molecular dynamics simulations at K with the Quantum Thermal Bath are in red. The fitted curves to the Vinet equation(Vinet et al., 1987) in the low-pressure (LP) and the high-pressure (HP) structures are drawn in blue and green, respectively; the corresponding values for the equilibrium volumes, the bulk moduli and their derivatives are: , , GPa, GPa and . Also shown are the values computed by Tsuchiya and coworkers at K in the and the structures.(Tsuchiya et al., 2002)
Figure 3: Average oxygen-oxygen distance computed from standard Langevin molecular dynamics and Quantum Thermal Bath simulations. The vertical bars indicate the widths of the probability distributions of O-O. The black symbols indicate experimental results.(Komatsu et al., 2006; Sano-Furukawa et al., 2008; Kuribayashi et al., 2014)
Figure 4: Oxygen-hydrogen pair correlation function for different pressures between and GPa, obtained from Quantum Thermal Bath (left) and standard Langevin (right) molecular dynamics simulations at ambient temperature. The two peaks indicate two different distances: the ionocovalent O-H bond length at approximately Å  and the HO hydrogen bond around Å. The two arrows in the left panel indicate neutron scattering measurements on -AlOOD at GPa.(Sano-Furukawa et al., 2008)
Figure 5: Proton transfer coordinate distribution , as a function of pressure in -AlOOH. The distributions are obtained through the time-averaged atomic position (blue and magenta) and the instanteneous atomic positions (red and green). Results are shown for both the Quantum Thermal Bath (Qu.) simulation, and classical Langevin (Cl.) simulations.
Figure 6: Pressure dependence of the hydrogen vibration modes frequencies in -AlOOH, obtained from Quantum Thermal Bath molecular dynamics simulations at ambient temperature. Upon compression, the O-H stretching modes soften drastically (red and green curves), while the O-H bending modes frequencies (in the plane in blue, and in the direction in purple) slightly increase. The black dots indicate Raman measurements on -AlOOH.(Ohtani et al., 2001) The vertical bars indicate the width of the peaks obtained via the Fourier Transform.


  1. Aoki, K., Yamawaki, H., Sakashita, M. & Fujihisa, H. (1996): Infrared absorption study of the hydrogen-bond symmetrization in ice to 110 GPa, Phys. Rev. B 54, 15673-15677
  2. Baroni, S., de Gironcoli, S., Dal Corso, A. & Giannozzi, P. (2001): Phonons and related crystal properties from density-functional perturbation theory, Rev. Modern Phys. 73, 515-562
  3. Basire, M., Borgis, D. & Vuilleumier, R. (2013): Computing Wigner distributions and time correlation functions using the quantum thermal bath method: Application to proton transfer spectroscopy, Phys. Chem. Chem. Phys. 15, 12591-12601
  4. Becke, A. D. (1988): Density-functional exchange-energy approximation with correct asymptotic behavior, Phys. Rev. A 38, 3098-3100
  5. Bedoya-Martinez, O. N., Barrat, J.-L. & Rodney, D. (2014): Computation of the thermal conductivity using methods based on classical and quantum molecular dynamics, Phys. Rev. B 89, 014303
  6. Benoit, M., Marx, D. & Parrinello, M. (1998): Tunnelling and zero-point motion in high-pressure ice, Nature 392, 258-261
  7. Benoit M. & Marx, D. (2005): The shapes of protons in hydrogen bonds depend on the bond length, Chem. Phys. Phys. Chem. 6, 1738-1741
  8. Berne, B. J.& Thirumalai, D. (1986): On the simulation of quantum systems: path integral methods, Ann. Rev. Phys. Chem. 37, 401-424
  9. Brieuc, F., Bronstein, Y., Dammak, H., Depondt, P., Finocchi, F. & Hayoun, M. (2016): Zero-point Energy Leakage in Quantum Thermal Bath Molecular Dynamics Simulations, J. Chem. Theory and Comput. 12, 5688-5697
  10. Bronstein, Y., Depondt, P., Finocchi, F. & Saitta A. M. (2014): Quantum-driven phase transition in ice described via an efficient Langevin approach, Phys. Rev. B 89, 214101
  11. Bronstein Y., Depondt, P., Bove, L. E., Gaal, R., Saitta, A. M. & Finocchi, F. (2016): Quantum versus classical protons in pure and salty ice under pressure, Phys. Rev. B 93, 024104
  12. Callen, H. B. & Welton, T. A. (1951): Irreversibility and generalized noise, Phys. Rev. 83, 34-40
  13. Calvo, F., Van-Oanh, N. T., Parneix, P. & Falvo, C. (2012): Vibrational spectra of polyatomic molecules assisted by quantum thermal baths, Phys. Chem. Chem. Phys. 14, 10503-10506
  14. Dakhnovskii, Y., Bursulaya, B. & Kim, H. J. (1995): Quantum tunneling in an anharmonic classical bath. Enhanced kinetic isotope effects in an Arrhenius region, J. Chem. Phys. 102, 7838
  15. Dammak, H., Chalopin, Y., Laroche, M., Hayoun, M. & Greffet, J. J. (2009): Quantum thermal bath for molecular dynamics simulation, Phys. Rev. Lett. 103, 190601
  16. Dammak, H., Antoshchenkova, E., Hayoun, M. & Finocchi, F. (2012): Isotope effects in lithium hydride and lithium deuteride crystals by molecular dynamics simulations, Journal of Physics: Condensed Matter 24, 435402
  17. Friedrich, A., Haussühl, E., Boehler, R., Morgenroth, W., Juarez-Arellano, E. A. & Winkler, B. (2007): Single-crystal structure refinement of diaspore at 50 GPa, Am. Mineral. 92, 1640-1644
  18. Giannozzi, P. et al. (2009): QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Cond. Matter 21, 395502
  19. Goncharov, A. F., Struzhkin, V. V., Somayazulu, M; S., Hemley, R. J. & Mao, H.-K. (1996): Compression of ice to 210 gigapascals: Infrared evidence for a symmetric hydrogen-bonded phase, Science 273, 218-220
  20. Goncharov, A. F., Struzhkin, V. V., Mao, H.-K. & Hemley, R. J. (1999): Raman spectroscopy of dense H O and the transition to symmetric hydrogen bonds, Phys. Rev. Lett. 83 1998-2001
  21. Hermann, H. & Mookerjee, M. (2016): High-pressure phase of brucite stable at earth’s mantle transition zone and lower mantle conditions, PNAS DOI: 10.1073/pnas/1611571113
  22. Hernandez-Rojas, J., Calvo, F. & Noya, E. G.: (2015) Applicability of Quantum Thermal Baths to Complex Many-Body Systems with Various Degrees of Anharmonicity, J. Chem. Theory Comput. 11, 861-870
  23. Hills, R. J. (1979): Crystal structure refinement and electron density distribution in diaspore, Phys. Chem. Minerals 5, 179-200
  24. Izaguirre, J. A., Catarello, D. P., Wozniak, J. M. & Skeel, R. D. (2001): Langevin stabilization of molecular dynamics, J. Chem. Phys. 114, 2090
  25. Jahn, S., Wunder, B., Koch-Müller, M., Tarrieu, L., Pöhle, M., Watenphul, A. & Taran, M. N. (2012): Pressure-induced hydrogen bond symmetrisation in guyanaite, β-CrOOH: evidence from spectroscopy and ab initio simulations, Eur. J. Mineral. 24, 839-850
  26. Johannsen, P. G. (1998): Vibrational states and optical transitions in hydrogen bonds, J. Physics: Cond. Matter 10, 2241-2260
  27. Komatsu, K., Kuribayashi, T., Sano, A., Ohtani, E. & Kudoh, Y. (2006): Redetermination of the high-pressure modification of AlOOH from single-crystal synchrotron data, Acta Cryst. E 62, i216-i218
  28. Kubo, R., Toda, M. & Hashitsume, N. (1978): Statistical Physics II: Nonequilibrium statistical mechanics, Springer
  29. Kuribayashi, T., Sano-Furukawa, A. & Nagase, T. (2014): Observation of pressure-induced phase transition of δ-AlOOH by using single-crystal synchrotron X-ray diffraction method, Phys. Chem. Minerals 41, 303-312
  30. Li, S., Ahuja, R. & Johansson, B. (2006): The elastic and optical properties of the high-pressure hydrous phase -AlOOH, Solid State Comm. 137, 101-105
  31. Lee, C., Yang, W. & Parr, R. G. (1988): Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density, Phys. Rev. B 37, 785-789
  32. Lin, L., Morrone, J. A. & Car, R. (2011): Correlated tunneling in hydrogen bonds, J. Stat. Phys. 145, 365-384
  33. Maradudin, A. A., Michel, T., McGurn A. R. & Mendez, E. R. (1990): Enhanced backscattering of light from a random grating, Annals of Physics (Leipzig) 203 255
  34. Marshall, W. & Lovesey, S. W. (1971): Theory of thermal neutron scattering: the use of neutrons for the investigation of condensed matter. Oxford University Press
  35. Marx, D. & Parrinello, M. (1996): Ab initio path integral molecular dynamics: Basic ideas, J. Chem. Phys. 104, 4077-4082
  36. Morrone, J. A., Lin, L. & Car, R. (2009): Tunneling and delocalization effects in hydrogen bonded systems: A study in position and momentum space, J. Chem. Phys. 130, 204511
  37. Nachimuthu, S., Gao, J. L. & Truhlar, D. G. (2012): A benchmark test suite for proton transfer energies and its use to test electronic structure model chemistries, Chemical Physics 400, 8-12
  38. Nielsen, O. H. & Martin, R. M. (1985): Quantum-mechanical theory of stress and force, Phys. Rev. B 32, 3792-3805
  39. Ohira, I., Ohtani, E., Sakai, T., Miyahara, M., Hirao, N., Ohishi Y. & Nishijima, M. (2014): Stability of a hydrous δ-phase, AlOOH–MgSiO 2 (OH) 2, and a mechanism for water transport into the base of lower mantle, Earth and Planetary Science Lett. 401, 12-17
  40. Ohtani, E., Litasov, K., Suzuki, A. & Kondo, T. (2001): Stability field of new hydrous phase, δ‐AlOOH, with implications for water transport into the deep mantle, Geophys. Research Lett. 28, 3991-3993
  41. Panero, W. R. & Stixrude, L. P. (2004): Hydrogen incorporation in stishovite at high pressure and symmetric hydrogen bonding in δ-AlOOH, Earth and Planetary Science Lett. 221, 421-431
  42. Perdew, J. P., Burke, K. & Ernzerhof, M. (1996): Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865
  43. Sano, A., E. Ohtani, Kubo, T. & Funakoshi, K. (2004): In situ X-ray observation of decomposition of hydrous aluminum silicate AlSiO OH and aluminum oxide hydroxide -AlOOH at high pressure and temperature, J. Phys. Chem. Solids 65, 1547-1554
  44. Sano, A., Ohtani, E., Kondo, T., Hirao, N., Sakai, T., Sata, N., Ohishi, Y. & Kikegawa, T. (2008): Aluminous hydrous mineral -AlOOH as a carrier of hydrogen into the core-mantle boundary, Geophys. Research Lett. 35, L03303
  45. Sano-Furukawa, A., Komatsu, K., Vanpeteghem, C. B. & Ohtani, E. (2008): Neutron diffraction study of -AlOOD at high pressure and its implication for symmetrization of the hydrogen bond, Am. Mineral. 93, 1558-1567
  46. Sano-Furukawa, A., Kagi, H., Nagai, T., Nakano, S., Ushijima, D., Iizuka, R., Ohtani, E. & Yagi, T. (2009): Change in compressibility of δ-AlOOH and δ-AlOOD at high pressure: a study of isotope effect and hydrogen-bond symmetrization, Am. Mineral. 94, 1255-1261
  47. Savin, A. V., Kosevich, Y. A. & Cantarero, A. (2012): Semiquantum molecular dynamics simulation of thermal properties and heat transport in low-dimensional nanostructures, Phys. Rev. B 86, 064305
  48. Struzhkin,V. V., Goncharov, A. F., Hemley, R. J. & Mao, H.-K. (1997): Cascading Fermi resonances and the soft mode in dense ice, Phys. Rev. Lett. 78 4446-4449
  49. Sugimura, E., Iiataka, T., Hirose, K., Kawamura, K., Sata, N. & Ohishi, Y. (2008): Compression of HO ice and implications for hydrogen-bond symmetrization: Synchrotron x-ray diffraction measurements and density-functional calculations, Phys. Rev. B 77 214103
  50. Suzuki, A., E. Ohtani, E. & Kamada, T. (2000): A new hydrous phase -AlOOH synthesized at 21 GPa and 1000 C, Phys. Chem. of Minerals 27,(2008) 689-693
  51. Teixeira J. & Luzar A. (1999): Physics of liquid water in Hydration Processes in Biology: Theoretical and Experimental Approaches, M.-C. Bellisent-Funel ed., IOS Press, NATO Science Studies, 35-65
  52. Tsuchiya, J., Tsuchiya, T., Tsuneyuki, T. & Yamanaka, T. (2002): First principles calculation of a high-pressure hydrous phase, -AlOOH, Geophys. Research Lett. 29, 1909
  53. Tsuchiya, J., Tsuchiya; T. & Wentzcovitch, R. M. (2008): Vibrational properties of -AlOOH under pressure, Am. Mineral. 93, 477-482
  54. Tsuchiya, J. & Tsuchiya, T. (2009): Elastic properties of -AlOOH under pressure: First principles investigation, Phys. Earth and Planetary Interiors 174, 122-127
  55. Tsuchiya, J. (2013): First principles prediction of a new high‐pressure phase of dense hydrous magnesium silicates in the lower mantle, Geophys. Research Lett. 40, 4570-4573
  56. Umemoto, K., Sugimura, E., de Gironcoli, S., Nakajima, Y., Hirose, K., Ohishi, Y. & Wentzcovitch, R. M. (2015): Nature of the Volume Isotope Effect in Ice, Phys. Rev. Lett. 115, 173005
  57. Vanpeteghem, C. B., Ohtani, E. & Kondo, K. (2002): Equation of state of the hydrous phase ‐AlOOH at room temperature up to 22.5 GPa, Geophys. Research Lett. 29, 1119-1121
  58. Vinet, P., Smith, J. R., Ferrante & Rose, J. H. (1987): Temperature effects on the universal equation of state of solids, Phys. Rev. B 35, 1945
  59. Winkler, B., Hytha, M., Pickard, P., Milman, V., Warren, M. & Segall, M. (2001): Theoretical investigation of bonding in diaspore, Eur. J. Mineral. 13, 343-349
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description